Summary: Trying to keep directest which seems simpler and more general

Author: Nicolas Brouard

diff --git a/src/imach.c b/src/imach.c
index 235b7f1f210112661971b1147c3d07a2e547b8bf..97868a63cee71ccb713998dc64832700d8c7c51b 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,11 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.181  2015/02/11 23:22:24  brouard
+  Summary: Comments on Powell added
+
+  Author:
+
   Revision 1.180  2015/02/11 17:33:45  brouard
   Summary: Finishing move from main to function (hpijx and prevalence_limit)
 
@@ -1491,7 +1496,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
       fprintf(ficlog,"\n");
 #endif
     } /* end i */
-    if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
+    if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */
 #ifdef DEBUG
       int k[2],l;
       k[0]=1;
@@ -1538,9 +1543,9 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
       /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del */
       /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
 
-      t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del);
+      t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line */
       t= t- del*SQR(fp-fptt);
-      directest = SQR(fp-2.0*(*fret)+fptt) - 2.0 * del; /* If del was big enough we change it for a new direction */
+      directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */
 #ifdef DEBUG
       printf("t1= %.12lf, t2= %.12lf, t=%.12lf  directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
       fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
@@ -1552,6 +1557,12 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
       fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
 #endif
 #ifdef POWELLDIRECT
+      if (directest*t < 0.0) { /* Contradiction between both tests */
+      printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt);
+      printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
+      fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt);
+      fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
+    } 
       if (directest < 0.0) { /* Then we use it for new direction */
 #else
       if (t < 0.0) { /* Then we use it for new direction */