Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 5594)
+++ /issm/trunk/src/c/Makefile.am	(revision 5595)
@@ -241,4 +241,5 @@
 					./shared/Numerics/GaussPoints.cpp\
 					./shared/Numerics/BrentSearch.cpp\
+					./shared/Numerics/OptimalSearch.cpp\
 					./shared/Numerics/OptFunc.cpp\
 					./shared/Numerics/extrema.cpp\
@@ -791,4 +792,5 @@
 					./shared/Numerics/norm.cpp\
 					./shared/Numerics/BrentSearch.cpp\
+					./shared/Numerics/OptimalSearch.cpp\
 					./shared/Numerics/OptFunc.cpp\
 					./shared/Numerics/extrema.cpp\
Index: /issm/trunk/src/c/shared/Numerics/OptimalSearch.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/OptimalSearch.cpp	(revision 5595)
+++ /issm/trunk/src/c/shared/Numerics/OptimalSearch.cpp	(revision 5595)
@@ -0,0 +1,93 @@
+/*!\file:  OptimalSearch.cpp
+ * \brief optimization algorithm
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./numerics.h"
+#include "../../objects/objects.h"
+#include "../../include/include.h"
+#include "../../shared/shared.h"
+#include <float.h>
+
+void OptimalSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs){
+
+	/* This routine is optimizing a given function*/
+
+	/*function values: */
+	double fx1,fx2,fxbest;
+	double x1,x2,xmin,xbest;
+	double distance;
+
+	/*tolerances: */
+	double tol1,tol2,seps,tolerance;
+	int    maxiter;
+
+	/*counters: */
+	int  iter=0;
+	bool loop;
+
+	/*Recover parameters:*/
+	xmin     =optpars->xmin;
+	x1       =optpars->xmin;
+	x2       =optpars->xmax;
+	tolerance=optpars->tolerance;
+	maxiter  =optpars->maxiter;
+	
+	//get the value of the function at the first boundary
+	fx1= (*f)(x1,optargs);
+	if isnan(fx1) ISSMERROR("Function evaluation returned NaN");
+	_printf_("\n        Iteration         x           f(x)       Tolerance\n\n");
+	_printf_("        %s    %12.6g  %12.6g  %s","   N/A",x1,fx1,"         N/A\n");
+
+	//update tolerances
+	seps=sqrt(DBL_EPSILON);
+	tol1=seps*sqrt(pow(x1,2))+tolerance/3.0;
+	tol2=2.0*tol1;
+
+	loop=true;
+	while(loop){
+
+		/*get f(x2)*/
+		iter++;
+		fx2 = (*f)(x2,optargs);
+		if isnan(fx2) ISSMERROR("Function evaluation returned NaN");
+		_printf_("         %5i    %12.6g  %12.6g  %12.6g\n",iter,x2,fx2,pow(pow(x2-x1,2),0.5));
+
+		//Stop the optimization?
+		tol1=seps*pow(pow(x2,2),0.5)+tolerance/3.0;
+		tol2=2.0*tol1;
+		if (sqrt(pow(x2-x1,2)) < (tol2-0.5*(distance))){
+			_printf_("      %s%g\n","optimization terminated: the current x satisfies the termination criteria using 'tolx' of" ,tolerance);
+			loop=false;
+		}
+		else if (iter>=maxiter){
+			_printf_("      %s\n","exiting: Maximum number of iterations has been exceeded  - increase 'maxiter'\n");
+			loop=false;
+		}
+		else{
+			//continue
+			loop=true;
+		}
+
+		/*Compute x2*/
+		if(fx2<fx1){
+			xbest=x2; fxbest=fx2;
+			x1=x2;
+			x2=xmin+1.1*(x2-xmin); 
+			fx1=fx2;
+		}
+		else{
+			xbest=x1; fxbest=fx1;
+			x2=xmin+0.5*(x2-xmin);
+		}
+	}//end while
+
+	/*Assign output pointers: */
+	*psearch_scalar=xbest;
+	*pJ=fxbest;
+}
Index: /issm/trunk/src/c/shared/Numerics/numerics.h
===================================================================
--- /issm/trunk/src/c/shared/Numerics/numerics.h	(revision 5594)
+++ /issm/trunk/src/c/shared/Numerics/numerics.h	(revision 5595)
@@ -18,4 +18,5 @@
 double OptFunc(double scalar, OptArgs* optargs);
 void   BrentSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs);
+void   OptimalSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs);
 void   cross(double* result,double* vector1,double* vector2);
 double norm(double* vector);