COMMIT : Method 7 BFGS with MT Added

header/optimc.h

@@ -104,7 +104,7 @@ struct custom_jacobian_set{
 
 void setnlsTOL(nls_object obj,double gtol,double ftol,double xtol);
 
-void summary(opt_object obj);
+void optsummary(opt_object obj);
 
 void setMaxIter(opt_object obj,int MaxIter);
 

src/conjgrad.c

@@ -180,7 +180,7 @@ int cgpr_mt(custom_function *funcpt, custom_gradient *funcgrad, double *xc, int
 
 	gfdcode = grad_fd(funcpt,funcgrad, xi, N, dx,eps2, jac);
 	if (gfdcode == 15) {
-		return 15;
+		rcode = 15;
 	}
 	for (i = 0; i < N; ++i) {
 		pk[i] = -jac[i];
@@ -191,7 +191,11 @@ int cgpr_mt(custom_function *funcpt, custom_gradient *funcgrad, double *xc, int
 
 	if (fxf >= DBL_MAX || fxf <= -DBL_MAX) {
 		printf("Program Exiting as the function value exceeds the maximum double value");
-		return 15;
+		rcode = 15;
+	}
+	if (fxf != fxf) {
+		printf("Program Exiting as the function returns NaN");
+		rcode = 15;
 	}
 
 	if (restart < N) {

src/lnsrchmp.c

@@ -19,44 +19,6 @@ A Copy of the original Fortran routine is available at
 */
 
 
-double pmax(double a, double b) {
-	if (a > b) {
-		return a;
-	} else {
-		return b;
-	}
-}
-
-double pmin(double a, double b) {
-	if (a < b) {
-		return a;
-	}
-	else {
-		return b;
-	}
-}
-
-double signx(double x) {
-	double sgn;
-	if (x >= 0.) {
-		sgn = 1.0;
-	}
-	else {
-		sgn = -1.0;
-	}
-	return sgn;
-}
-
-double l2norm(double *vec, int N) {
-	double l2, sum;
-	int i;
-	sum = 0.;
-	for (i = 0; i < N; ++i) {
-		sum += vec[i] * vec[i];
-	}
-	l2 = sqrt(sum);
-	return l2;
-}
 
 int grad_fd(custom_function *funcpt, custom_gradient *funcgrad, double *x, int N, double *dx,
 		double eps2, double *f) {
@@ -113,12 +75,20 @@ int grad_calc2(custom_function *funcpt, double *x, int N, double *dx, double eps
 			printf("Program Exiting as the function value exceeds the maximum double value");
 			return 15;
 		}
+		if (fp != fp) {
+			printf("Program Exiting as the function returns NaN");
+			return 15;
+		}
 		x[j] = temp - stepsize;
 		fm = FUNCPT_EVAL(funcpt,x,N);
 		if (fm >= DBL_MAX || fm <= -DBL_MAX) {
 			printf("Program Exiting as the function value exceeds the maximum double value");
 			return 15;
 		}
+		if (fm != fm) {
+			printf("Program Exiting as the function returns NaN");
+			return 15;
+		}
 		f[j] = (fp - fm)/ (2 * stepsize);
 		x[j] = temp;
 	}
@@ -150,6 +120,12 @@ int grad_calc(custom_function *funcpt, double *x, int N, double *dx, double eps2
 		f[i] = (FUNCPT_EVAL(funcpt, xi, N) - FUNCPT_EVAL(funcpt, x, N)) / step;
 		if (f[i] >= DBL_MAX || f[i] <= -DBL_MAX) {
 			printf("Program Exiting as the function value exceeds the maximum double value");
+			free(xi);
+			return 15;
+		}
+		if (f[i] != f[i]) {
+			printf("Program Exiting as the function returns NaN");
+			free(xi);
 			return 15;
 		}
 		//xi[i] -= step;
@@ -607,13 +583,19 @@ int cvsrch(custom_function *funcpt, custom_gradient *funcgrad, double *x, double
 		*f = FUNCPT_EVAL(funcpt,x, N);
 		if (*f >= DBL_MAX || *f <= -DBL_MAX) {
 			printf("Program Exiting as the function value exceeds the maximum double value");
+			free(rcheck);
+			free(wa);
+			return 15;
+		}
+		if (*f != *f) {
+			printf("Program Exiting as the function returns NaN");
+			free(rcheck);
+			free(wa);
 			return 15;
 		}
 		grad_cd(funcpt,funcgrad, x, N, dx, eps2,g);
 		nfev++;
 
-		//printf("ITER %d stp %g \n", nfev,stp);
-		//mdisplay(x, 1, N);
 
 		dg = 0.0;
 		for (j = 0; j < N; ++j) {
@@ -650,6 +632,8 @@ int cvsrch(custom_function *funcpt, custom_gradient *funcgrad, double *x, double
 		if (stage1 == 1 && *f <= ftest1 && dg >= pmin(ftol, gtol)*dginit) {
 			stage1 = 0;
 		}
+
+
 		/*
 		A modified function is used to predict the step only if
         we have not obtained a step for which the modified

src/lnsrchmp.h

@@ -54,14 +54,6 @@ struct custom_jacobian_set{
 
 #define JACOBIAN_EVAL(F,x,M,N,jac) (*((F)->jacobian))(x,M,N,(jac),(F)->params)
 
-double pmax(double a, double b);
-
-double pmin(double a, double b);
-
-double signx(double x);
-
-double l2norm(double *vec, int N);
-
 int stopcheck_mt(double fx, int N, double *xc, double *xf, double *jac, double *dx, double fsval, double gtol, double stol, int retval);
 
 int stopcheck2_mt(double fx, int N, double fo, double *jac, double *dx, double eps,double stoptol, double functol, int retval);

src/matrix.c

@@ -28,6 +28,64 @@ double macheps() {
 	return macheps;
 }
 
+double pmax(double a, double b) {
+	if (a > b) {
+		return a;
+	}
+	else {
+		return b;
+	}
+}
+
+double pmin(double a, double b) {
+	if (a < b) {
+		return a;
+	}
+	else {
+		return b;
+	}
+}
+
+int imax(int a, int b) {
+	if (a > b) {
+		return a;
+	}
+	else {
+		return b;
+	}
+}
+
+int imin(int a, int b) {
+	if (a < b) {
+		return a;
+	}
+	else {
+		return b;
+	}
+}
+
+double signx(double x) {
+	double sgn;
+	if (x >= 0.) {
+		sgn = 1.0;
+	}
+	else {
+		sgn = -1.0;
+	}
+	return sgn;
+}
+
+double l2norm(double *vec, int N) {
+	double l2, sum;
+	int i;
+	sum = 0.;
+	for (i = 0; i < N; ++i) {
+		sum += vec[i] * vec[i];
+	}
+	l2 = sqrt(sum);
+	return l2;
+}
+
 int compare (const void* ind1, const void* ind2)
 {
 	if (*((vipair *)ind1)->a > *((vipair *)ind2)->a)

src/matrix.h

@@ -28,6 +28,18 @@ extern "C" {
 
 double macheps();
 
+double pmax(double a, double b);
+
+double pmin(double a, double b);
+
+int imax(int a, int b);
+
+int imin(int a, int b);
+
+double signx(double x);
+
+double l2norm(double *vec, int N);
+
 int compare (const void* ind1, const void* ind2);
 
 void sort1d(double* v,int N, int* pos);