COPASI API  4.16.103
Classes | Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
CInternalSolver Class Reference

#include <CInternalSolver.h>

Inheritance diagram for CInternalSolver:
Inheritance graph
[legend]
Collaboration diagram for CInternalSolver:
Collaboration graph
[legend]

Classes

struct  State
 

Public Member Functions

void enablePrint (const bool &print=true)
 
void resetState ()
 
void saveState ()
 
void setOstream (std::ostream &os)
 
 ~CInternalSolver ()
 

Protected Member Functions

 CInternalSolver ()
 
C_INT dintdy_ (double *t, const C_INT *k, double *yh, C_INT *nyh, double *dky, C_INT *iflag)
 
C_INT dprja_ (C_INT *neq, double *y, double *yh, C_INT *nyh, double *ewt, double *ftem, double *savf, double *wm, C_INT *iwm, evalF f, evalJ jac)
 
C_INT drchek_ (const C_INT *job, evalG g, C_INT *neq, double *y, double *yh, C_INT *nyh, double *g0, double *g1, double *gx, C_INT *jroot, C_INT *irt)
 
C_INT droots_ (C_INT *ng, double *hmin, C_INT *jflag, double *x0, double *x1, double *g0, double *g1, double *gx, double *x, C_INT *jroot)
 
C_INT dsolsy_ (double *wm, C_INT *iwm, double *x, double *tem)
 
C_INT dstoda_ (C_INT *neq, double *y, double *yh, C_INT *nyh, double *yh1, double *ewt, double *savf, double *acor, double *wm, C_INT *iwm, evalF f, evalJ jac, PJAC *pjac, SLVS *slvs)
 

Protected Attributes

dls001 mdls001_
 
dlsa01 mdlsa01_
 
dlsr01 mdlsr01_
 
State mState
 
Cxerrwd mxerrwd
 

Detailed Description

Definition at line 25 of file CInternalSolver.h.

Constructor & Destructor Documentation

CInternalSolver::CInternalSolver ( )
protected

Definition at line 12 of file CInternalSolver.cpp.

12  :
13  mxerrwd(true),
14  mdls001_(),
15  mdlsa01_(),
16  mdlsr01_()
17 {}
CInternalSolver::~CInternalSolver ( )

Definition at line 19 of file CInternalSolver.cpp.

20 {}

Member Function Documentation

C_INT CInternalSolver::dintdy_ ( double *  t,
const C_INT k,
double *  yh,
C_INT nyh,
double *  dky,
C_INT iflag 
)
protected

Definition at line 53 of file dintdy.cpp.

References c__0, c__1, c__2, c__30, c__51, c__52, c__60, c_b34, C_INT, d_sign(), dls001_1, if(), mxerrwd, and pow_di.

Referenced by drchek_(), CLSODA::operator()(), and CLSODAR::operator()().

55 {
56  /* System generated locals */
57  C_INT yh_dim1, yh_offset, i__1, i__2;
58  double d__1;
59 
60  /* Local variables */
61  double c__;
62  C_INT i__, j;
63  double r__, s;
64  C_INT ic, jb, jj;
65  double tp;
66  C_INT jb2, jj1, jp1;
67  std::string msg;
68 
69  /* ***BEGIN PROLOGUE DINTDY */
70  /* ***SUBSIDIARY */
71  /* ***PURPOSE Interpolate solution derivatives. */
72  /* ***TYPE DOUBLE PRECISION (SINTDY-S, DINTDY-D) */
73  /* ***AUTHOR Hindmarsh, Alan C., (LLNL) */
74  /* ***DESCRIPTION */
75 
76  /* DINTDY computes interpolated values of the K-th derivative of the */
77  /* dependent variable vector y, and stores it in DKY. This routine */
78  /* is called within the package with K = 0 and T = TOUT, but may */
79  /* also be called by the user for any K up to the current order. */
80  /* (See detailed instructions in the usage documentation.) */
81 
82  /* The computed values in DKY are gotten by interpolation using the */
83  /* Nordsieck history array YH. This array corresponds uniquely to a */
84  /* vector-valued polynomial of degree NQCUR or less, and DKY is set */
85  /* to the K-th derivative of this polynomial at T. */
86  /* The formula for DKY is: */
87  /* q */
88  /* DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1) */
89  /* j=K */
90  /* where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR. */
91  /* The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are */
92  /* communicated by COMMON. The above sum is done in reverse order. */
93  /* IFLAG is returned negative if either K or T is out of bounds. */
94 
95  /* ***SEE ALSO DLSODE */
96  /* ***ROUTINES CALLED XERRWD */
97  /* ***COMMON BLOCKS DLS001 */
98  /* ***REVISION HISTORY (YYMMDD) */
99  /* 791129 DATE WRITTEN */
100  /* 890501 Modified prologue to SLATEC/LDOC format. (FNF) */
101  /* 890503 Minor cosmetic changes. (FNF) */
102  /* 930809 Renamed to allow single/double precision versions. (ACH) */
103  /* 010418 Reduced size of Common block /DLS001/. (ACH) */
104  /* 031105 Restored 'own' variables to Common block /DLS001/, to */
105  /* enable interrupt/restart feature. (ACH) */
106  /* 050427 Corrected roundoff decrement in TP. (ACH) */
107  /* ***END PROLOGUE DINTDY */
108  /* **End */
109 
110  /* ***FIRST EXECUTABLE STATEMENT DINTDY */
111  /* Parameter adjustments */
112  yh_dim1 = *nyh;
113  yh_offset = 1 + yh_dim1;
114  yh -= yh_offset;
115  --dky;
116 
117  /* Function Body */
118  *iflag = 0;
119 
120  if (*k < 0 || *k > dls001_1.nq)
121  {
122  goto L80;
123  }
124 
125  d__1 = fabs(dls001_1.tn) + fabs(dls001_1.hu);
126  tp = dls001_1.tn - dls001_1.hu - dls001_1.uround * 100. * d_sign(d__1, dls001_1.hu);
127 
128  if ((*t - tp) * (*t - dls001_1.tn) > 0.)
129  {
130  goto L90;
131  }
132 
133  s = (*t - dls001_1.tn) / dls001_1.h__;
134  ic = 1;
135 
136  if (*k == 0)
137  {
138  goto L15;
139  }
140 
141  jj1 = dls001_1.l - *k;
142  i__1 = dls001_1.nq;
143 
144  for (jj = jj1; jj <= i__1; ++jj)
145  {
146  /* L10: */
147  ic *= jj;
148  }
149 
150 L15:
151  c__ = (double) ic;
152  i__1 = dls001_1.n;
153 
154  for (i__ = 1; i__ <= i__1; ++i__)
155  {
156  /* L20: */
157  dky[i__] = c__ * yh[i__ + dls001_1.l * yh_dim1];
158  }
159 
160  if (*k == dls001_1.nq)
161  {
162  goto L55;
163  }
164 
165  jb2 = dls001_1.nq - *k;
166  i__1 = jb2;
167 
168  for (jb = 1; jb <= i__1; ++jb)
169  {
170  j = dls001_1.nq - jb;
171  jp1 = j + 1;
172  ic = 1;
173 
174  if (*k == 0)
175  {
176  goto L35;
177  }
178 
179  jj1 = jp1 - *k;
180  i__2 = j;
181 
182  for (jj = jj1; jj <= i__2; ++jj)
183  {
184  /* L30: */
185  ic *= jj;
186  }
187 
188 L35:
189  c__ = (double) ic;
190  i__2 = dls001_1.n;
191 
192  for (i__ = 1; i__ <= i__2; ++i__)
193  {
194  /* L40: */
195  dky[i__] = c__ * yh[i__ + jp1 * yh_dim1] + s * dky[i__];
196  }
197 
198  /* L50: */
199  }
200 
201  if (*k == 0)
202  {
203  return 0;
204  }
205 
206 L55:
207  i__1 = -(*k);
208  r__ = pow_di(&dls001_1.h__, &i__1);
209  i__1 = dls001_1.n;
210 
211  for (i__ = 1; i__ <= i__1; ++i__)
212  {
213  /* L60: */
214  dky[i__] = r__ * dky[i__];
215  }
216 
217  return 0;
218 
219 L80:
220  msg = "DINTDY- K (=I1) illegal ";
221  mxerrwd(msg, &c__30, &c__51, &c__0, &c__1, k, &c__0, &c__0, &c_b34, &
222  c_b34, (C_INT)80);
223  *iflag = -1;
224  return 0;
225 L90:
226  msg = "DINTDY- T (=R1) illegal ";
227  mxerrwd(msg, &c__30, &c__52, &c__0, &c__0, &c__0, &c__0, &c__1, t, &c_b34,
228  (C_INT)80);
229  msg = " T not in interval TCUR - HU (= R1) to TCUR (=R2) ";
230  mxerrwd(msg, &c__60, &c__52, &c__0, &c__0, &c__0, &c__0, &c__2, &tp, &
231  dls001_1.tn, (C_INT)80);
232  *iflag = -2;
233  return 0;
234  /* ----------------------- END OF SUBROUTINE DINTDY ---------------------- */
235 } /* dintdy_ */
#define C_INT
Definition: copasi.h:115
static C_INT c__0
Definition: dintdy.cpp:41
static C_INT c__30
Definition: dintdy.cpp:44
static C_INT c__60
Definition: dintdy.cpp:47
#define dls001_1
Definition: dintdy.cpp:29
#define pow_di(__x, __n)
Definition: dintdy.cpp:49
static C_INT c__2
Definition: dintdy.cpp:43
double d_sign(const double &a, const double &b)
Definition: CLSODA.cpp:2471
static C_INT c__52
Definition: dintdy.cpp:46
static C_INT c__51
Definition: dintdy.cpp:45
static double c_b34
Definition: dintdy.cpp:39
static C_INT c__1
Definition: dintdy.cpp:42
if(!yymsg) yymsg
C_INT CInternalSolver::dprja_ ( C_INT neq,
double *  y,
double *  yh,
C_INT nyh,
double *  ewt,
double *  ftem,
double *  savf,
double *  wm,
C_INT iwm,
evalF  f,
evalJ  jac 
)
protected

Definition at line 47 of file dprja.cpp.

References c__0, C_INT, dbnorm_(), dfnorm_(), dgbfa_, dgefa_, dls001_1, dlsa01_2, dmnorm_(), max, and min.

Referenced by CLSODA::CLSODA(), and CLSODAR::CLSODAR().

