home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / c_lsode.zip / LSODE.H < prev    next >
Text File  |  1994-01-15  |  22KB  |  305 lines
  1.  
  2.  
  3. /* --------------------------------------------------------------------- */
  4. /*  this is the may 18, 1993 version of clsode. (c version of lsode)     */
  5. /*  this version is 100% ansi c compatible.                              */
  6. /*                                                                       */
  7. /*  for more information about clsode or mathpack++,                     */
  8. /*     contact..   Woon-Chul Choi                                        */
  9. /*                 choi2@chd1s0.engr.ccny.cuny.edu                       */
  10. /*                                                                       */
  11. /*  clsode was translated and modified from hindmarsh's lsode written    */
  12. /*  in fortran.                                                          */
  13. /*                                                                       */
  14. /*  I did my best to avoid and, later, to find a bug, but there must     */
  15. /*  be bugs as always.  if you ever find any, please let me know.        */
  16. /*                                                                       */
  17. /* --------------------------------------------------------------------- */
  18. /*  file list                                                            */
  19. /*      ls0001.h  - header file for fortran common block ls0001          */
  20. /*      lsode.h   - main header for clsode                               */
  21. /*      cfode.c   - c version of cfode                                   */
  22. /*      cutil.c   - utility routine module for odepack                   */
  23. /*      d1mach.c  - c version of d1mach                                  */
  24. /*      daxpy.c   - c version of daxpy                                   */
  25. /*      ddot.c    - c version of ddot                                    */
  26. /*      dgbfa.c   - c version of dgbfa                                   */
  27. /*      dgbsl.c   - c version of dgbsl                                   */
  28. /*      dgefa.c   - c version of dgefa                                   */
  29. /*      dgesl.c   - c version of dgesl                                   */
  30. /*      dscal.c   - c version of dscal                                   */
  31. /*      ewset.c   - c version of ewset                                   */
  32. /*      idamax.c  - c version of idamax                                  */
  33. /*      intdy.c   - c version of intdy                                   */
  34. /*      lsode.c   - c version of lsode                                   */
  35. /*      prepj.c   - c version of prepj                                   */
  36. /*      solsy.c   - c version of solsy                                   */
  37. /*      stode.c   - c version of stode                                   */
  38. /*      vnorm.c   - c version of vnorm                                   */
  39. /*      Xermsg.c  - message handler for mathpack++                       */
  40. /*                                                                       */
  41. /* --------------------------------------------------------------------- */
  42. /*  argument transfer in clsode                                          */
  43. /*      following the type of the argument transfer (i.e. a scalar to a  */
  44. /*      vector array, or a vector to a matrix array) is not unusual in   */
  45. /*      fortran version of lsode.  but it is not allowed in c language   */
  46. /*      as far as i know. thus i applied the array arithmatic to solve   */
  47. /*      the problem. (in ansi fortran 77, a scalar is not the same as    */
  48. /*      an array of length 1, too.)                                      */
  49. /*                                                                       */
  50. /*      example :                                                        */
  51. /*                                                                       */
  52. /*         subroutine lsode(rwork, ....)                                 */
  53. /*         dimension rwork(1)                                            */
  54. /*         :                                                             */
  55. /*         ldf = 21                                                      */
  56. /*         nr = 2                                                        */
  57. /*         call stode(rwork(ldf), nr,...)  <-- a vector array in         */
  58. /*         :                                                             */
  59. /*         :                                                             */
  60. /*         end                                                           */
  61. /*                                                                       */
  62. /*         subroutine stode(yh, nr,.....)                                */
  63. /*         dimension yh(nr,1)              <-- become a matrix array     */
  64. /*         :                                                             */
  65. /*         :                                                             */
  66. /*         end                                                           */
  67. /*                                                                       */
  68. /*      in fortran, argument connection is following                     */
  69. /*         yh(1,1) = rwork(21)                                           */
  70. /*         yh(2,1) = rwork(22)                                           */
  71. /*         yh(1,2) = rwork(23)                                           */
  72. /*         yh(2,2) = rwork(24)                                           */
  73. /*                                                                       */
  74. /*      in c, using array arithmatic: yh(i,j) = yh[i+nr(j-1)].           */
  75. /*      again, matrix(i,j) = vector(i+nr*(j-1))                          */
  76. /*      where nr is number of row of a matrix.                           */
  77. /*          yh[1] = yh(1,1) = rwork[21]                                  */
  78. /*          yh[2] = yh(2,1) = rwork[22]                                  */
  79. /*          yh[3] = yh(1,2) = rwork[23]                                  */
  80. /*          yh[4] = yh(2,2) = rwork[24]                                  */
  81. /*                                                                       */
  82. /*      unfortunately, it is not hidden to the end user of clsode.       */
  83. /*      unlike pd(i,j) in fortran jac (routine for jacobian matrix),     */
  84. /*      in clsode, pd is a vector, double pd[].  use following macro     */
  85. /*      to avoid the confusion.                                          */
  86. /*                                                                       */
  87. /*      array arithmatic macro                                           */
  88. /*            #define  pd(a,b)  pd[a+nr*(b-1)]                           */
  89. /*                                                                       */
  90. /*            in fortran               in c                              */
  91. /*              pd(1,1)               pd(1,1) --> pd[1]                  */
  92. /*              pd(2,1)               pd(2,1) --> pd[2]                  */
  93. /*              pd(1,2)               pd(1,2) --> pd[1+nr]               */
  94. /*                                                                       */
  95. /*  certain compiler may issue a warning for pd(i,1), such as            */
  96. /*  "code has no effect", but don't worry about it.  it is just from     */
  97. /*  the macro, pd(a,nr*(1-1)).                                           */
  98. /*                                                                       */
  99. /*  clsode has one more argument than flsode. the last argument, which   */
  100. /*  is "stdout" in example program, is a predefined stream in c.         */
  101. /*  stdout is generally ok.  but if you want to print the error          */
  102. /*  messages elsewhere, you could open a file, and connect to lsode.     */
  103. /*  if you open it, don't forget to close it.                            */
  104. /*                                                                       */
  105. /*  possible choices                                                     */
  106. /*      stdout     a standard output device.                             */
  107. /*      stderr     a standard error output device.                       */
  108. /*      FILE *     a file handle.                                        */
  109. /*                                                                       */
  110. /*  in clsode, two of lsode()'s argument carry the address operators,    */
  111. /*  &, to receive the the values, independent variable and istate.       */
  112. /*                                                                       */
  113. /*  fortran: lsode(fcn,neq,y, t,tout,itol,rtol,atol, itask, istate,      */
  114. /*  c(*):    lsode(fcn,neq,y,&T,tout,itol,rtol,ATOLE,itask,&ISTATE,      */
  115. /*  fortran:           iopt,rwork,lrw,iwork,liw,jac,mf)                  */
  116. /*  c(*):              iopt,rwork,lrw,iwork,liw,jac,mf,STDOUT);          */
  117. /*                                                                       */
  118. /*       *. uppercases are used just for the purpose of emphasizing.     */
  119. /*                                                                       */
  120. /* --------------------------------------------------------------------- */
  121. /*  rtol   = a relative error tolerance parameter, in fortran either a   */
  122. /*           scalar or an array of length neq.  in c, it is always a     */
  123. /*           vector array, but size is depended on itol, either neq[1]   */
  124. /*           or 1.  (actually 2 including a zeroth component(*).)        */
  125. /*                                                                       */
  126. /*  atole  = an absolute error tolerance parameter, in fortran either    */
  127. /*           a scalar or an array of length neq.  in c, it is always a   */
  128. /*           vector array, but size is depended on itol, either 1 or     */
  129. /*           neq[1].                                                     */
  130. /*                                                                       */
  131. /*           certain compiler generates an error with atol[], because    */
  132. /*           atol() is a c function defined in stdlib.h.  hence, i       */
  133. /*           replaced atol[] by atole[].  this will be applied to the    */
  134. /*           routines, lsode, ewset, and your main program.              */
  135. /*                                                                       */
  136. /*  neq    = the size of the ode system (number of first order           */
  137. /*           ordinary differential equations).  used only for input.     */
  138. /*                                                                       */
  139. /*           always, neq is a vector array, with neq[1] set to the       */
  140. /*           system size. (the clsode package accesses only neq[1].)     */
  141. /*           this parameter is passed as the neq argument in all calls   */
  142. /*           to f and jac.  hence, locations neq[2],... may be used to   */
  143. /*           store other integer data and pass it to f and/or jac.       */
  144. /*           function fcn and/or jac must include neq[] instead of neq   */
  145. /*           in a function argument. minimum size of neq array is 1.     */
  146. /*           (actually 2 including a zeroth component(*).)               */
  147. /*                                                                       */
  148. /*           *. all elements of the array are referenced from 0          */
  149. /*              instead of 1. (i.e. neq[0], neq[1], ...., neq[n]).       */
  150. /*                                                                       */
  151. /*  ***. do not optimize the common subexpression.                       */
  152. /* --------------------------------------------------------------------- */
  153. /*                                                                       */
  154. /*  example program.                                                     */
  155. /*    - it is the same one that f-lsode has in the comment statements.   */
  156. /*                                                                       */
  157. /*  the following is a simple example problem, with the coding           */
  158. /*  needed for its solution by lsode.  the problem is from chemical      */
  159. /*  kinetics, and consists of the following three rate equations..       */
  160. /*      dy1/dt = -.04*y1 + 1.e4*y2*y3                                    */
  161. /*      dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2                        */
  162. /*      dy3/dt = 3.e7*y2**2                                              */
  163. /*  on the interval from t = 0.0 to t = 4.e10, with initial conditions   */
  164. /*  y1 = 1.0, y2 = y3 = 0.  the problem is stiff.                        */
  165. /*                                                                       */
  166. /*  the following coding solves this problem with lsode, using mf = 21   */
  167. /*  and printing results at t = .4, 4., ..., 4.e10.  it uses             */
  168. /*  itol = 2 and atol much smaller for y2 than y1 or y3 because          */
  169. /*  y2 has much smaller values.                                          */
  170. /*  at the end of the run, statistical quantities of interest are        */
  171. /*  printed (see optional outputs in the full description).              */
  172. /*                                                                       */
  173. /* --------------------------------------------------------------------- */
  174. /*  #include <stdlib.h>                                                  */
  175. /*  #include <stdio.h>                                                   */
  176. /*  #include "lsode.h"                                                   */
  177. /*                                                                       */
  178. /*                                                                       */
  179. /*  void fcn(int[], double, double[], double[]);                         */
  180. /*  void jac(int[], double, double[], int, int, double[], int);          */
  181. /*                                                                       */
  182. /*                                                                       */
  183. /*  int main()                                                           */
  184. /*  {                                                                    */
  185. /*      double atole[3+1],rtol[2],*rwork,t,tout,*y ;                     */
  186. /*      int itol,itask,istate,iopt,lrw,liw,mf,*iwork,neq[2];             */
  187. /*      int i;                                                           */
  188. /*                                                                       */
  189. /*        itol = 2;                                                      */
  190. /*        lrw = 58;                                                      */
  191. /*        liw = 23;                                                      */
  192. /*        neq[1] = 3;                                                    */
  193. /*                                                                       */
  194. /*        iwork = ivector(liw);                                          */
  195. /*        rwork = dvector(lrw);                                          */
  196. /*        y     = dvector(neq[1]);                                       */
  197. /*                                                                       */
  198. /*        y[1] = 1.0e0;                                                  */
  199. /*        y[2] = 0.0e0;                                                  */
  200. /*        y[3] = 0.0e0;                                                  */
  201. /*        t    = 0.0e0;                                                  */
  202. /*        tout = 0.40e0;                                                 */
  203. /*                                                                       */
  204. /*        rtol[1]  = 1.e-4;                                              */
  205. /*        atole[1] = 1.e-6;                                              */
  206. /*        atole[2] = 1.e-10;                                             */
  207. /*        atole[3] = 1.e-6;                                              */
  208. /*                                                                       */
  209. /*        itask  = 1;                                                    */
  210. /*        istate = 1;                                                    */
  211. /*        iopt   = 0;                                                    */
  212. /*        mf     = 21;                                                   */
  213. /*                                                                       */
  214. /*        for(i=1;i<=12;i++) {                                           */
  215. /*            lsode(fcn,neq,y,&t,tout,itol,rtol,atole,itask,&istate,     */
  216. /*                    iopt,rwork,lrw,iwork,liw,jac,mf,stdout);           */
  217. /*            printf("\n at t =%11.4e  y =%14.6e  %14.6e %14.6e",        */
  218. /*                     t, y[1], y[2], y[3]);                             */
  219. /*            if (istate < 0) {                                          */
  220. /*                printf("\n\n error halt.. istate =%3d",istate);        */
  221. /*                break;                                                 */
  222. /*            }                                                          */
  223. /*            tout *= 10.0e0;                                            */
  224. /*       }                                                               */
  225. /*                                                                       */
  226. /*       printf("\n\n no. steps = %4d  no. f-s = %4d  no. j-s =%4d\n",   */
  227. /*                     iwork[11], iwork[12], iwork[13]);                 */
  228. /*       free(iwork);                                                    */
  229. /*       free(rwork);                                                    */
  230. /*       free(y);                                                        */
  231. /*       return 1;                                                       */
  232. /*  }                                                                    */
  233. /*                                                                       */
  234. /*                                                                       */
  235. /*  void fcn(int neq[], double t, double y[], double ydot[])             */
  236. /*  {                                                                    */
  237. /*       ydot[1] = -.04e0*y[1] + 1.e4*y[2]*y[3];                         */
  238. /*       ydot[3] = 3.e7*y[2]*y[2];                                       */
  239. /*       ydot[2] = -ydot[1] - ydot[3];                                   */
  240. /*       return;                                                         */
  241. /*  }                                                                    */
  242. /*                                                                       */
  243. /*                                                                       */
  244. /*  #define   pd(a,b)  pd[nr*(b-1)+a]                                    */
  245. /*  void jac(int neq[], double t, double y[], int ml,                    */
  246. /*                  int mu, double pd[], int nr)                         */
  247. /*  {                                                                    */
  248. /*       pd(1,1) = -0.04e0;                                              */
  249. /*       pd(1,2) = 1.e4*y[3];                                            */
  250. /*       pd(1,3) = 1.e4*y[2];                                            */
  251. /*       pd(2,1) = 0.04e0;                                               */
  252. /*       pd(2,3) = -pd(1,3);                                             */
  253. /*       pd(3,2) = 6.e7*y[2];                                            */
  254. /*       pd(2,2) = -pd(1,2) - pd(3,2);                                   */
  255. /*       return;                                                         */
  256. /*  }                                                                    */
  257. /*                                                                       */
  258. /* --------------------------------------------------------------------- */
  259. /*  the output(*) of this program is as follows..                        */
  260. /*                                                                       */
  261. /*  at t = 4.0000e-01  y =  9.851726e-01   3.386406e-05   1.479357e-02   */
  262. /*  at t = 4.0000e+00  y =  9.055142e-01   2.240418e-05   9.446344e-02   */
  263. /*  at t = 4.0000e+01  y =  7.158050e-01   9.184616e-06   2.841858e-01   */
  264. /*  at t = 4.0000e+02  y =  4.504846e-01   3.222434e-06   5.495122e-01   */
  265. /*  at t = 4.0000e+03  y =  1.831701e-01   8.940379e-07   8.168290e-01   */
  266. /*  at t = 4.0000e+04  y =  3.897016e-02   1.621193e-07   9.610297e-01   */
  267. /*  at t = 4.0000e+05  y =  4.935213e-03   1.983756e-08   9.950648e-01   */
  268. /*  at t = 4.0000e+06  y =  5.159269e-04   2.064759e-09   9.994841e-01   */
  269. /*  at t = 4.0000e+07  y =  5.306413e-05   2.122677e-10   9.999469e-01   */
  270. /*  at t = 4.0000e+08  y =  5.494530e-06   2.197824e-11   9.999945e-01   */
  271. /*  at t = 4.0000e+09  y =  5.129458e-07   2.051784e-12   9.999995e-01   */
  272. /*  at t = 4.0000e+10  y = -7.170561e-08  -2.868224e-13   1.000000e+00   */
  273. /*                                                                       */
  274. /*  no. steps =  330  no. f-s =  405  no. j-s =  69                      */
  275. /*                                                                       */
  276. /*  *. the output is the same as that from f-lsode. (see teste.doc)      */
  277. /*     however, it may vary slightly from compiler to compiler.          */
  278. /*                                                                       */
  279. /* --------------------------------------------------------------------- */
  280.  
  281. #include <stdio.h>
  282.  
  283. #define  NOREF(a)    (a=a)
  284. #define  FALSE 0
  285. #define  TRUE  1
  286. #define  max(a,b)    (((a) > (b)) ? (a) : (b))
  287. #define  min(a,b)    (((a) < (b)) ? (a) : (b))
  288.  
  289. extern void Xermsg(char*, char*, FILE*, int, int, char*, ...);
  290. extern double d1mach(void);
  291. extern 
  292. int lsode (void (*)(int [], double, double [], double []),
  293.            int [], double [], double*, double, int, double [],
  294.            double [], int, int*, int, double [], int, int [], int,
  295.            void (*)(int [],double,double [],int,int,double [],int),
  296.            int, FILE *);
  297.  
  298. extern double dsign(double, double);
  299. extern void tol_array_size(int, int, int*, int*);
  300. extern double* dvector(int);
  301. extern int* ivector(int);
  302.  
  303. /*-------------------------- end of lsode.h ---------------------------- */
  304.  
  305.