--- imach/src/imach.c	2014/09/02 09:24:05	1.160
+++ imach/src/imach.c	2014/09/15 20:41:41	1.161
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.160 2014/09/02 09:24:05 brouard Exp $
+/* $Id: imach.c,v 1.161 2014/09/15 20:41:41 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.161  2014/09/15 20:41:41  brouard
+  Summary: Problem with macro SQR on Intel compiler
+
   Revision 1.160  2014/09/02 09:24:05  brouard
   *** empty log message ***
 
@@ -547,11 +550,11 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.160 2014/09/02 09:24:05 brouard Exp $ */
+/* $Id: imach.c,v 1.161 2014/09/15 20:41:41 brouard Exp $ */
 /* $State: Exp $ */
 
 char version[]="Imach version 0.98nX, August 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";
-char fullversion[]="$Revision: 1.160 $ $Date: 2014/09/02 09:24:05 $"; 
+char fullversion[]="$Revision: 1.161 $ $Date: 2014/09/15 20:41:41 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1400,23 +1403,43 @@ void powell(double p[], double **xi, int
       return; 
     } 
     if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); 
-    for (j=1;j<=n;j++) { 
+    for (j=1;j<=n;j++) { /* Computes an extrapolated point */
       ptt[j]=2.0*p[j]-pt[j]; 
       xit[j]=p[j]-pt[j]; 
       pt[j]=p[j]; 
     } 
     fptt=(*func)(ptt); 
-    if (fptt < fp) { 
-      t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); 
-      if (t < 0.0) { 
-	linmin(p,xit,n,fret,func); 
+    if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
+      /* x1 f1=fp x2 f2=*fret x3 f3=fptt, xm fm */
+      /* From x1 (P0) distance of x2 is at h and x3 is 2h */
+      /* Let f"(x2) be the 2nd derivative equal everywhere. Then the parabolic through (x1,f1), (x2,f2) and (x3,f3)
+	 will reach at f3 = fm + h^2/2 f''m  ; f" = (f1 -2f2 +f3 ) / h**2 */
+      /* f1-f3 = delta(2h) = 2 h**2 f'' = 2(f1- 2f2 +f3) */
+      /* Thus we compare delta(2h) with observed f1-f3 */
+      /* or best gain on one ancient line 'del' with total gain f1-f2 = f1 - f2 - 'del' with 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= t- del*SQR(fp-fptt);
+      printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);
+      fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);
+#ifdef DEBUG
+      printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
+	     (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
+      fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
+	     (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
+      printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
+      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
+      if (t < 0.0) { /* Then we use it for last direction */
+	linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/
 	for (j=1;j<=n;j++) { 
-	  xi[j][ibig]=xi[j][n]; 
-	  xi[j][n]=xit[j]; 
+	  xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */
+	  xi[j][n]=xit[j];      /* and nth direction by the extrapolated */
 	}
+	printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+	fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction :\n",n,ibig);
+
 #ifdef DEBUG
-	printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
-	fprintf(ficlog,"Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
 	for(j=1;j<=n;j++){
 	  printf(" %.12e",xit[j]);
 	  fprintf(ficlog," %.12e",xit[j]);
@@ -6569,7 +6592,8 @@ Interval (in months) between two waves:
 
 	for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
 	  oldm=oldms;savm=savms; /* Segmentation fault */
-	  varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart);
+	  cptcod= 0; /* To be deleted */
+	  varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
 	  fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
 	  if(vpopbased==1)
 	    fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);