COPASI API  4.16.103
dsolsy.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 <algorithm>
21 
22 #include "copasi.h"
23 
24 #include "CInternalSolver.h"
25 
26 #define dls001_1 (mdls001_._1)
27 #define dls001_2 (mdls001_._2)
28 #define dls001_3 (mdls001_._3)
29 
30 #define dlsa01_1 (mdlsa01_._1)
31 #define dlsa01_2 (mdlsa01_._2)
32 #define dlsa01_3 (mdlsa01_._3)
33 
34 static C_INT c__0 = 0;
35 
36 #include "dgbsl.h"
37 #include "dgesl.h"
38 
39 #include "lapack/lapackwrap.h"
40 
41 /* DECK DSOLSY */
42 /* Subroutine */
43 C_INT CInternalSolver::dsolsy_(double *wm, C_INT *iwm, double *x,
44  double *tem)
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
C_INT dsolsy_(double *wm, C_INT *iwm, double *x, double *tem)
Definition: dsolsy.cpp:43
#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