50 {
51  /* System generated locals */
52  C_INT yh_dim1, yh_offset, i__1, i__2, i__3, i__4;
53  double d__1, d__2;
54 
55  /* Local variables */
56  C_INT i__, j;
57  double r__;
58  C_INT i1, i2, j1;
59  double r0;
60  C_INT ii, jj, ml, mu;
61  double yi, yj, hl0;
62  C_INT ml3, np1;
63  double fac;
64  C_INT mba, ier;
65  double con, yjj;
66  C_INT meb1, lenp;
67  double srur;
68  C_INT mband, meband;
69 
70  /* ----------------------------------------------------------------------- */
71  /* DPRJA is called by DSTODA to compute and process the matrix */
72  /* P = I - H*EL(1)*J , where J is an approximation to the Jacobian. */
73  /* Here J is computed by the user-supplied routine JAC if */
74  /* MITER = 1 or 4 or by finite differencing if MITER = 2 or 5. */
75  /* J, scaled by -H*EL(1), is stored in WM. Then the norm of J (the */
76  /* matrix norm consistent with the weighted max-norm on vectors given */
77  /* by DMNORM) is computed, and J is overwritten by P. P is then */
78  /* subjected to LU decomposition in preparation for later solution */
79  /* of linear systems with P as coefficient matrix. This is done */
80  /* by DGEFA if MITER = 1 or 2, and by DGBFA if MITER = 4 or 5. */
81 
82  /* In addition to variables described previously, communication */
83  /* with DPRJA uses the following: */
84  /* Y = array containing predicted values on entry. */
85  /* FTEM = work array of length N (ACOR in DSTODA). */
86  /* SAVF = array containing f evaluated at predicted y. */
87  /* WM = real work space for matrices. On output it contains the */
88  /* LU decomposition of P. */
89  /* Storage of matrix elements starts at WM(3). */
90  /* WM also contains the following matrix-related data: */
91  /* WM(1) = SQRT(UROUND), used in numerical Jacobian increments. */
92  /* IWM = integer work space containing pivot information, starting at */
93  /* IWM(21). IWM also contains the band parameters */
94  /* ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. */
95  /* EL0 = EL(1) (input). */
96  /* PDNORM= norm of Jacobian matrix. (Output). */
97  /* IERPJ = output error flag, = 0 if no trouble, .gt. 0 if */
98  /* P matrix found to be singular. */
99  /* JCUR = output flag = 1 to indicate that the Jacobian matrix */
100  /* (or approximation) is now current. */
101  /* This routine also uses the Common variables EL0, H, TN, UROUND, */
102  /* MITER, N, NFE, and NJE. */
103  /* ----------------------------------------------------------------------- */
104  /* Parameter adjustments */
105  --neq;
106  --y;
107  yh_dim1 = *nyh;
108  yh_offset = 1 + yh_dim1;
109  yh -= yh_offset;
110  --ewt;
111  --ftem;
112  --savf;
113  --wm;
114  --iwm;
115 
116  /* Function Body */
117  ++dls001_1.nje;
118  dls001_1.ierpj = 0;
119  dls001_1.jcur = 1;
120  hl0 = dls001_1.h__ * dls001_1.el0;
121 
122  switch (dls001_1.miter)
123  {
124  case 1: goto L100;
125  case 2: goto L200;
126  case 3: goto L300;
127  case 4: goto L400;
128  case 5: goto L500;
129  }
130 
131  /* If MITER = 1, call JAC and multiply by scalar. ----------------------- */
132 L100:
133  lenp = dls001_1.n * dls001_1.n;
134  i__1 = lenp;
135 
136  for (i__ = 1; i__ <= i__1; ++i__)
137  {
138  /* L110: */
139  wm[i__ + 2] = 0.;
140  }
141 
142  jac(&neq[1], &dls001_1.tn, &y[1], &c__0, &c__0, &wm[3], &dls001_1.n);
143  con = -hl0;
144  i__1 = lenp;
145 
146  for (i__ = 1; i__ <= i__1; ++i__)
147  {
148  /* L120: */
149  wm[i__ + 2] *= con;
150  }
151 
152  goto L240;
153  /* If MITER = 2, make N calls to F to approximate J. -------------------- */
154 L200:
155  fac = dmnorm_(&dls001_1.n, &savf[1], &ewt[1]);
156  r0 = fabs(dls001_1.h__) * 1e3 * dls001_1.uround * dls001_1.n * fac;
157 
158  if (r0 == 0.)
159  {
160  r0 = 1.;
161  }
162 
163  srur = wm[1];
164  j1 = 2;
165  i__1 = dls001_1.n;
166 
167  for (j = 1; j <= i__1; ++j)
168  {
169  yj = y[j];
170  /* Computing MAX */
171  d__1 = srur * fabs(yj), d__2 = r0 / ewt[j];
172  r__ = std::max(d__1, d__2);
173  y[j] += r__;
174  fac = -hl0 / r__;
175  f(&neq[1], &dls001_1.tn, &y[1], &ftem[1]);
176  i__2 = dls001_1.n;
177 
178  for (i__ = 1; i__ <= i__2; ++i__)
179  {
180  /* L220: */
181  wm[i__ + j1] = (ftem[i__] - savf[i__]) * fac;
182  }
183 
184  y[j] = yj;
185  j1 += dls001_1.n;
186  /* L230: */
187  }
188 
189  dls001_1.nfe += dls001_1.n;
190 L240:
191  /* Compute norm of Jacobian. -------------------------------------------- */
192  dlsa01_2.pdnorm = dfnorm_(&dls001_1.n, &wm[3], &ewt[1]) / fabs(hl0);
193  /* Add identity matrix. ------------------------------------------------- */
194  j = 3;
195  np1 = dls001_1.n + 1;
196  i__1 = dls001_1.n;
197 
198  for (i__ = 1; i__ <= i__1; ++i__)
199  {
200  wm[j] += 1.;
201  /* L250: */
202  j += np1;
203  }
204 
205  /* Do LU decomposition on P. -------------------------------------------- */
206  dgefa_(&wm[3], &dls001_1.n, &dls001_1.n, &iwm[21], &ier);
207 
208  if (ier != 0)
209  {
210  dls001_1.ierpj = 1;
211  }
212 
213  return 0;
214  /* Dummy block only, since MITER is never 3 in this routine. ------------ */
215 L300:
216  return 0;
217  /* If MITER = 4, call JAC and multiply by scalar. ----------------------- */
218 L400:
219  ml = iwm[1];
220  mu = iwm[2];
221  ml3 = ml + 3;
222  mband = ml + mu + 1;
223  meband = mband + ml;
224  lenp = meband * dls001_1.n;
225  i__1 = lenp;
226 
227  for (i__ = 1; i__ <= i__1; ++i__)
228  {
229  /* L410: */
230  wm[i__ + 2] = 0.;
231  }
232 
233  jac(&neq[1], &dls001_1.tn, &y[1], &ml, &mu, &wm[ml3], &meband);
234  con = -hl0;
235  i__1 = lenp;
236 
237  for (i__ = 1; i__ <= i__1; ++i__)
238  {
239  /* L420: */
240  wm[i__ + 2] *= con;
241  }
242 
243  goto L570;
244  /* If MITER = 5, make MBAND calls to F to approximate J. ---------------- */
245 L500:
246  ml = iwm[1];
247  mu = iwm[2];
248  mband = ml + mu + 1;
249  mba = std::min(mband, dls001_1.n);
250  meband = mband + ml;
251  meb1 = meband - 1;
252  srur = wm[1];
253  fac = dmnorm_(&dls001_1.n, &savf[1], &ewt[1]);
254  r0 = fabs(dls001_1.h__) * 1e3 * dls001_1.uround * dls001_1.n * fac;
255 
256  if (r0 == 0.)
257  {
258  r0 = 1.;
259  }
260 
261  i__1 = mba;
262 
263  for (j = 1; j <= i__1; ++j)
264  {
265  i__2 = dls001_1.n;
266  i__3 = mband;
267 
268  for (i__ = j; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__3)
269  {
270  yi = y[i__];
271  /* Computing MAX */
272  d__1 = srur * fabs(yi), d__2 = r0 / ewt[i__];
273  r__ = std::max(d__1, d__2);
274  /* L530: */
275  y[i__] += r__;
276  }
277 
278  f(&neq[1], &dls001_1.tn, &y[1], &ftem[1]);
279  i__3 = dls001_1.n;
280  i__2 = mband;
281 
282  for (jj = j; i__2 < 0 ? jj >= i__3 : jj <= i__3; jj += i__2)
283  {
284  y[jj] = yh[jj + yh_dim1];
285  yjj = y[jj];
286  /* Computing MAX */
287  d__1 = srur * fabs(yjj), d__2 = r0 / ewt[jj];
288  r__ = std::max(d__1, d__2);
289  fac = -hl0 / r__;
290  /* Computing MAX */
291  i__4 = jj - mu;
292  i1 = std::max(i__4, (C_INT)1);
293  /* Computing MIN */
294  i__4 = jj + ml;
295  i2 = std::min(i__4, dls001_1.n);
296  ii = jj * meb1 - ml + 2;
297  i__4 = i2;
298 
299  for (i__ = i1; i__ <= i__4; ++i__)
300  {
301  /* L540: */
302  wm[ii + i__] = (ftem[i__] - savf[i__]) * fac;
303  }
304 
305  /* L550: */
306  }
307 
308  /* L560: */
309  }
310 
311  dls001_1.nfe += mba;
312 L570:
313  /* Compute norm of Jacobian. -------------------------------------------- */
314  dlsa01_2.pdnorm = dbnorm_(&dls001_1.n, &wm[ml + 3], &meband, &ml, &mu, &
315  ewt[1]) / fabs(hl0);
316  /* Add identity matrix. ------------------------------------------------- */
317  ii = mband + 2;
318  i__1 = dls001_1.n;
319 
320  for (i__ = 1; i__ <= i__1; ++i__)
321  {
322  wm[ii] += 1.;
323  /* L580: */
324  ii += meband;
325  }
326 
327  /* Do LU decomposition of P. -------------------------------------------- */
328  dgbfa_(&wm[3], &meband, &dls001_1.n, &ml, &mu, &iwm[21], &ier);
329 
330  if (ier != 0)
331  {
332  dls001_1.ierpj = 1;
333  }
334 
335  return 0;
336  /* ----------------------- End of Subroutine DPRJA ----------------------- */
337 } /* dprja_ */
#define C_INT
Definition: copasi.h:115
static C_INT c__0
Definition: dprja.cpp:37
double dbnorm_(C_INT *n, double *a, C_INT *nra, C_INT *ml, C_INT *mu, double *w)
Definition: dbnorm.cpp:31
#define dls001_1
Definition: dprja.cpp:29
#define dlsa01_2
Definition: dprja.cpp:34
#define dgefa_(__a, __lda, __n, __ipvt, __info)
Definition: dgefa.h:20
#define dgbfa_(__ab, __lda, __n, __ml, __mu, __ipvt, __info)
Definition: dgbfa.h:20
double dfnorm_(C_INT *n, double *a, double *w)
Definition: dfnorm.cpp:31
double dmnorm_(C_INT *n, double *v, double *w)
Definition: dmnorm.cpp:35
#define min(a, b)
Definition: f2c.h:175
#define max(a, b)
Definition: f2c.h:176
C_INT CInternalSolver::drchek_ ( const C_INT job,
evalG  g,
C_INT neq,
double *  y,
double *  yh,
C_INT nyh,
double *  g0,
double *  g1,
double *  gx,
C_INT jroot,
C_INT irt 
)
protected

Definition at line 46 of file drcheck.cpp.

References c__0, c__1, C_INT, d_sign(), dcopy_(), dintdy_(), dls001_1, dlsr01_1, droots_(), max, and switch().

Referenced by CLSODAR::operator()().

