--- imach/src/imach.c	2015/02/11 23:22:24	1.181
+++ imach/src/imach.c	2015/02/12 08:19:57	1.182
@@ -1,6 +1,10 @@
-/* $Id: imach.c,v 1.181 2015/02/11 23:22:24 brouard Exp $
+/* $Id: imach.c,v 1.182 2015/02/12 08:19:57 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.182  2015/02/12 08:19:57  brouard
+  Summary: Trying to keep directest which seems simpler and more general
+  Author: Nicolas Brouard
+
   Revision 1.181  2015/02/11 23:22:24  brouard
   Summary: Comments on Powell added
 
@@ -644,11 +648,11 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.181 2015/02/11 23:22:24 brouard Exp $ */
+/* $Id: imach.c,v 1.182 2015/02/12 08:19:57 brouard Exp $ */
 /* $State: Exp $ */
 
 char version[]="Imach version 0.98p, Février 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
-char fullversion[]="$Revision: 1.181 $ $Date: 2015/02/11 23:22:24 $"; 
+char fullversion[]="$Revision: 1.182 $ $Date: 2015/02/12 08:19:57 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1496,7 +1500,7 @@ void powell(double p[], double **xi, int
       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;
@@ -1543,9 +1547,9 @@ void powell(double p[], double **xi, int
       /* 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);
@@ -1557,6 +1561,12 @@ void powell(double p[], double **xi, int
       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 */