COPASI API  4.16.103
Functions
dcfode.cpp File Reference
#include "copasi.h"
#include "dcfode.h"
Include dependency graph for dcfode.cpp:

Go to the source code of this file.

Functions

C_INT dcfode_ (C_INT *meth, double *elco, double *tesco)
 

Function Documentation

C_INT dcfode_ ( C_INT meth,
double *  elco,
double *  tesco 
)

Definition at line 22 of file dcfode.cpp.

References C_INT.

Referenced by CInternalSolver::dstoda_().

23 {
24  /* System generated locals */
25  C_INT i__1;
26 
27  /* Local variables */
28  C_INT i__, ib;
29  double pc[12];
30  C_INT nq;
31  double fnq;
32  C_INT nqm1, nqp1;
33  double ragq, pint, xpin, fnqm1, agamq, rqfac, tsign, rq1fac;
34 
35  /* ***BEGIN PROLOGUE DCFODE */
36  /* ***SUBSIDIARY */
37  /* ***PURPOSE Set ODE integrator coefficients. */
38  /* ***TYPE DOUBLE PRECISION (SCFODE-S, DCFODE-D) */
39  /* ***AUTHOR Hindmarsh, Alan C., (LLNL) */
40  /* ***DESCRIPTION */
41 
42  /* DCFODE is called by the integrator routine to set coefficients */
43  /* needed there. The coefficients for the current method, as */
44  /* given by the value of METH, are set for all orders and saved. */
45  /* The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2. */
46  /* (A smaller value of the maximum order is also allowed.) */
47  /* DCFODE is called once at the beginning of the problem, */
48  /* and is not called again unless and until METH is changed. */
49 
50  /* The ELCO array contains the basic method coefficients. */
51  /* The coefficients el(i), 1 .le. i .le. nq+1, for the method of */
52  /* order nq are stored in ELCO(i,nq). They are given by a genetrating */
53  /* polynomial, i.e., */
54  /* l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq. */
55  /* For the implicit Adams methods, l(x) is given by */
56  /* dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0. */
57  /* For the BDF methods, l(x) is given by */
58  /* l(x) = (x+1)*(x+2)* ... *(x+nq)/K, */
59  /* where K = factorial(nq)*(1 + 1/2 + ... + 1/nq). */
60 
61  /* The TESCO array contains test constants used for the */
62  /* local error test and the selection of step size and/or order. */
63  /* At order nq, TESCO(k,nq) is used for the selection of step */
64  /* size at order nq - 1 if k = 1, at order nq if k = 2, and at order */
65  /* nq + 1 if k = 3. */
66 
67  /* ***SEE ALSO DLSODE */
68  /* ***ROUTINES CALLED (NONE) */
69  /* ***REVISION HISTORY (YYMMDD) */
70  /* 791129 DATE WRITTEN */
71  /* 890501 Modified prologue to SLATEC/LDOC format. (FNF) */
72  /* 890503 Minor cosmetic changes. (FNF) */
73  /* 930809 Renamed to allow single/double precision versions. (ACH) */
74  /* ***END PROLOGUE DCFODE */
75  /* **End */
76 
77  /* ***FIRST EXECUTABLE STATEMENT DCFODE */
78  /* Parameter adjustments */
79  tesco -= 4;
80  elco -= 14;
81 
82  /* Function Body */
83  switch (*meth)
84  {
85  case 1: goto L100;
86  case 2: goto L200;
87  }
88 
89 L100:
90  elco[14] = 1.;
91  elco[15] = 1.;
92  tesco[4] = 0.;
93  tesco[5] = 2.;
94  tesco[7] = 1.;
95  tesco[39] = 0.;
96  pc[0] = 1.;
97  rqfac = 1.;
98  for (nq = 2; nq <= 12; ++nq)
99  {
100  /* ----------------------------------------------------------------------- */
101  /* The PC array will contain the coefficients of the polynomial */
102  /* p(x) = (x+1)*(x+2)*...*(x+nq-1). */
103  /* Initially, p(x) = 1. */
104  /* ----------------------------------------------------------------------- */
105  rq1fac = rqfac;
106  rqfac /= nq;
107  nqm1 = nq - 1;
108  fnqm1 = (double) nqm1;
109  nqp1 = nq + 1;
110  /* Form coefficients of p(x)*(x+nq-1). ---------------------------------- */
111  pc[nq - 1] = 0.;
112  i__1 = nqm1;
113  for (ib = 1; ib <= i__1; ++ib)
114  {
115  i__ = nqp1 - ib;
116  /* L110: */
117  pc[i__ - 1] = pc[i__ - 2] + fnqm1 * pc[i__ - 1];
118  }
119  pc[0] = fnqm1 * pc[0];
120  /* Compute integral, -1 to 0, of p(x) and x*p(x). ----------------------- */
121  pint = pc[0];
122  xpin = pc[0] / 2.;
123  tsign = 1.;
124  i__1 = nq;
125  for (i__ = 2; i__ <= i__1; ++i__)
126  {
127  tsign = -tsign;
128  pint += tsign * pc[i__ - 1] / i__;
129  /* L120: */
130  xpin += tsign * pc[i__ - 1] / (i__ + 1);
131  }
132  /* Store coefficients in ELCO and TESCO. -------------------------------- */
133  elco[nq * 13 + 1] = pint * rq1fac;
134  elco[nq * 13 + 2] = 1.;
135  i__1 = nq;
136  for (i__ = 2; i__ <= i__1; ++i__)
137  {
138  /* L130: */
139  elco[i__ + 1 + nq * 13] = rq1fac * pc[i__ - 1] / i__;
140  }
141  agamq = rqfac * xpin;
142  ragq = 1. / agamq;
143  tesco[nq * 3 + 2] = ragq;
144  if (nq < 12)
145  {
146  tesco[nqp1 * 3 + 1] = ragq * rqfac / nqp1;
147  }
148  tesco[nqm1 * 3 + 3] = ragq;
149  /* L140: */
150  }
151  return 0;
152 
153 L200:
154  pc[0] = 1.;
155  rq1fac = 1.;
156  for (nq = 1; nq <= 5; ++nq)
157  {
158  /* ----------------------------------------------------------------------- */
159  /* The PC array will contain the coefficients of the polynomial */
160  /* p(x) = (x+1)*(x+2)*...*(x+nq). */
161  /* Initially, p(x) = 1. */
162  /* ----------------------------------------------------------------------- */
163  fnq = (double) nq;
164  nqp1 = nq + 1;
165  /* Form coefficients of p(x)*(x+nq). ------------------------------------ */
166  pc[nqp1 - 1] = 0.;
167  i__1 = nq;
168  for (ib = 1; ib <= i__1; ++ib)
169  {
170  i__ = nq + 2 - ib;
171  /* L210: */
172  pc[i__ - 1] = pc[i__ - 2] + fnq * pc[i__ - 1];
173  }
174  pc[0] = fnq * pc[0];
175  /* Store coefficients in ELCO and TESCO. -------------------------------- */
176  i__1 = nqp1;
177  for (i__ = 1; i__ <= i__1; ++i__)
178  {
179  /* L220: */
180  elco[i__ + nq * 13] = pc[i__ - 1] / pc[1];
181  }
182  elco[nq * 13 + 2] = 1.;
183  tesco[nq * 3 + 1] = rq1fac;
184  tesco[nq * 3 + 2] = nqp1 / elco[nq * 13 + 1];
185  tesco[nq * 3 + 3] = (nq + 2) / elco[nq * 13 + 1];
186  rq1fac /= fnq;
187  /* L230: */
188  }
189  return 0;
190  /* ----------------------- END OF SUBROUTINE DCFODE ---------------------- */
191 } /* dcfode_ */
#define C_INT
Definition: copasi.h:115