COPASI API  4.16.103
dintdy.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 #include <string>
21 
22 #include <algorithm>
23 
24 #include "copasi.h"
25 
26 #include "CInternalSolver.h"
27 #include "Cxerrwd.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 double d_sign(const double & a, const double & b);
38 
39 static double c_b34 = 0.0;
40 
41 static C_INT c__0 = 0;
42 static C_INT c__1 = 1;
43 static C_INT c__2 = 2;
44 static C_INT c__30 = 30;
45 static C_INT c__51 = 51;
46 static C_INT c__52 = 52;
47 static C_INT c__60 = 60;
48 
49 #define pow_di(__x, __n) pow(*__x, (double) *__n)
50 
51 /* DECK DINTDY */
52 /* Subroutine */
53 C_INT CInternalSolver::dintdy_(double *t, const C_INT *k, double *yh,
54  C_INT *nyh, double *dky, C_INT *iflag)
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
C_INT dintdy_(double *t, const C_INT *k, double *yh, C_INT *nyh, double *dky, C_INT *iflag)
Definition: dintdy.cpp:53
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