Changeset 2028 for branches/newmole/source/newton_step.cpp

Timestamp:

05/06/08 16:35:36 (8 months ago)

Author:

rjrw

Message:

Further refactor newton_step to more purely numerical function.

Remove more dead or leftover code and data.

Files:

• branches/newmole/source/newton_step.cpp (modified) (11 diffs)

branches/newmole/source/newton_step.cpp

@@ -11 +11 @@
 #define ABSLIM  1e-12
 #define ERRLIM  1e-12
-#define SMALLABUND 1e-24
 #       ifdef MAT
 #               undef MAT
@@ -21 +20 @@
 enum {PRINTSOL = false};
 
-extern void funjac(valarray<double> &b2vec, double *ervals, double *amat, bool lgJac);
-
 /* mole_newton_step -- improve balance in chemical network along
  * descent direction, step limited to ensure improvement */
-void newton_step(valarray<double> &b2vec, int *nBad, realnum *error, long n)
+void newton_step(const valarray<double> &b0vec, valarray<double> &b2vec, realnum *error, const long n,
+                        void (*jacobn)(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac))
 {
         long int i, 
@@ -31 +29 @@
                 iworst;
 
-        int printsol = PRINTSOL;
-
         double 
                 etmp,
                 emax,
-                fracneg,
                 fstep,
                 flast,
@@ -45 +40 @@
         valarray<double> 
                 amat(n*n),
-                b0vec(n),
                 b1vec(n),
                 escale(n),
@@ -51 +45 @@
         int32 merror;
 
-        DEBUG_ENTRY( "mole_newton_step()" );
+        DEBUG_ENTRY( "newton_step()" );
 
         for( i=0; i < n; i++ )
         {
-                b0vec[i] = b2vec[i];
+                b2vec[i] = b0vec[i];
         }
 
@@ -61 +55 @@
         error1 = grad = 0.;
         fstep = flast = 1.;                     
-        *nBad = 0; 
 
         for (loop=0;loop<20;loop++) 
@@ -67 +60 @@
                 flast = fstep;
                 
-                funjac(b2vec,&ervals[0],&amat[0],loop==0);
+                jacobn(b2vec,&ervals[0],&amat[0],loop==0);
 
                 if (loop == 0) 
@@ -89 +82 @@
                                 iworst = i;
                         }
-                        
-                        if(printsol)
-                                fprintf(ioQQQ,"%15.7g\t",etmp);
                         error0 += etmp;
                 }
-
-                // fprintf(stdout,"Loop %ld error %15.8g fstep %15.8g; worst species %s abund %g  error %g\n",
-                //                              loop,error0,fstep,groupspecies[iworst]->label,groupspecies[iworst]->den,,emax);
 
                 if (loop == 0) 
@@ -107 +94 @@
                         error1 = error0;
                         grad = -2*error1;
-                        /* b1vec is now descent direction */
-
-                        //fprintf(ioQQQ,"Step factor %g bad species %s abund %g step %g\n",
-                        //                              fstep,iworst!=-1?groupspecies[iworst]->label:"None",
-                        //                              iworst!=-1?b0vec[iworst]:0.,
-                        //                              iworst!=-1?-b1vec[iworst]:0.); 
                 } else {
-                        //fprintf(ioQQQ,"Step factor %15.8g bad species %d %s abund %15.8g step %15.8g\n",
-                        //                              fstep,iworst,iworst!=-1?groupspecies[iworst]->label:"None",
-                        //                              iworst!=-1?b0vec[iworst]:0.,
-                        //                              iworst!=-1?-b1vec[iworst]:0.); 
-                        //fprintf(stdout,"Grads %15.8g %15.8g (%15.8g %15.8g)\n",grad,(error0-error1)/fstep,error0,error1);
                         if (error0 < (1-2e-4*fstep)*error1 || error0 < 1e-20)
                         {
@@ -124 +100 @@
                         }
                         fstep = -0.5*grad*fstep*fstep/SDIV(error0-error1-grad*fstep);
-                        //fprintf(stdout,"Errs1 loop %ld %g now %g fac %g was %g pred %g\n",loop,
-                        //                              error1,error0,fstep,flast,
-                        //                              error1+grad*fstep+(fstep*fstep/(flast*flast))*(error0-(error1+grad*flast)));
                         if (fstep > 0.9*flast || fstep < 0.)
                                 fstep = 0.5*flast;
@@ -133 +106 @@
                 }
 
-                //if (iworst != -1)
-                //      fprintf(stdout,"Loop %ld, worst error for %s: %g of %g, step %g of %g %g\n",                                                    
-                //                                      loop, groupspecies[iworst]->label,emax,error0,-b1vec[iworst]*fstep,b0vec[iworst],fstep);
-                
-                //fprintf(ioQQQ,"Loops %ld %g %g\n",loop,error0,(1-2e-4*fstep)*error1);
-
                 /*
                  * b1vec[] is descent direction
-                 * new density in b2vec[]
+                 * new densities in b2vec[]
                  */
-
+                
                 for( i=0; i < n; i++ )
                 {
                         b2vec[i] = b0vec[i]-b1vec[i]*fstep;                             
                 }
-                
-                /* check for negative populations */
-                *nBad = 0; 
-                fracneg = 0.;
-                iworst = -1;
-                for( i=0; i < n; i++ )
-                {
-                        if( b2vec[i] < 0. ) 
-                        {
-                                // fprintf(stdout,"Species %s negative at %g from %g\n",groupspecies[i]->label,b2vec[i],b0vec[i]);
-                                /* largest relative change in solution for neg soln */
-                                /* this can only occur for negative solutions since fracneg starts
-                                 * as zero */
-                                if(-b2vec[i] > fracneg*SDIV(b0vec[i]) && b0vec[i] > SMALLABUND) 
-                                {
-                                        fracneg = -b2vec[i]/SDIV(b0vec[i]);
-                                        iworst = i;
-                                }
-                                if (b0vec[i] > SMALLABUND || b2vec[i] < -SMALLABUND)
-                                        ++*nBad;
-                        }
-                }
-                                                                        
-                /* Expand species groups into mole.list[]->den -- either for
-                 * re-evaluation of network error, or for communication with
-                 * outside code */
-                
         }
-        //fprintf(ioQQQ,"Final fstep %g\n",fstep);
 
         *error = (realnum) MIN2(error0,1e30);