49 {
50  /* System generated locals */
51  C_INT yh_dim1, yh_offset, i__1;
52  double d__1;
53 
54  /* Local variables */
55  C_INT i__;
56  double x, t1, temp1, temp2;
57  C_INT iflag, jflag;
58  double hming;
59  bool zroot;
60 
61  /* ----------------------------------------------------------------------- */
62  /* This routine checks for the presence of a root in the vicinity of */
63  /* the current T, in a manner depending on the input flag JOB. It calls */
64  /* Subroutine DROOTS to locate the root as precisely as possible. */
65 
66  /* In addition to variables described previously, DRCHEK */
67  /* uses the following for communication: */
68  /* JOB = integer flag indicating type of call: */
69  /* JOB = 1 means the problem is being initialized, and DRCHEK */
70  /* is to look for a root at or very near the initial T. */
71  /* JOB = 2 means a continuation call to the solver was just */
72  /* made, and DRCHEK is to check for a root in the */
73  /* relevant part of the step last taken. */
74  /* JOB = 3 means a successful step was just taken, and DRCHEK */
75  /* is to look for a root in the interval of the step. */
76  /* G0 = array of length NG, containing the value of g at T = T0. */
77  /* G0 is input for JOB .ge. 2, and output in all cases. */
78  /* G1,GX = arrays of length NG for work space. */
79  /* IRT = completion flag: */
80  /* IRT = 0 means no root was found. */
81  /* IRT = -1 means JOB = 1 and a root was found too near to T. */
82  /* IRT = 1 means a legitimate root was found (JOB = 2 or 3). */
83  /* On return, T0 is the root location, and Y is the */
84  /* corresponding solution vector. */
85  /* T0 = value of T at one endpoint of interval of interest. Only */
86  /* roots beyond T0 in the direction of integration are sought. */
87  /* T0 is input if JOB .ge. 2, and output in all cases. */
88  /* T0 is updated by DRCHEK, whether a root is found or not. */
89  /* TLAST = last value of T returned by the solver (input only). */
90  /* TOUTC = copy of TOUT (input only). */
91  /* IRFND = input flag showing whether the last step taken had a root. */
92  /* IRFND = 1 if it did, = 0 if not. */
93  /* ITASKC = copy of ITASK (input only). */
94  /* NGC = copy of NG (input only). */
95  /* ----------------------------------------------------------------------- */
96  /* Parameter adjustments */
97  --neq;
98  --y;
99  yh_dim1 = *nyh;
100  yh_offset = 1 + yh_dim1;
101  yh -= yh_offset;
102  --g0;
103  --g1;
104  --gx;
105  --jroot;
106 
107  /* Function Body */
108  *irt = 0;
109  i__1 = dlsr01_1.ngc;
110 
111  for (i__ = 1; i__ <= i__1; ++i__)
112  {
113  /* L10: */
114  jroot[i__] = 0;
115  }
116 
117  hming = (fabs(dls001_1.tn) + fabs(dls001_1.h__)) * dls001_1.uround * 100.;
118 
119  switch (*job)
120  {
121  case 1: goto L100;
122  case 2: goto L200;
123  case 3: goto L300;
124  }
125 
126  /* Evaluate g at initial T, and check for zero values. ------------------ */
127 L100:
128  dlsr01_1.t0 = dls001_1.tn;
129  (*g)(&neq[1], &dlsr01_1.t0, &y[1], &dlsr01_1.ngc, &g0[1]);
130  dlsr01_1.nge = 1;
131  zroot = false;
132  i__1 = dlsr01_1.ngc;
133 
134  for (i__ = 1; i__ <= i__1; ++i__)
135  {
136  /* L110: */
137  if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
138  {
139  zroot = true;
140  }
141  }
142 
143  if (! zroot)
144  {
145  goto L190;
146  }
147 
148  /* g has a zero at T. Look at g at T + (small increment). -------------- */
149  /* Computing MAX */
150  d__1 = hming / fabs(dls001_1.h__);
151  temp2 = std::max(d__1, .1);
152  temp1 = temp2 * dls001_1.h__;
153  dlsr01_1.t0 += temp1;
154  i__1 = dls001_1.n;
155 
156  for (i__ = 1; i__ <= i__1; ++i__)
157  {
158  /* L120: */
159  y[i__] += temp2 * yh[i__ + (yh_dim1 << 1)];
160  }
161 
162  (*g)(&neq[1], &dlsr01_1.t0, &y[1], &dlsr01_1.ngc, &g0[1]);
163  ++dlsr01_1.nge;
164  zroot = false;
165  i__1 = dlsr01_1.ngc;
166 
167  for (i__ = 1; i__ <= i__1; ++i__)
168  {
169  /* L130: */
170  if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
171  {
172  zroot = true;
173  }
174  }
175 
176  if (! zroot)
177  {
178  goto L190;
179  }
180 
181  /* g has a zero at T and also close to T. Take error return. ----------- */
182  *irt = -1;
183  return 0;
184 
185 L190:
186  return 0;
187 
188 L200:
189 
190  if (dlsr01_1.irfnd == 0)
191  {
192  goto L260;
193  }
194 
195  /* If a root was found on the previous step, evaluate G0 = g(T0). ------- */
196  dintdy_(&dlsr01_1.t0, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
197  (*g)(&neq[1], &dlsr01_1.t0, &y[1], &dlsr01_1.ngc, &g0[1]);
198  ++dlsr01_1.nge;
199  zroot = false;
200  i__1 = dlsr01_1.ngc;
201 
202  for (i__ = 1; i__ <= i__1; ++i__)
203  {
204  /* L210: */
205  if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
206  {
207  zroot = true;
208  }
209  }
210 
211  if (! zroot)
212  {
213  goto L260;
214  }
215 
216  /* g has a zero at T0. Look at g at T + (small increment). ------------- */
217  temp1 = d_sign(hming, dls001_1.h__);
218  dlsr01_1.t0 += temp1;
219 
220  if ((dlsr01_1.t0 - dls001_1.tn) * dls001_1.h__ < 0.)
221  {
222  goto L230;
223  }
224 
225  temp2 = temp1 / dls001_1.h__;
226  i__1 = dls001_1.n;
227 
228  for (i__ = 1; i__ <= i__1; ++i__)
229  {
230  /* L220: */
231  y[i__] += temp2 * yh[i__ + (yh_dim1 << 1)];
232  }
233 
234  goto L240;
235 L230:
236  dintdy_(&dlsr01_1.t0, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
237 L240:
238  (*g)(&neq[1], &dlsr01_1.t0, &y[1], &dlsr01_1.ngc, &g0[1]);
239  ++dlsr01_1.nge;
240  zroot = false;
241  i__1 = dlsr01_1.ngc;
242 
243  for (i__ = 1; i__ <= i__1; ++i__)
244  {
245  if ((d__1 = g0[i__], fabs(d__1)) > 0.)
246  {
247  goto L250;
248  }
249 
250  jroot[i__] = 1;
251  zroot = true;
252 L250:
253  ;
254  }
255 
256  if (! zroot)
257  {
258  goto L260;
259  }
260 
261  /* g has a zero at T0 and also close to T0. Return root. --------------- */
262  *irt = 1;
263  return 0;
264  /* G0 has no zero components. Proceed to check relevant interval. ------ */
265 L260:
266 
267  if (dls001_1.tn == dlsr01_1.tlast)
268  {
269  goto L390;
270  }
271 
272 L300:
273 
274  /* Set T1 to TN or TOUTC, whichever comes first, and get g at T1. ------- */
275  if (dlsr01_1.itaskc == 2 || dlsr01_1.itaskc == 3 || dlsr01_1.itaskc == 5)
276  {
277  goto L310;
278  }
279 
280  if ((dlsr01_1.toutc - dls001_1.tn) * dls001_1.h__ >= 0.)
281  {
282  goto L310;
283  }
284 
285  t1 = dlsr01_1.toutc;
286 
287  if ((t1 - dlsr01_1.t0) * dls001_1.h__ <= 0.)
288  {
289  goto L390;
290  }
291 
292  dintdy_(&t1, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
293  goto L330;
294 L310:
295  t1 = dls001_1.tn;
296  i__1 = dls001_1.n;
297 
298  for (i__ = 1; i__ <= i__1; ++i__)
299  {
300  /* L320: */
301  y[i__] = yh[i__ + yh_dim1];
302  }
303 
304 L330:
305  (*g)(&neq[1], &t1, &y[1], &dlsr01_1.ngc, &g1[1]);
306  ++dlsr01_1.nge;
307  /* Call DROOTS to search for root in interval from T0 to T1. ------------ */
308  jflag = 0;
309 L350:
310  droots_(&dlsr01_1.ngc, &hming, &jflag, &dlsr01_1.t0, &t1, &g0[1], &g1[1],
311  &gx[1], &x, &jroot[1]);
312 
313  if (jflag > 1)
314  {
315  goto L360;
316  }
317 
318  dintdy_(&x, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
319  (*g)(&neq[1], &x, &y[1], &dlsr01_1.ngc, &gx[1]);
320  ++dlsr01_1.nge;
321  goto L350;
322 L360:
323  dlsr01_1.t0 = x;
324  dcopy_(&dlsr01_1.ngc, &gx[1], &c__1, &g0[1], &c__1);
325 
326  if (jflag == 4)
327  {
328  goto L390;
329  }
330 
331  /* Found a root. Interpolate to X and return. -------------------------- */
332  dintdy_(&x, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
333  *irt = 1;
334  return 0;
335 
336 L390:
337  return 0;
338  /* ----------------------- End of Subroutine DRCHEK ---------------------- */
339 } /* drchek_ */
#define C_INT
Definition: copasi.h:115
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
C_INT droots_(C_INT *ng, double *hmin, C_INT *jflag, double *x0, double *x1, double *g0, double *g1, double *gx, double *x, C_INT *jroot)
Definition: drcheck.cpp:343
static C_INT c__1
Definition: drcheck.cpp:42
C_INT dintdy_(double *t, const C_INT *k, double *yh, C_INT *nyh, double *dky, C_INT *iflag)
Definition: dintdy.cpp:53
#define dlsr01_1
Definition: drcheck.cpp:33
switch(yytype)
#define dls001_1
Definition: drcheck.cpp:29
double d_sign(const double &a, const double &b)
Definition: CLSODA.cpp:2471
static const C_INT c__0
Definition: drcheck.cpp:41
#define max(a, b)
Definition: f2c.h:176
C_INT CInternalSolver::droots_ ( C_INT ng,
double *  hmin,
C_INT jflag,
double *  x0,
double *  x1,
double *  g0,
double *  g1,
double *  gx,
double *  x,
C_INT jroot 
)
protected

Definition at line 343 of file drcheck.cpp.

References c__1, c_b3, C_INT, d_sign(), dcopy_(), and dlsr01_2.

Referenced by drchek_().

346 {
347  /* Initialized data */
348 
349  static double zero = 0.;
350  static double half = .5;
351  static double tenth = .1;
352  static double five = 5.;
353 
354  /* System generated locals */
355  C_INT i__1;
356  double d__1;
357 
358  /* Local variables */
359  C_INT i__;
360  double t2, tmax;
361  bool xroot, zroot, sgnchg;
362  C_INT imxold, nxlast;
363  double fracsub, fracint;
364 
365  /* ----------------------------------------------------------------------- */
366  /* This subroutine finds the leftmost root of a set of arbitrary */
367  /* functions gi(x) (i = 1,...,NG) in an interval (X0,X1). Only roots */
368  /* of odd multiplicity (i.e. changes of sign of the gi) are found. */
369  /* Here the sign of X1 - X0 is arbitrary, but is constant for a given */
370  /* problem, and -leftmost- means nearest to X0. */
371  /* The values of the vector-valued function g(x) = (gi, i=1...NG) */
372  /* are communicated through the call sequence of DROOTS. */
373  /* The method used is the Illinois algorithm. */
374 
375  /* Reference: */
376  /* Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined */
377  /* Output Points for Solutions of ODEs, Sandia Report SAND80-0180, */
378  /* February 1980. */
379 
380  /* Description of parameters. */
381 
382  /* NG = number of functions gi, or the number of components of */
383  /* the vector valued function g(x). Input only. */
384 
385  /* HMIN = resolution parameter in X. Input only. When a root is */
386  /* found, it is located only to within an error of HMIN in X. */
387  /* Typically, HMIN should be set to something on the order of */
388  /* 100 * UROUND * MAX(ABS(X0),ABS(X1)), */
389  /* where UROUND is the unit roundoff of the machine. */
390 
391  /* JFLAG = integer flag for input and output communication. */
392 
393  /* On input, set JFLAG = 0 on the first call for the problem, */
394  /* and leave it unchanged until the problem is completed. */
395  /* (The problem is completed when JFLAG .ge. 2 on return.) */
396 
397  /* On output, JFLAG has the following values and meanings: */
398  /* JFLAG = 1 means DROOTS needs a value of g(x). Set GX = g(X) */
399  /* and call DROOTS again. */
400  /* JFLAG = 2 means a root has been found. The root is */
401  /* at X, and GX contains g(X). (Actually, X is the */
402  /* rightmost approximation to the root on an interval */
403  /* (X0,X1) of size HMIN or less.) */
404  /* JFLAG = 3 means X = X1 is a root, with one or more of the gi */
405  /* being zero at X1 and no sign changes in (X0,X1). */
406  /* GX contains g(X) on output. */
407  /* JFLAG = 4 means no roots (of odd multiplicity) were */
408  /* found in (X0,X1) (no sign changes). */
409 
410  /* X0,X1 = endpoints of the interval where roots are sought. */
411  /* X1 and X0 are input when JFLAG = 0 (first call), and */
412  /* must be left unchanged between calls until the problem is */
413  /* completed. X0 and X1 must be distinct, but X1 - X0 may be */
414  /* of either sign. However, the notion of -left- and -right- */
415  /* will be used to mean nearer to X0 or X1, respectively. */
416  /* When JFLAG .ge. 2 on return, X0 and X1 are output, and */
417  /* are the endpoints of the relevant interval. */
418 
419  /* G0,G1 = arrays of length NG containing the vectors g(X0) and g(X1), */
420  /* respectively. When JFLAG = 0, G0 and G1 are input and */
421  /* none of the G0(i) should be zero. */
422  /* When JFLAG .ge. 2 on return, G0 and G1 are output. */
423 
424  /* GX = array of length NG containing g(X). GX is input */
425  /* when JFLAG = 1, and output when JFLAG .ge. 2. */
426 
427  /* X = independent variable value. Output only. */
428  /* When JFLAG = 1 on output, X is the point at which g(x) */
429  /* is to be evaluated and loaded into GX. */
430  /* When JFLAG = 2 or 3, X is the root. */
431  /* When JFLAG = 4, X is the right endpoint of the interval, X1. */
432 
433  /* JROOT = integer array of length NG. Output only. */
434  /* When JFLAG = 2 or 3, JROOT indicates which components */
435  /* of g(x) have a root at X. JROOT(i) is 1 if the i-th */
436  /* component has a root, and JROOT(i) = 0 otherwise. */
437  /* ----------------------------------------------------------------------- */
438  /* Parameter adjustments */
439  --jroot;
440  --gx;
441  --g1;
442  --g0;
443 
444  /* Function Body */
445 
446  if (*jflag == 1)
447  {
448  goto L200;
449  }
450 
451  /* JFLAG .ne. 1. Check for change in sign of g or zero at X1. ---------- */
452  dlsr01_2.imax = 0;
453  tmax = zero;
454  zroot = false;
455  i__1 = *ng;
456 
457  for (i__ = 1; i__ <= i__1; ++i__)
458  {
459  if ((d__1 = g1[i__], fabs(d__1)) > zero)
460  {
461  goto L110;
462  }
463 
464  zroot = true;
465  goto L120;
466  /* At this point, G0(i) has been checked and cannot be zero. ------------ */
467 L110:
468 
469  if (d_sign(c_b3, g0[i__]) == d_sign(c_b3, g1[i__]))
470  {
471  goto L120;
472  }
473 
474  t2 = (d__1 = g1[i__] / (g1[i__] - g0[i__]), fabs(d__1));
475 
476  if (t2 <= tmax)
477  {
478  goto L120;
479  }
480 
481  tmax = t2;
482  dlsr01_2.imax = i__;
483 L120:
484  ;
485  }
486 
487  if (dlsr01_2.imax > 0)
488  {
489  goto L130;
490  }
491 
492  sgnchg = false;
493  goto L140;
494 L130:
495  sgnchg = true;
496 L140:
497 
498  if (! sgnchg)
499  {
500  goto L400;
501  }
502 
503  /* There is a sign change. Find the first root in the interval. -------- */
504  xroot = false;
505  nxlast = 0;
506  dlsr01_2.last = 1;
507 
508  /* Repeat until the first root in the interval is found. Loop point. --- */
509 L150:
510 
511  if (xroot)
512  {
513  goto L300;
514  }
515 
516  if (nxlast == dlsr01_2.last)
517  {
518  goto L160;
519  }
520 
521  dlsr01_2.alpha = 1.;
522  goto L180;
523 L160:
524 
525  if (dlsr01_2.last == 0)
526  {
527  goto L170;
528  }
529 
530  dlsr01_2.alpha *= .5;
531  goto L180;
532 L170:
533  dlsr01_2.alpha *= 2.;
534 L180:
535  dlsr01_2.x2 = *x1 - (*x1 - *x0) * g1[dlsr01_2.imax] / (g1[dlsr01_2.imax]
536  - dlsr01_2.alpha * g0[dlsr01_2.imax]);
537 
538  /* If X2 is too close to X0 or X1, adjust it inward, by a fractional ---- */
539 
540  /* distance that is between 0.1 and 0.5. -------------------------------- */
541  if ((d__1 = dlsr01_2.x2 - *x0, fabs(d__1)) < half * *hmin)
542  {
543  fracint = (d__1 = *x1 - *x0, fabs(d__1)) / *hmin;
544  fracsub = tenth;
545 
546  if (fracint <= five)
547  {
548  fracsub = half / fracint;
549  }
550 
551  dlsr01_2.x2 = *x0 + fracsub * (*x1 - *x0);
552  }
553 
554  if ((d__1 = *x1 - dlsr01_2.x2, fabs(d__1)) < half * *hmin)
555  {
556  fracint = (d__1 = *x1 - *x0, fabs(d__1)) / *hmin;
557  fracsub = tenth;
558 
559  if (fracint <= five)
560  {
561  fracsub = half / fracint;
562  }
563 
564  dlsr01_2.x2 = *x1 - fracsub * (*x1 - *x0);
565  }
566 
567  *jflag = 1;
568  *x = dlsr01_2.x2;
569  /* Return to the calling routine to get a value of GX = g(X). ----------- */
570  return 0;
571  /* Check to see in which interval g changes sign. ----------------------- */
572 L200:
573  imxold = dlsr01_2.imax;
574  dlsr01_2.imax = 0;
575  tmax = zero;
576  zroot = false;
577  i__1 = *ng;
578 
579  for (i__ = 1; i__ <= i__1; ++i__)
580  {
581  if ((d__1 = gx[i__], fabs(d__1)) > zero)
582  {
583  goto L210;
584  }
585 
586  zroot = true;
587  goto L220;
588  /* Neither G0(i) nor GX(i) can be zero at this point. ------------------- */
589 L210:
590 
591  if (d_sign(c_b3, g0[i__]) == d_sign(c_b3, gx[i__]))
592  {
593  goto L220;
594  }
595 
596  t2 = (d__1 = gx[i__] / (gx[i__] - g0[i__]), fabs(d__1));
597 
598  if (t2 <= tmax)
599  {
600  goto L220;
601  }
602 
603  tmax = t2;
604  dlsr01_2.imax = i__;
605 L220:
606  ;
607  }
608 
609  if (dlsr01_2.imax > 0)
610  {
611  goto L230;
612  }
613 
614  sgnchg = false;
615  dlsr01_2.imax = imxold;
616  goto L240;
617 L230:
618  sgnchg = true;
619 L240:
620  nxlast = dlsr01_2.last;
621 
622  if (! sgnchg)
623  {
624  goto L250;
625  }
626 
627  /* Sign change between X0 and X2, so replace X1 with X2. ---------------- */
628  *x1 = dlsr01_2.x2;
629  dcopy_(ng, &gx[1], &c__1, &g1[1], &c__1);
630  dlsr01_2.last = 1;
631  xroot = false;
632  goto L270;
633 L250:
634 
635  if (! zroot)
636  {
637  goto L260;
638  }
639 
640  /* Zero value at X2 and no sign change in (X0,X2), so X2 is a root. ----- */
641  *x1 = dlsr01_2.x2;
642  dcopy_(ng, &gx[1], &c__1, &g1[1], &c__1);
643  xroot = true;
644  goto L270;
645  /* No sign change between X0 and X2. Replace X0 with X2. --------------- */
646 L260:
647  dcopy_(ng, &gx[1], &c__1, &g0[1], &c__1);
648  *x0 = dlsr01_2.x2;
649  dlsr01_2.last = 0;
650  xroot = false;
651 L270:
652 
653  if ((d__1 = *x1 - *x0, fabs(d__1)) <= *hmin)
654  {
655  xroot = true;
656  }
657 
658  goto L150;
659 
660  /* Return with X1 as the root. Set JROOT. Set X = X1 and GX = G1. ----- */
661 L300:
662  *jflag = 2;
663  *x = *x1;
664  dcopy_(ng, &g1[1], &c__1, &gx[1], &c__1);
665  i__1 = *ng;
666 
667  for (i__ = 1; i__ <= i__1; ++i__)
668  {
669  jroot[i__] = 0;
670 
671  if ((d__1 = g1[i__], fabs(d__1)) > zero)
672  {
673  goto L310;
674  }
675 
676  jroot[i__] = 1;
677  goto L320;
678 L310:
679 
680  if (d_sign(c_b3, g0[i__]) != d_sign(c_b3, g1[i__]))
681  {
682  jroot[i__] = 1;
683  }
684 
685 L320:
686  ;
687  }
688 
689  return 0;
690 
691  /* No sign change in the interval. Check for zero at right endpoint. --- */
692 L400:
693 
694  if (! zroot)
695  {
696  goto L420;
697  }
698 
699  /* Zero value at X1 and no sign change in (X0,X1). Return JFLAG = 3. --- */
700  *x = *x1;
701  dcopy_(ng, &g1[1], &c__1, &gx[1], &c__1);
702  i__1 = *ng;
703 
704  for (i__ = 1; i__ <= i__1; ++i__)
705  {
706  jroot[i__] = 0;
707 
708  if ((d__1 = g1[i__], fabs(d__1)) <= zero)
709  {
710  jroot[i__] = 1;
711  }
712 
713  /* L410: */
714  }
715 
716  *jflag = 3;
717  return 0;
718 
719  /* No sign changes in this interval. Set X = X1, return JFLAG = 4. ----- */
720 L420:
721  dcopy_(ng, &g1[1], &c__1, &gx[1], &c__1);
722  *x = *x1;
723  *jflag = 4;
724  return 0;
725  /* ----------------------- End of Subroutine DROOTS ---------------------- */
726 } /* droots_ */
#define C_INT
Definition: copasi.h:115
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
static C_INT c__1
Definition: drcheck.cpp:42
#define dlsr01_2
Definition: drcheck.cpp:34
static double c_b3
Definition: drcheck.cpp:39
double d_sign(const double &a, const double &b)
Definition: CLSODA.cpp:2471
C_INT CInternalSolver::dsolsy_ ( double *  wm,
C_INT iwm,
double *  x,
double *  tem 
)
protected

Definition at line 43 of file dsolsy.cpp.

References c__0, C_INT, dgbsl_, dgesl_, and dls001_1.

Referenced by CLSODA::CLSODA(), and CLSODAR::CLSODAR().

45 {
46  /* System generated locals */
47  C_INT i__1;
48 
49  /* Local variables */
50  C_INT i__;
51  double r__, di;
52  C_INT ml, mu;
53  double hl0, phl0;
54  C_INT meband;
55 
56  /* ***BEGIN PROLOGUE DSOLSY */
57  /* ***SUBSIDIARY */
58  /* ***PURPOSE ODEPACK linear system solver. */
59  /* ***TYPE DOUBLE PRECISION (SSOLSY-S, DSOLSY-D) */
60  /* ***AUTHOR Hindmarsh, Alan C., (LLNL) */
61  /* ***DESCRIPTION */
62 
63  /* This routine manages the solution of the linear system arising from */
64  /* a chord iteration. It is called if MITER .ne. 0. */
65  /* If MITER is 1 or 2, it calls DGESL to accomplish this. */
66  /* If MITER = 3 it updates the coefficient h*EL0 in the diagonal */
67  /* matrix, and then computes the solution. */
68  /* If MITER is 4 or 5, it calls DGBSL. */
69  /* Communication with DSOLSY uses the following variables: */
70  /* WM = real work space containing the inverse diagonal matrix if */
71  /* MITER = 3 and the LU decomposition of the matrix otherwise. */
72  /* Storage of matrix elements starts at WM(3). */
73  /* WM also contains the following matrix-related data: */
74  /* WM(1) = SQRT(UROUND) (not used here), */
75  /* WM(2) = HL0, the previous value of h*EL0, used if MITER = 3. */
76  /* IWM = integer work space containing pivot information, starting at */
77  /* IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band */
78  /* parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. */
79  /* X = the right-hand side vector on input, and the solution vector */
80  /* on output, of length N. */
81  /* TEM = vector of work space of length N, not used in this version. */
82  /* IERSL = output flag (in COMMON). IERSL = 0 if no trouble occurred. */
83  /* IERSL = 1 if a singular matrix arose with MITER = 3. */
84  /* This routine also uses the COMMON variables EL0, H, MITER, and N. */
85 
86  /* ***SEE ALSO DLSODE */
87  /* ***ROUTINES CALLED DGBSL, DGESL */
88  /* ***COMMON BLOCKS DLS001 */
89  /* ***REVISION HISTORY (YYMMDD) */
90  /* 791129 DATE WRITTEN */
91  /* 890501 Modified prologue to SLATEC/LDOC format. (FNF) */
92  /* 890503 Minor cosmetic changes. (FNF) */
93  /* 930809 Renamed to allow single/double precision versions. (ACH) */
94  /* 010418 Reduced size of Common block /DLS001/. (ACH) */
95  /* 031105 Restored 'own' variables to Common block /DLS001/, to */
96  /* enable interrupt/restart feature. (ACH) */
97  /* ***END PROLOGUE DSOLSY */
98  /* **End */
99 
100  /* ***FIRST EXECUTABLE STATEMENT DSOLSY */
101  /* Parameter adjustments */
102  --tem;
103  --x;
104  --iwm;
105  --wm;
106 
107  /* Function Body */
108  dls001_1.iersl = 0;
109 
110  switch (dls001_1.miter)
111  {
112  case 1: goto L100;
113  case 2: goto L100;
114  case 3: goto L300;
115  case 4: goto L400;
116  case 5: goto L400;
117  }
118 
119 L100:
120  dgesl_(&wm[3], &dls001_1.n, &dls001_1.n, &iwm[21], &x[1], &c__0);
121  return 0;
122 
123 L300:
124  phl0 = wm[2];
125  hl0 = dls001_1.h__ * dls001_1.el0;
126  wm[2] = hl0;
127 
128  if (hl0 == phl0)
129  {
130  goto L330;
131  }
132 
133  r__ = hl0 / phl0;
134  i__1 = dls001_1.n;
135 
136  for (i__ = 1; i__ <= i__1; ++i__)
137  {
138  di = 1. - r__ * (1. - 1. / wm[i__ + 2]);
139 
140  if (fabs(di) == 0.)
141  {
142  goto L390;
143  }
144 
145  /* L320: */
146  wm[i__ + 2] = 1. / di;
147  }
148 
149 L330:
150  i__1 = dls001_1.n;
151 
152  for (i__ = 1; i__ <= i__1; ++i__)
153  {
154  /* L340: */
155  x[i__] = wm[i__ + 2] * x[i__];
156  }
157 
158  return 0;
159 L390:
160  dls001_1.iersl = 1;
161  return 0;
162 
163 L400:
164  ml = iwm[1];
165  mu = iwm[2];
166  meband = (ml << 1) + mu + 1;
167  dgbsl_(&wm[3], &meband, &dls001_1.n, &ml, &mu, &iwm[21], &x[1], &c__0);
168  return 0;
169  /* ----------------------- END OF SUBROUTINE DSOLSY ---------------------- */
170 } /* dsolsy_ */
#define C_INT
Definition: copasi.h:115
#define dgesl_(__a, __lda, __n, __ipvt, __b, __job)
Definition: dgesl.h:20
static C_INT c__0
Definition: dsolsy.cpp:34
#define dls001_1
Definition: dsolsy.cpp:26
#define dgbsl_(__ab, __lda, __n, __ml, __mu, __ipvt, __b, __job)
Definition: dgbsl.h:20
C_INT CInternalSolver::dstoda_ ( C_INT neq,
double *  y,
double *  yh,
C_INT nyh,
double *  yh1,
double *  ewt,
double *  savf,
double *  acor,
double *  wm,
C_INT iwm,
evalF  f,
evalJ  jac,
PJAC pjac,
SLVS slvs 
)
protected

Definition at line 45 of file dstoda.cpp.

References c__1, c__2, C_INT, dcfode_(), dls001_3, dlsa01_1, dmnorm_(), max, min, and pow_dd.

Referenced by CLSODA::operator()(), and CLSODAR::operator()().

49 {
50  /* Initialized data */
51 
52  static const double sm1[12] =
53  {
54  .5, .575, .55, .45, .35, .25, .2, .15, .1, .075, .05,
55  .025
56  };
57 
58  /* System generated locals */
59  C_INT yh_dim1, yh_offset, i__1, i__2;
60  double d__1, d__2, d__3;
61 
62  /* Local variables */
63  C_INT i__, j, m;
64  double r__;
65  C_INT i1, jb;
66  double rh = 0.0, rm, dm1, dm2;
67  C_INT lm1, lm2;
68  double rh1, rh2, del, ddn;
69  C_INT ncf;
70  double pdh = 0.0, dsm, dup, exm1, exm2;
71  C_INT nqm1, nqm2;
72  double dcon, delp;
73  C_INT lm1p1, lm2p1;
74  double exdn, rhdn;
75  C_INT iret;
76  double told, rhsm;
77  C_INT newq;
78  double exsm, rhup, rate, exup, rh1it, alpha;
79  C_INT iredo = 0;
80  double pnorm;
81 
82  /* Parameter adjustments */
83  --neq;
84  --y;
85  yh_dim1 = *nyh;
86  yh_offset = 1 + yh_dim1;
87  yh -= yh_offset;
88  --yh1;
89  --ewt;
90  --savf;
91  --acor;
92  --wm;
93  --iwm;
94 
95  /* Function Body */
96  /* ----------------------------------------------------------------------- */
97  /* DSTODA performs one step of the integration of an initial value */
98  /* problem for a system of ordinary differential equations. */
99  /* Note: DSTODA is independent of the value of the iteration method */
100  /* indicator MITER, when this is .ne. 0, and hence is independent */
101  /* of the type of chord method used, or the Jacobian structure. */
102  /* Communication with DSTODA is done with the following variables: */
103 
104  /* Y = an array of length .ge. N used as the Y argument in */
105  /* all calls to F and JAC. */
106  /* NEQ = integer array containing problem size in NEQ(1), and */
107  /* passed as the NEQ argument in all calls to F and JAC. */
108  /* YH = an NYH by LMAX array containing the dependent variables */
109  /* and their approximate scaled derivatives, where */
110  /* LMAX = MAXORD + 1. YH(i,j+1) contains the approximate */
111  /* j-th derivative of y(i), scaled by H**j/factorial(j) */
112  /* (j = 0,1,...,NQ). On entry for the first step, the first */
113  /* two columns of YH must be set from the initial values. */
114  /* NYH = a constant integer .ge. N, the first dimension of YH. */
115  /* YH1 = a one-dimensional array occupying the same space as YH. */
116  /* EWT = an array of length N containing multiplicative weights */
117  /* for local error measurements. Local errors in y(i) are */
118  /* compared to 1.0/EWT(i) in various error tests. */
119  /* SAVF = an array of working storage, of length N. */
120  /* ACOR = a work array of length N, used for the accumulated */
121  /* corrections. On a successful return, ACOR(i) contains */
122  /* the estimated one-step local error in y(i). */
123  /* WM,IWM = real and integer work arrays associated with matrix */
124  /* operations in chord iteration (MITER .ne. 0). */
125  /* PJAC = name of routine to evaluate and preprocess Jacobian matrix */
126  /* and P = I - H*EL0*Jac, if a chord method is being used. */
127  /* It also returns an estimate of norm(Jac) in PDNORM. */
128  /* SLVS = name of routine to solve linear system in chord iteration. */
129  /* CCMAX = maximum relative change in H*EL0 before PJAC is called. */
130  /* H = the step size to be attempted on the next step. */
131  /* H is altered by the error control algorithm during the */
132  /* problem. H can be either positive or negative, but its */
133  /* sign must remain constant throughout the problem. */
134  /* HMIN = the minimum absolute value of the step size H to be used. */
135  /* HMXI = inverse of the maximum absolute value of H to be used. */
136  /* HMXI = 0.0 is allowed and corresponds to an infinite HMAX. */
137  /* HMIN and HMXI may be changed at any time, but will not */
138  /* take effect until the next change of H is considered. */
139  /* TN = the independent variable. TN is updated on each step taken. */
140  /* JSTART = an integer used for input only, with the following */
141  /* values and meanings: */
142  /* 0 perform the first step. */
143  /* .gt.0 take a new step continuing from the last. */
144  /* -1 take the next step with a new value of H, */
145  /* N, METH, MITER, and/or matrix parameters. */
146  /* -2 take the next step with a new value of H, */
147  /* but with other inputs unchanged. */
148  /* On return, JSTART is set to 1 to facilitate continuation. */
149  /* KFLAG = a completion code with the following meanings: */
150  /* 0 the step was succesful. */
151  /* -1 the requested error could not be achieved. */
152  /* -2 corrector convergence could not be achieved. */
153  /* -3 fatal error in PJAC or SLVS. */
154  /* A return with KFLAG = -1 or -2 means either */
155  /* ABS(H) = HMIN or 10 consecutive failures occurred. */
156  /* On a return with KFLAG negative, the values of TN and */
157  /* the YH array are as of the beginning of the last */
158  /* step, and H is the last step size attempted. */
159  /* MAXORD = the maximum order of integration method to be allowed. */
160  /* MAXCOR = the maximum number of corrector iterations allowed. */
161  /* MSBP = maximum number of steps between PJAC calls (MITER .gt. 0). */
162  /* MXNCF = maximum number of convergence failures allowed. */
163  /* METH = current method. */
164  /* METH = 1 means Adams method (nonstiff) */
165  /* METH = 2 means BDF method (stiff) */
166  /* METH may be reset by DSTODA. */
167  /* MITER = corrector iteration method. */
168  /* MITER = 0 means functional iteration. */
169  /* MITER = JT .gt. 0 means a chord iteration corresponding */
170  /* to Jacobian type JT. (The DLSODA/DLSODAR argument JT is */
171  /* communicated here as JTYP, but is not used in DSTODA */
172  /* except to load MITER following a method switch.) */
173  /* MITER may be reset by DSTODA. */
174  /* N = the number of first-order differential equations. */
175  /* ----------------------------------------------------------------------- */
176  dls001_3.kflag = 0;
177  told = dls001_3.tn;
178  ncf = 0;
179  dls001_3.ierpj = 0;
180  dls001_3.iersl = 0;
181  dls001_3.jcur = 0;
182  dls001_3.icf = 0;
183  delp = 0.;
184 
185  if (dls001_3.jstart > 0)
186  {
187  goto L200;
188  }
189 
190  if (dls001_3.jstart == -1)
191  {
192  goto L100;
193  }
194 
195  if (dls001_3.jstart == -2)
196  {
197  goto L160;
198  }
199 
200  /* ----------------------------------------------------------------------- */
201  /* On the first call, the order is set to 1, and other variables are */
202  /* initialized. RMAX is the maximum ratio by which H can be increased */
203  /* in a single step. It is initially 1.E4 to compensate for the small */
204  /* initial H, but then is normally equal to 10. If a failure */
205  /* occurs (in corrector convergence or error test), RMAX is set at 2 */
206  /* for the next increase. */
207  /* DCFODE is called to get the needed coefficients for both methods. */
208  /* ----------------------------------------------------------------------- */
209  dls001_3.lmax = dls001_3.maxord + 1;
210  dls001_3.nq = 1;
211  dls001_3.l = 2;
212  dls001_3.ialth = 2;
213  dls001_3.rmax = 1e4;
214  dls001_3.rc = 0.;
215  dls001_3.el0 = 1.;
216  dls001_3.crate = .7;
217  dls001_3.hold = dls001_3.h__;
218  dls001_3.nslp = 0;
219  dls001_3.ipup = dls001_3.miter;
220  iret = 3;
221  /* Initialize switching parameters. METH = 1 is assumed initially. ----- */
222  dlsa01_1.icount = 20;
223  dlsa01_1.irflag = 0;
224  dlsa01_1.pdest = 0.;
225  dlsa01_1.pdlast = 0.;
226  dlsa01_1.ratio = 5.;
227  dcfode_(&c__2, dls001_3.elco, dls001_3.tesco);
228 
229  for (i__ = 1; i__ <= 5; ++i__)
230  {
231  /* L10: */
232  dlsa01_1.cm2[i__ - 1] = dls001_3.tesco[i__ * 3 - 2] * dls001_3.elco[
233  i__ + 1 + i__ * 13 - 14];
234  }
235 
236  dcfode_(&c__1, dls001_3.elco, dls001_3.tesco);
237 
238  for (i__ = 1; i__ <= 12; ++i__)
239  {
240  /* L20: */
241  dlsa01_1.cm1[i__ - 1] = dls001_3.tesco[i__ * 3 - 2] * dls001_3.elco[
242  i__ + 1 + i__ * 13 - 14];
243  }
244 
245  goto L150;
246  /* ----------------------------------------------------------------------- */
247  /* The following block handles preliminaries needed when JSTART = -1. */
248  /* IPUP is set to MITER to force a matrix update. */
249  /* If an order increase is about to be considered (IALTH = 1), */
250  /* IALTH is reset to 2 to postpone consideration one more step. */
251  /* If the caller has changed METH, DCFODE is called to reset */
252  /* the coefficients of the method. */
253  /* If H is to be changed, YH must be rescaled. */
254  /* If H or METH is being changed, IALTH is reset to L = NQ + 1 */
255  /* to prevent further changes in H for that many steps. */
256  /* ----------------------------------------------------------------------- */
257 L100:
258  dls001_3.ipup = dls001_3.miter;
259  dls001_3.lmax = dls001_3.maxord + 1;
260 
261  if (dls001_3.ialth == 1)
262  {
263  dls001_3.ialth = 2;
264  }
265 
266  if (dls001_3.meth == dlsa01_1.mused)
267  {
268  goto L160;
269  }
270 
271  dcfode_(&dls001_3.meth, dls001_3.elco, dls001_3.tesco);
272  dls001_3.ialth = dls001_3.l;
273  iret = 1;
274  /* ----------------------------------------------------------------------- */
275  /* The el vector and related constants are reset */
276  /* whenever the order NQ is changed, or at the start of the problem. */
277  /* ----------------------------------------------------------------------- */
278 L150:
279  i__1 = dls001_3.l;
280 
281  for (i__ = 1; i__ <= i__1; ++i__)
282  {
283  /* L155: */
284  dls001_3.el[i__ - 1] = dls001_3.elco[i__ + dls001_3.nq * 13 - 14];
285  }
286 
287  dls001_3.nqnyh = dls001_3.nq * *nyh;
288  dls001_3.rc = dls001_3.rc * dls001_3.el[0] / dls001_3.el0;
289  dls001_3.el0 = dls001_3.el[0];
290  dls001_3.conit = .5 / (dls001_3.nq + 2);
291 
292  switch (iret)
293  {
294  case 1: goto L160;
295  case 2: goto L170;
296  case 3: goto L200;
297  }
298 
299  /* ----------------------------------------------------------------------- */
300  /* If H is being changed, the H ratio RH is checked against */
301  /* RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to */
302  /* L = NQ + 1 to prevent a change of H for that many steps, unless */
303  /* forced by a convergence or error test failure. */
304  /* ----------------------------------------------------------------------- */
305 L160:
306 
307  if (dls001_3.h__ == dls001_3.hold)
308  {
309  goto L200;
310  }
311 
312  rh = dls001_3.h__ / dls001_3.hold;
313  dls001_3.h__ = dls001_3.hold;
314  iredo = 3;
315  goto L175;
316 L170:
317  /* Computing MAX */
318  d__1 = rh, d__2 = dls001_3.hmin / fabs(dls001_3.h__);
319  rh = std::max(d__1, d__2);
320 L175:
321  rh = std::min(rh, dls001_3.rmax);
322  /* Computing MAX */
323  d__1 = 1., d__2 = fabs(dls001_3.h__) * dls001_3.hmxi * rh;
324  rh /= std::max(d__1, d__2);
325 
326  /* ----------------------------------------------------------------------- */
327 
328  /* If METH = 1, also restrict the new step size by the stability region. */
329 
330  /* If this reduces H, set IRFLAG to 1 so that if there are roundoff */
331 
332  /* problems later, we can assume that is the cause of the trouble. */
333 
334  /* ----------------------------------------------------------------------- */
335  if (dls001_3.meth == 2)
336  {
337  goto L178;
338  }
339 
340  dlsa01_1.irflag = 0;
341  /* Computing MAX */
342  d__1 = fabs(dls001_3.h__) * dlsa01_1.pdlast;
343  pdh = std::max(d__1, 1e-6);
344 
345  if (rh * pdh * 1.00001 < sm1[dls001_3.nq - 1])
346  {
347  goto L178;
348  }
349 
350  rh = sm1[dls001_3.nq - 1] / pdh;
351  dlsa01_1.irflag = 1;
352 L178:
353  r__ = 1.;
354  i__1 = dls001_3.l;
355 
356  for (j = 2; j <= i__1; ++j)
357  {
358  r__ *= rh;
359  i__2 = dls001_3.n;
360 
361  for (i__ = 1; i__ <= i__2; ++i__)
362  {
363  /* L180: */
364  yh[i__ + j * yh_dim1] *= r__;
365  }
366  }
367 
368  dls001_3.h__ *= rh;
369  dls001_3.rc *= rh;
370  dls001_3.ialth = dls001_3.l;
371 
372  if (iredo == 0)
373  {
374  goto L690;
375  }
376 
377  /* ----------------------------------------------------------------------- */
378  /* This section computes the predicted values by effectively */
379  /* multiplying the YH array by the Pascal triangle matrix. */
380  /* RC is the ratio of new to old values of the coefficient H*EL(1). */
381  /* When RC differs from 1 by more than CCMAX, IPUP is set to MITER */
382  /* to force PJAC to be called, if a Jacobian is involved. */
383  /* In any case, PJAC is called at least every MSBP steps. */
384  /* ----------------------------------------------------------------------- */
385 L200:
386 
387  if ((d__1 = dls001_3.rc - 1., fabs(d__1)) > dls001_3.ccmax)
388  {
389  dls001_3.ipup = dls001_3.miter;
390  }
391 
392  if (dls001_3.nst >= dls001_3.nslp + dls001_3.msbp)
393  {
394  dls001_3.ipup = dls001_3.miter;
395  }
396 
397  dls001_3.tn += dls001_3.h__;
398  i1 = dls001_3.nqnyh + 1;
399  i__2 = dls001_3.nq;
400 
401  for (jb = 1; jb <= i__2; ++jb)
402  {
403  i1 -= *nyh;
404  /* DIR$ IVDEP */
405  i__1 = dls001_3.nqnyh;
406 
407  for (i__ = i1; i__ <= i__1; ++i__)
408  {
409  /* L210: */
410  yh1[i__] += yh1[i__ + *nyh];
411  }
412 
413  /* L215: */
414  }
415 
416  pnorm = dmnorm_(&dls001_3.n, &yh1[1], &ewt[1]);
417  /* ----------------------------------------------------------------------- */
418  /* Up to MAXCOR corrector iterations are taken. A convergence test is */
419  /* made on the RMS-norm of each correction, weighted by the error */
420  /* weight vector EWT. The sum of the corrections is accumulated in the */
421  /* vector ACOR(i). The YH array is not altered in the corrector loop. */
422  /* ----------------------------------------------------------------------- */
423 L220:
424  m = 0;
425  rate = 0.;
426  del = 0.;
427  i__2 = dls001_3.n;
428 
429  for (i__ = 1; i__ <= i__2; ++i__)
430  {
431  /* L230: */
432  y[i__] = yh[i__ + yh_dim1];
433  }
434 
435  f(&neq[1], &dls001_3.tn, &y[1], &savf[1]);
436  ++dls001_3.nfe;
437 
438  if (dls001_3.ipup <= 0)
439  {
440  goto L250;
441  }
442 
443  /* ----------------------------------------------------------------------- */
444  /* If indicated, the matrix P = I - H*EL(1)*J is reevaluated and */
445  /* preprocessed before starting the corrector iteration. IPUP is set */
446  /* to 0 as an indicator that this has been done. */
447  /* ----------------------------------------------------------------------- */
448  (*pjac)(&neq[1], &y[1], &yh[yh_offset], nyh, &ewt[1], &acor[1], &savf[1],
449  &wm[1], &iwm[1], f, jac);
450  dls001_3.ipup = 0;
451  dls001_3.rc = 1.;
452  dls001_3.nslp = dls001_3.nst;
453  dls001_3.crate = .7;
454 
455  if (dls001_3.ierpj != 0)
456  {
457  goto L430;
458  }
459 
460 L250:
461  i__2 = dls001_3.n;
462 
463  for (i__ = 1; i__ <= i__2; ++i__)
464  {
465  /* L260: */
466  acor[i__] = 0.;
467  }
468 
469 L270:
470 
471  if (dls001_3.miter != 0)
472  {
473  goto L350;
474  }
475 
476  /* ----------------------------------------------------------------------- */
477  /* In the case of functional iteration, update Y directly from */
478  /* the result of the last function evaluation. */
479  /* ----------------------------------------------------------------------- */
480  i__2 = dls001_3.n;
481 
482  for (i__ = 1; i__ <= i__2; ++i__)
483  {
484  savf[i__] = dls001_3.h__ * savf[i__] - yh[i__ + (yh_dim1 << 1)];
485  /* L290: */
486  y[i__] = savf[i__] - acor[i__];
487  }
488 
489  del = dmnorm_(&dls001_3.n, &y[1], &ewt[1]);
490  i__2 = dls001_3.n;
491 
492  for (i__ = 1; i__ <= i__2; ++i__)
493  {
494  y[i__] = yh[i__ + yh_dim1] + dls001_3.el[0] * savf[i__];
495  /* L300: */
496  acor[i__] = savf[i__];
497  }
498 
499  goto L400;
500  /* ----------------------------------------------------------------------- */
501  /* In the case of the chord method, compute the corrector error, */
502  /* and solve the linear system with that as right-hand side and */
503  /* P as coefficient matrix. */
504  /* ----------------------------------------------------------------------- */
505 L350:
506  i__2 = dls001_3.n;
507 
508  for (i__ = 1; i__ <= i__2; ++i__)
509  {
510  /* L360: */
511  y[i__] = dls001_3.h__ * savf[i__] - (yh[i__ + (yh_dim1 << 1)] + acor[
512  i__]);
513  }
514 
515  (*slvs)(&wm[1], &iwm[1], &y[1], &savf[1]);
516 
517  if (dls001_3.iersl < 0)
518  {
519  goto L430;
520  }
521 
522  if (dls001_3.iersl > 0)
523  {
524  goto L410;
525  }
526 
527  del = dmnorm_(&dls001_3.n, &y[1], &ewt[1]);
528  i__2 = dls001_3.n;
529 
530  for (i__ = 1; i__ <= i__2; ++i__)
531  {
532  acor[i__] += y[i__];
533  /* L380: */
534  y[i__] = yh[i__ + yh_dim1] + dls001_3.el[0] * acor[i__];
535  }
536 
537  /* ----------------------------------------------------------------------- */
538  /* Test for convergence. If M .gt. 0, an estimate of the convergence */
539  /* rate constant is stored in CRATE, and this is used in the test. */
540 
541  /* We first check for a change of iterates that is the size of */
542  /* roundoff error. If this occurs, the iteration has converged, and a */
543  /* new rate estimate is not formed. */
544  /* In all other cases, force at least two iterations to estimate a */
545  /* local Lipschitz constant estimate for Adams methods. */
546  /* On convergence, form PDEST = local maximum Lipschitz constant */
547  /* estimate. PDLAST is the most recent nonzero estimate. */
548  /* ----------------------------------------------------------------------- */
549 L400:
550 
551  if (del <= pnorm * 100. * dls001_3.uround)
552  {
553  goto L450;
554  }
555 
556  if (m == 0 && dls001_3.meth == 1)
557  {
558  goto L405;
559  }
560 
561  if (m == 0)
562  {
563  goto L402;
564  }
565 
566  rm = 1024.;
567 
568  if (del <= delp * 1024.)
569  {
570  rm = del / delp;
571  }
572 
573  rate = std::max(rate, rm);
574  /* Computing MAX */
575  d__1 = dls001_3.crate * .2;
576  dls001_3.crate = std::max(d__1, rm);
577 L402:
578  /* Computing MIN */
579  d__1 = 1., d__2 = dls001_3.crate * 1.5;
580  dcon = del * std::min(d__1, d__2) / (dls001_3.tesco[dls001_3.nq * 3 - 2] *
581  dls001_3.conit);
582 
583  if (dcon > 1.)
584  {
585  goto L405;
586  }
587 
588  /* Computing MAX */
589  d__2 = dlsa01_1.pdest, d__3 = rate / (d__1 = dls001_3.h__ * dls001_3.el[0]
590  , fabs(d__1));
591  dlsa01_1.pdest = std::max(d__2, d__3);
592 
593  if (dlsa01_1.pdest != 0.)
594  {
595  dlsa01_1.pdlast = dlsa01_1.pdest;
596  }
597 
598  goto L450;
599 L405:
600  ++m;
601 
602  if (m == dls001_3.maxcor)
603  {
604  goto L410;
605  }
606 
607  if (m >= 2 && del > delp * 2.)
608  {
609  goto L410;
610  }
611 
612  delp = del;
613  f(&neq[1], &dls001_3.tn, &y[1], &savf[1]);
614  ++dls001_3.nfe;
615  goto L270;
616  /* ----------------------------------------------------------------------- */
617  /* The corrector iteration failed to converge. */
618  /* If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for */
619  /* the next try. Otherwise the YH array is retracted to its values */
620  /* before prediction, and H is reduced, if possible. If H cannot be */
621  /* reduced or MXNCF failures have occurred, exit with KFLAG = -2. */
622  /* ----------------------------------------------------------------------- */
623 L410:
624 
625  if (dls001_3.miter == 0 || dls001_3.jcur == 1)
626  {
627  goto L430;
628  }
629 
630  dls001_3.icf = 1;
631  dls001_3.ipup = dls001_3.miter;
632  goto L220;
633 L430:
634  dls001_3.icf = 2;
635  ++ncf;
636  dls001_3.rmax = 2.;
637  dls001_3.tn = told;
638  i1 = dls001_3.nqnyh + 1;
639  i__2 = dls001_3.nq;
640 
641  for (jb = 1; jb <= i__2; ++jb)
642  {
643  i1 -= *nyh;
644  /* DIR$ IVDEP */
645  i__1 = dls001_3.nqnyh;
646 
647  for (i__ = i1; i__ <= i__1; ++i__)
648  {
649  /* L440: */
650  yh1[i__] -= yh1[i__ + *nyh];
651  }
652 
653  /* L445: */
654  }
655 
656  if (dls001_3.ierpj < 0 || dls001_3.iersl < 0)
657  {
658  goto L680;
659  }
660 
661  if (fabs(dls001_3.h__) <= dls001_3.hmin * 1.00001)
662  {
663  goto L670;
664  }
665 
666  if (ncf == dls001_3.mxncf)
667  {
668  goto L670;
669  }
670 
671  rh = .25;
672  dls001_3.ipup = dls001_3.miter;
673  iredo = 1;
674  goto L170;
675  /* ----------------------------------------------------------------------- */
676  /* The corrector has converged. JCUR is set to 0 */
677  /* to signal that the Jacobian involved may need updating later. */
678  /* The local error test is made and control passes to statement 500 */
679  /* if it fails. */
680  /* ----------------------------------------------------------------------- */
681 L450:
682  dls001_3.jcur = 0;
683 
684  if (m == 0)
685  {
686  dsm = del / dls001_3.tesco[dls001_3.nq * 3 - 2];
687  }
688 
689  if (m > 0)
690  {
691  dsm = dmnorm_(&dls001_3.n, &acor[1], &ewt[1]) / dls001_3.tesco[
692  dls001_3.nq * 3 - 2];
693  }
694 
695  if (dsm > 1.)
696  {
697  goto L500;
698  }
699 
700  /* ----------------------------------------------------------------------- */
701  /* After a successful step, update the YH array. */
702  /* Decrease ICOUNT by 1, and if it is -1, consider switching methods. */
703  /* If a method switch is made, reset various parameters, */
704  /* rescale the YH array, and exit. If there is no switch, */
705  /* consider changing H if IALTH = 1. Otherwise decrease IALTH by 1. */
706  /* If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for */
707  /* use in a possible order increase on the next step. */
708  /* If a change in H is considered, an increase or decrease in order */
709  /* by one is considered also. A change in H is made only if it is by a */
710  /* factor of at least 1.1. If not, IALTH is set to 3 to prevent */
711  /* testing for that many steps. */
712  /* ----------------------------------------------------------------------- */
713  dls001_3.kflag = 0;
714  iredo = 0;
715  ++dls001_3.nst;
716  dls001_3.hu = dls001_3.h__;
717  dls001_3.nqu = dls001_3.nq;
718  dlsa01_1.mused = dls001_3.meth;
719  i__2 = dls001_3.l;
720 
721  for (j = 1; j <= i__2; ++j)
722  {
723  i__1 = dls001_3.n;
724 
725  for (i__ = 1; i__ <= i__1; ++i__)
726  {
727  /* L460: */
728  yh[i__ + j * yh_dim1] += dls001_3.el[j - 1] * acor[i__];
729  }
730  }
731 
732  --dlsa01_1.icount;
733 
734  if (dlsa01_1.icount >= 0)
735  {
736  goto L488;
737  }
738 
739  if (dls001_3.meth == 2)
740  {
741  goto L480;
742  }
743 
744  /* ----------------------------------------------------------------------- */
745  /* We are currently using an Adams method. Consider switching to BDF. */
746  /* If the current order is greater than 5, assume the problem is */
747  /* not stiff, and skip this section. */
748  /* If the Lipschitz constant and error estimate are not polluted */
749  /* by roundoff, go to 470 and perform the usual test. */
750  /* Otherwise, switch to the BDF methods if the last step was */
751  /* restricted to insure stability (irflag = 1), and stay with Adams */
752  /* method if not. When switching to BDF with polluted error estimates, */
753  /* in the absence of other information, double the step size. */
754 
755  /* When the estimates are OK, we make the usual test by computing */
756 
757  /* the step size we could have (ideally) used on this step, */
758 
759  /* with the current (Adams) method, and also that for the BDF. */
760 
761  /* If NQ .gt. MXORDS, we consider changing to order MXORDS on switching. */
762 
763  /* Compare the two step sizes to decide whether to switch. */
764 
765  /* The step size advantage must be at least RATIO = 5 to switch. */
766 
767  /* ----------------------------------------------------------------------- */
768  if (dls001_3.nq > 5)
769  {
770  goto L488;
771  }
772 
773  if (dsm > pnorm * 100. * dls001_3.uround && dlsa01_1.pdest != 0.)
774  {
775  goto L470;
776  }
777 
778  if (dlsa01_1.irflag == 0)
779  {
780  goto L488;
781  }
782 
783  rh2 = 2.;
784  nqm2 = std::min(dls001_3.nq, dlsa01_1.mxords);
785  goto L478;
786 L470:
787  exsm = 1. / dls001_3.l;
788  rh1 = 1. / (pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
789  rh1it = rh1 * 2.;
790  pdh = dlsa01_1.pdlast * fabs(dls001_3.h__);
791 
792  if (pdh * rh1 > 1e-5)
793  {
794  rh1it = sm1[dls001_3.nq - 1] / pdh;
795  }
796 
797  rh1 = std::min(rh1, rh1it);
798 
799  if (dls001_3.nq <= dlsa01_1.mxords)
800  {
801  goto L474;
802  }
803 
804  nqm2 = dlsa01_1.mxords;
805  lm2 = dlsa01_1.mxords + 1;
806  exm2 = 1. / lm2;
807  lm2p1 = lm2 + 1;
808  dm2 = dmnorm_(&dls001_3.n, &yh[lm2p1 * yh_dim1 + 1], &ewt[1]) /
809  dlsa01_1.cm2[dlsa01_1.mxords - 1];
810  rh2 = 1. / (pow_dd(&dm2, &exm2) * 1.2 + 1.2e-6);
811  goto L476;
812 L474:
813  dm2 = dsm * (dlsa01_1.cm1[dls001_3.nq - 1] / dlsa01_1.cm2[dls001_3.nq - 1]
814  );
815  rh2 = 1. / (pow_dd(&dm2, &exsm) * 1.2 + 1.2e-6);
816  nqm2 = dls001_3.nq;
817 L476:
818 
819  if (rh2 < dlsa01_1.ratio * rh1)
820  {
821  goto L488;
822  }
823 
824  /* THE SWITCH TEST PASSED. RESET RELEVANT QUANTITIES FOR BDF. ---------- */
825 L478:
826  rh = rh2;
827  dlsa01_1.icount = 20;
828  dls001_3.meth = 2;
829  dls001_3.miter = dlsa01_1.jtyp;
830  dlsa01_1.pdlast = 0.;
831  dls001_3.nq = nqm2;
832  dls001_3.l = dls001_3.nq + 1;
833  goto L170;
834  /* ----------------------------------------------------------------------- */
835  /* We are currently using a BDF method. Consider switching to Adams. */
836  /* Compute the step size we could have (ideally) used on this step, */
837  /* with the current (BDF) method, and also that for the Adams. */
838  /* If NQ .gt. MXORDN, we consider changing to order MXORDN on switching. */
839  /* Compare the two step sizes to decide whether to switch. */
840  /* The step size advantage must be at least 5/RATIO = 1 to switch. */
841  /* If the step size for Adams would be so small as to cause */
842  /* roundoff pollution, we stay with BDF. */
843  /* ----------------------------------------------------------------------- */
844 L480:
845  exsm = 1. / dls001_3.l;
846 
847  if (dlsa01_1.mxordn >= dls001_3.nq)
848  {
849  goto L484;
850  }
851 
852  nqm1 = dlsa01_1.mxordn;
853  lm1 = dlsa01_1.mxordn + 1;
854  exm1 = 1. / lm1;
855  lm1p1 = lm1 + 1;
856  dm1 = dmnorm_(&dls001_3.n, &yh[lm1p1 * yh_dim1 + 1], &ewt[1]) /
857  dlsa01_1.cm1[dlsa01_1.mxordn - 1];
858  rh1 = 1. / (pow_dd(&dm1, &exm1) * 1.2 + 1.2e-6);
859  goto L486;
860 L484:
861  dm1 = dsm * (dlsa01_1.cm2[dls001_3.nq - 1] / dlsa01_1.cm1[dls001_3.nq - 1]
862  );
863  rh1 = 1. / (pow_dd(&dm1, &exsm) * 1.2 + 1.2e-6);
864  nqm1 = dls001_3.nq;
865  exm1 = exsm;
866 L486:
867  rh1it = rh1 * 2.;
868  pdh = dlsa01_1.pdnorm * fabs(dls001_3.h__);
869 
870  if (pdh * rh1 > 1e-5)
871  {
872  rh1it = sm1[nqm1 - 1] / pdh;
873  }
874 
875  rh1 = std::min(rh1, rh1it);
876  rh2 = 1. / (pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
877 
878  if (rh1 * dlsa01_1.ratio < rh2 * 5.)
879  {
880  goto L488;
881  }
882 
883  alpha = std::max(.001, rh1);
884  dm1 = pow_dd(&alpha, &exm1) * dm1;
885 
886  if (dm1 <= dls001_3.uround * 1e3 * pnorm)
887  {
888  goto L488;
889  }
890 
891  /* The switch test passed. Reset relevant quantities for Adams. -------- */
892  rh = rh1;
893  dlsa01_1.icount = 20;
894  dls001_3.meth = 1;
895  dls001_3.miter = 0;
896  dlsa01_1.pdlast = 0.;
897  dls001_3.nq = nqm1;
898  dls001_3.l = dls001_3.nq + 1;
899  goto L170;
900 
901  /* No method switch is being made. Do the usual step/order selection. -- */
902 L488:
903  --dls001_3.ialth;
904 
905  if (dls001_3.ialth == 0)
906  {
907  goto L520;
908  }
909 
910  if (dls001_3.ialth > 1)
911  {
912  goto L700;
913  }
914 
915  if (dls001_3.l == dls001_3.lmax)
916  {
917  goto L700;
918  }
919 
920  i__1 = dls001_3.n;
921 
922  for (i__ = 1; i__ <= i__1; ++i__)
923  {
924  /* L490: */
925  yh[i__ + dls001_3.lmax * yh_dim1] = acor[i__];
926  }
927 
928  goto L700;
929  /* ----------------------------------------------------------------------- */
930  /* The error test failed. KFLAG keeps track of multiple failures. */
931  /* Restore TN and the YH array to their previous values, and prepare */
932  /* to try the step again. Compute the optimum step size for this or */
933  /* one lower order. After 2 or more failures, H is forced to decrease */
934  /* by a factor of 0.2 or less. */
935  /* ----------------------------------------------------------------------- */
936 L500:
937  --dls001_3.kflag;
938  dls001_3.tn = told;
939  i1 = dls001_3.nqnyh + 1;
940  i__1 = dls001_3.nq;
941 
942  for (jb = 1; jb <= i__1; ++jb)
943  {
944  i1 -= *nyh;
945  /* DIR$ IVDEP */
946  i__2 = dls001_3.nqnyh;
947 
948  for (i__ = i1; i__ <= i__2; ++i__)
949  {
950  /* L510: */
951  yh1[i__] -= yh1[i__ + *nyh];
952  }
953 
954  /* L515: */
955  }
956 
957  dls001_3.rmax = 2.;
958 
959  if (fabs(dls001_3.h__) <= dls001_3.hmin * 1.00001)
960  {
961  goto L660;
962  }
963 
964  if (dls001_3.kflag <= -3)
965  {
966  goto L640;
967  }
968 
969  iredo = 2;
970  rhup = 0.;
971  goto L540;
972  /* ----------------------------------------------------------------------- */
973  /* Regardless of the success or failure of the step, factors */
974  /* RHDN, RHSM, and RHUP are computed, by which H could be multiplied */
975  /* at order NQ - 1, order NQ, or order NQ + 1, respectively. */
976  /* In the case of failure, RHUP = 0.0 to avoid an order increase. */
977  /* The largest of these is determined and the new order chosen */
978  /* accordingly. If the order is to be increased, we compute one */
979  /* additional scaled derivative. */
980  /* ----------------------------------------------------------------------- */
981 L520:
982  rhup = 0.;
983 
984  if (dls001_3.l == dls001_3.lmax)
985  {
986  goto L540;
987  }
988 
989  i__1 = dls001_3.n;
990 
991  for (i__ = 1; i__ <= i__1; ++i__)
992  {
993  /* L530: */
994  savf[i__] = acor[i__] - yh[i__ + dls001_3.lmax * yh_dim1];
995  }
996 
997  dup = dmnorm_(&dls001_3.n, &savf[1], &ewt[1]) / dls001_3.tesco[
998  dls001_3.nq * 3 - 1];
999  exup = 1. / (dls001_3.l + 1);
1000  rhup = 1. / (pow_dd(&dup, &exup) * 1.4 + 1.4e-6);
1001 L540:
1002  exsm = 1. / dls001_3.l;
1003  rhsm = 1. / (pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
1004  rhdn = 0.;
1005 
1006  if (dls001_3.nq == 1)
1007  {
1008  goto L550;
1009  }
1010 
1011  ddn = dmnorm_(&dls001_3.n, &yh[dls001_3.l * yh_dim1 + 1], &ewt[1]) /
1012  dls001_3.tesco[dls001_3.nq * 3 - 3];
1013  exdn = 1. / dls001_3.nq;
1014  rhdn = 1. / (pow_dd(&ddn, &exdn) * 1.3 + 1.3e-6);
1015  /* If METH = 1, limit RH according to the stability region also. -------- */
1016 L550:
1017 
1018  if (dls001_3.meth == 2)
1019  {
1020  goto L560;
1021  }
1022 
1023  /* Computing MAX */
1024  d__1 = fabs(dls001_3.h__) * dlsa01_1.pdlast;
1025  pdh = std::max(d__1, 1e-6);
1026 
1027  if (dls001_3.l < dls001_3.lmax)
1028  {
1029  /* Computing MIN */
1030  d__1 = rhup, d__2 = sm1[dls001_3.l - 1] / pdh;
1031  rhup = std::min(d__1, d__2);
1032  }
1033 
1034  /* Computing MIN */
1035  d__1 = rhsm, d__2 = sm1[dls001_3.nq - 1] / pdh;
1036  rhsm = std::min(d__1, d__2);
1037 
1038  if (dls001_3.nq > 1)
1039  {
1040  /* Computing MIN */
1041  d__1 = rhdn, d__2 = sm1[dls001_3.nq - 2] / pdh;
1042  rhdn = std::min(d__1, d__2);
1043  }
1044 
1045  dlsa01_1.pdest = 0.;
1046 L560:
1047 
1048  if (rhsm >= rhup)
1049  {
1050  goto L570;
1051  }
1052 
1053  if (rhup > rhdn)
1054  {
1055  goto L590;
1056  }
1057 
1058  goto L580;
1059 L570:
1060 
1061  if (rhsm < rhdn)
1062  {
1063  goto L580;
1064  }
1065 
1066  newq = dls001_3.nq;
1067  rh = rhsm;
1068  goto L620;
1069 L580:
1070  newq = dls001_3.nq - 1;
1071  rh = rhdn;
1072 
1073  if (dls001_3.kflag < 0 && rh > 1.)
1074  {
1075  rh = 1.;
1076  }
1077 
1078  goto L620;
1079 L590:
1080  newq = dls001_3.l;
1081  rh = rhup;
1082 
1083  if (rh < 1.1)
1084  {
1085  goto L610;
1086  }
1087 
1088  r__ = dls001_3.el[dls001_3.l - 1] / dls001_3.l;
1089  i__1 = dls001_3.n;
1090 
1091  for (i__ = 1; i__ <= i__1; ++i__)
1092  {
1093  /* L600: */
1094  yh[i__ + (newq + 1) * yh_dim1] = acor[i__] * r__;
1095  }
1096 
1097  goto L630;
1098 L610:
1099  dls001_3.ialth = 3;
1100  goto L700;
1101  /* If METH = 1 and H is restricted by stability, bypass 10 percent test. */
1102 L620:
1103 
1104  if (dls001_3.meth == 2)
1105  {
1106  goto L622;
1107  }
1108 
1109  if (rh * pdh * 1.00001 >= sm1[newq - 1])
1110  {
1111  goto L625;
1112  }
1113 
1114 L622:
1115 
1116  if (dls001_3.kflag == 0 && rh < 1.1)
1117  {
1118  goto L610;
1119  }
1120 
1121 L625:
1122 
1123  if (dls001_3.kflag <= -2)
1124  {
1125  rh = std::min(rh, .2);
1126  }
1127 
1128  /* ----------------------------------------------------------------------- */
1129 
1130  /* If there is a change of order, reset NQ, L, and the coefficients. */
1131 
1132  /* In any case H is reset according to RH and the YH array is rescaled. */
1133 
1134  /* Then exit from 690 if the step was OK, or redo the step otherwise. */
1135 
1136  /* ----------------------------------------------------------------------- */
1137  if (newq == dls001_3.nq)
1138  {
1139  goto L170;
1140  }
1141 
1142 L630:
1143  dls001_3.nq = newq;
1144  dls001_3.l = dls001_3.nq + 1;
1145  iret = 2;
1146  goto L150;
1147  /* ----------------------------------------------------------------------- */
1148  /* Control reaches this section if 3 or more failures have occured. */
1149  /* If 10 failures have occurred, exit with KFLAG = -1. */
1150  /* It is assumed that the derivatives that have accumulated in the */
1151  /* YH array have errors of the wrong order. Hence the first */
1152  /* derivative is recomputed, and the order is set to 1. Then */
1153  /* H is reduced by a factor of 10, and the step is retried, */
1154  /* until it succeeds or H reaches HMIN. */
1155  /* ----------------------------------------------------------------------- */
1156 L640:
1157 
1158  if (dls001_3.kflag == -10)
1159  {
1160  goto L660;
1161  }
1162 
1163  rh = .1;
1164  /* Computing MAX */
1165  d__1 = dls001_3.hmin / fabs(dls001_3.h__);
1166  rh = std::max(d__1, rh);
1167  dls001_3.h__ *= rh;
1168  i__1 = dls001_3.n;
1169 
1170  for (i__ = 1; i__ <= i__1; ++i__)
1171  {
1172  /* L645: */
1173  y[i__] = yh[i__ + yh_dim1];
1174  }
1175 
1176  f(&neq[1], &dls001_3.tn, &y[1], &savf[1]);
1177  ++dls001_3.nfe;
1178  i__1 = dls001_3.n;
1179 
1180  for (i__ = 1; i__ <= i__1; ++i__)
1181  {
1182  /* L650: */
1183  yh[i__ + (yh_dim1 << 1)] = dls001_3.h__ * savf[i__];
1184  }
1185 
1186  dls001_3.ipup = dls001_3.miter;
1187  dls001_3.ialth = 5;
1188 
1189  if (dls001_3.nq == 1)
1190  {
1191  goto L200;
1192  }
1193 
1194  dls001_3.nq = 1;
1195  dls001_3.l = 2;
1196  iret = 3;
1197  goto L150;
1198  /* ----------------------------------------------------------------------- */
1199  /* All returns are made through this section. H is saved in HOLD */
1200  /* to allow the caller to change H on the next step. */
1201  /* ----------------------------------------------------------------------- */
1202 L660:
1203  dls001_3.kflag = -1;
1204  goto L720;
1205 L670:
1206  dls001_3.kflag = -2;
1207  goto L720;
1208 L680:
1209  dls001_3.kflag = -3;
1210  goto L720;
1211 L690:
1212  dls001_3.rmax = 10.;
1213 L700:
1214  r__ = 1. / dls001_3.tesco[dls001_3.nqu * 3 - 2];
1215  i__1 = dls001_3.n;
1216 
1217  for (i__ = 1; i__ <= i__1; ++i__)
1218  {
1219  /* L710: */
1220  acor[i__] *= r__;
1221  }
1222 
1223 L720:
1224  dls001_3.hold = dls001_3.h__;
1225  dls001_3.jstart = 1;
1226  return 0;
1227  /* ----------------------- End of Subroutine DSTODA ---------------------- */
1228 } /* dstoda_ */
#define C_INT
Definition: copasi.h:115
#define pow_dd(__x, __y)
Definition: dstoda.cpp:38
static C_INT c__2
Definition: dstoda.cpp:36
C_INT dcfode_(C_INT *meth, double *elco, double *tesco)
Definition: dcfode.cpp:22
static C_INT c__1
Definition: dstoda.cpp:35
#define dls001_3
Definition: dstoda.cpp:29
#define dlsa01_1
Definition: dstoda.cpp:31
double dmnorm_(C_INT *n, double *v, double *w)
Definition: dmnorm.cpp:35
#define min(a, b)
Definition: f2c.h:175
#define max(a, b)
Definition: f2c.h:176
void CInternalSolver::enablePrint ( const bool &  print = true)

Definition at line 27 of file CInternalSolver.cpp.

References Cxerrwd::enablePrint(), and mxerrwd.

28 {
29  mxerrwd.enablePrint(print);
30 }
void enablePrint(const bool &print=true)
Definition: Cxerrwd.cpp:35
void CInternalSolver::resetState ( )
void CInternalSolver::saveState ( )
void CInternalSolver::setOstream ( std::ostream &  os)

Member Data Documentation

dls001 CInternalSolver::mdls001_
protected

Definition at line 74 of file CInternalSolver.h.

Referenced by resetState(), and saveState().

dlsa01 CInternalSolver::mdlsa01_
protected

Definition at line 75 of file CInternalSolver.h.

Referenced by resetState(), and saveState().

dlsr01 CInternalSolver::mdlsr01_
protected

Definition at line 76 of file CInternalSolver.h.

Referenced by resetState(), and saveState().

State CInternalSolver::mState
protected

Definition at line 77 of file CInternalSolver.h.

Referenced by resetState(), and saveState().

Cxerrwd CInternalSolver::mxerrwd
protected

The documentation for this class was generated from the following files: