COPASI API  4.16.103
dprja.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2014 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2006 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 //
16 // This C++ code is based on an f2c conversion of the Fortran
17 // library ODEPACK available at: http://www.netlib.org/odepack/
18 
19 #include <cmath>
20 
21 #include <algorithm>
22 
23 #include "copasi.h"
24 
25 #include "CInternalSolver.h"
26 
27 #include "lapack/lapackwrap.h"
28 
29 #define dls001_1 (mdls001_._1)
30 #define dls001_2 (mdls001_._2)
31 #define dls001_3 (mdls001_._3)
32 
33 #define dlsa01_1 (mdlsa01_._1)
34 #define dlsa01_2 (mdlsa01_._2)
35 #define dlsa01_3 (mdlsa01_._3)
36 
37 static C_INT c__0 = 0;
38 
39 #include "dgbfa.h"
40 #include "dgefa.h"
41 #include "dbnorm.h"
42 #include "dfnorm.h"
43 #include "dmnorm.h"
44 
45 /* DECK DPRJA */
46 /* Subroutine */
47 C_INT CInternalSolver::dprja_(C_INT *neq, double *y, double *yh,
48  C_INT *nyh, double *ewt, double *ftem, double *savf,
49  double *wm, C_INT *iwm, evalF f, evalJ jac)
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
void(* evalF)(const C_INT *, const double *, const double *, double *)
Definition: common.h:29
double dbnorm_(C_INT *n, double *a, C_INT *nra, C_INT *ml, C_INT *mu, double *w)
Definition: dbnorm.cpp:31
void(* evalJ)(const C_INT *, const double *, const double *, const C_INT *, const C_INT *, double *, const C_INT *)
Definition: common.h:30
#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
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)
Definition: dprja.cpp:47
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