Summary: Problem with macro SQR on Intel compiler

diff --git a/src/imach.c b/src/imach.c
index fb56f5f1d22885e4e6fa42418b4a5704f85c457a..72055c0c305fbbcc1d0ff2af29cf6743b8a07a2e 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.160  2014/09/02 09:24:05  brouard
+  *** empty log message ***
+
   Revision 1.159  2014/09/01 10:34:10  brouard
   Summary: WIN32
   Author: Brouard
@@ -1397,23 +1400,43 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
       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]);
@@ -6566,7 +6589,8 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
 
        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);