31 realnum *eqerror,
realnum *error,
const long n,
double *rlimit,
double *rmax,
32 valarray<double> &escale,
34 const valarray<double> &b2vec,
35 double *
const ervals,
double *
const amat,
36 const bool lgJac,
bool *lgConserve))
62 for( i=0; i < n; i++ )
69 pred = error0 = error2 = erroreq = grad = f1 = f2 = 0.;
72 const int LOOPMAX = 40;
73 for (loop=0;loop<LOOPMAX;loop++)
80 for( i=0; i < n; i++ )
97 *rlimit = 1e-19 * (*rmax);
99 else if (*rlimit > *rmax)
104 for( i=0; i < n; i++ )
118 for( i=0; i < n; i++ )
120 double etmp = ervals[i]/(
SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
121 erroreq += etmp*etmp;
125 for (i=0; i < n; i++)
127 ervals1[i] = ervals[i];
129 ervals1[i] -= (*rlimit)*(b2vec[i]-b0vec[i]);
134 double emax0 = 0.0, emax1 = 0.0;
135 for( i=0; i < n; i++ )
140 double etmp = ervals1[i]/(
SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
141 double etmp0 = ervals0[i]/(
SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
144 if (fabs(etmp-etmp0) > fabs(emax1-emax0) )
156 for( i=0; i < n; i++ )
158 ervals0[i] = ervals1[i];
159 b1vec[i] = ervals1[i];
164 *eqerror = *error = 1e30f;
175 if (error1 < (1-2e-4*f1)*error0 || error1 < 1e-20)
183 f1 *= -0.5*f1*grad/(error1-error0-f1*grad);
184 pred = error0+0.5*grad*f1;
190 fprintf(
ioQQQ,
"Newt %ld stp %11.4g err %11.4g erreq %11.4g rlimit %11.4g grd %11.4g fdg %11.4g e0 %11.4g de %11.4g pred %11.4g\n",
191 loop,f1,error1,erroreq,*rlimit,grad,(error1-error0)/f1,error0,error1-error0,pred
194 fprintf(
ioQQQ,
"Maxi %ld %s emax from %11.4g of %11.4g -> %11.4g of %11.4g, diff %11.4g\n",
195 maxi,
groupspecies[maxi]->label.c_str(),emax0,error0, emax1,error1,emax1-emax0);
199 double a = (error1-error0)/(f1*f1) - (error2-error0)/(f2*f2);
200 a += grad * (1./f2 - 1./f1);
202 double b = -f2*(error1-error0)/(f1*f1) + f1*(error2-error0)/(f2*f2);
203 b += grad * (f2/f1 - f1/f2);
212 f1 = 1.-3.*(a/b)*(grad/b)*(ft-f2);
214 f1 = b/(3.*a)*(sqrt(f1)-1.);
223 pred = error0+f1*(grad+f1/(ft-f2)*(b+a*f1));
227 if (f1 > 0.5*f2 || f1 < 0.)
229 else if (f1 < 0.03*f2)
242 for( i=0; i < n; i++ )
244 b2vec[i] = b0vec[i]-b1vec[i]*f1;
251 for( i=0; i < n; i++ )
258 if (0 && LOOPMAX == loop)
260 double rvmax = 0., rval;
262 for( i=0; i < n; ++i )
265 rval = b1vec[i]/b0vec[i];
268 fprintf(
ioQQQ,
"%7s %11.4g %11.4g\n",
270 if (fabs(rval) > rvmax)
280 for( j=0; j < n; j++ )
286 for( i=0; i < n; i++ )
288 for( j=0; j < n; j++ )
292 double db = 1e-3*fabs(b0vec[i])+1e-9;
294 db = b1vec[i]-b0vec[i];
296 for( j=0; j < n; j++ )
299 double e2 = (escale[j]-ervals1[j])/db;
300 if (fabs(e1-
e2) > 0.01*fabs(e1+
e2))
301 fprintf(
ioQQQ,
"%7s %7s %11.4g %11.4g %11.4g\n",
365 const valarray<double> &a,
const valarray<double> &b)
367 fprintf(
ioQQQ,
" CO_solve getrf_wrapper error %ld",(
long int) merror );
368 if( merror > 0 && merror <= n )
370 fprintf(
ioQQQ,
" - problem with species %s\n\n",
groupspecies[merror-1]->label.c_str() );
371 fprintf(
ioQQQ,
"index \t Row A(i,%li)\t Col A(%li,j) \t B \t Species\n", merror, merror );
372 for(
long index=1; index<=n; index++ )
374 fprintf(
ioQQQ,
"%li\t%+.4e\t%+.4e\t%+.4e\t%s\n", index,
375 a[(merror-1)*n + index - 1],
376 a[(index -1)*n + merror- 1],
384 fprintf(
ioQQQ,
"\n" );
391 valarray<int32> ipiv(n);
395 valarray<double> lufac(n*n),oldb(n),err(n);
397 ASSERT(a.size() ==
size_t(n*n));
398 ASSERT(b.size() ==
size_t(n));
413 if (error_print != NULL)
414 error_print(n,merror,a,b);
416 fprintf(
ioQQQ,
"Singular matrix in solve_system\n");
424 fprintf(
ioQQQ,
" solve_system: dgetrs finds singular or ill-conditioned matrix\n" );
428 for (k=0;k<nrefine;k++)
438 err[i] -= a[i+j*n]*b[j];
446 double maxb=0., maxe=0.;
451 if (fabs(err[i])>maxe)
456 fprintf(
ioQQQ,
"Max error %g, Condition %g\n",maxe,1e16*maxe/maxb);
bool newton_step(GroupMap &MoleMap, const valarray< double > &b0vec, valarray< double > &b2vec, realnum *eqerror, realnum *error, const long n, double *rlimit, double *rmax, valarray< double > &escale, void(*jacobn)(GroupMap &MoleMap, const valarray< double > &b2vec, double *const ervals, double *const amat, const bool lgJac, bool *lgConserve))
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)