1c-----------------------------------------------------------------------
2c Demonstration program for the DLSODE package.
3c This is the version of 14 June 2001.
4c
5c This version is in double precision.
6c
7c The package is used to solve two simple problems,
8c one with a full Jacobian, the other with a banded Jacobian,
9c with all 8 of the appropriate values of mf in each case.
10c If the errors are too large, or other difficulty occurs,
11c a warning message is printed.  All output is on unit lout = 6.
12c-----------------------------------------------------------------------
13      external f1, jac1
14      integer i, iopar, iopt, iout, istate, itask, itol, iwork,
15     1   leniw, lenrw, liw, lout, lrw, mband, meth, mf, miter,
16     2   ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst
17      double precision atol, dtout, er, erm, ero, hu, rtol, rwork, t,
18     1   tout, tout1, y
19      dimension y(25), rwork(697), iwork(45)
20      data lout/6/, tout1/1.39283880203d0/, dtout/2.214773875d0/
21c
22      nerr = 0
23      itol = 1
24      rtol = 0.0d0
25      atol = 1.0d-6
26      lrw = 697
27      liw = 45
28      iopt = 0
29c
30c First problem
31c
32      neq = 2
33      nout = 4
34      write (lout,110) neq,itol,rtol,atol
35 110  format(/' Demonstration program for DLSODE package'///
36     1  ' Problem 1:  Van der Pol oscillator:'/
37     2  '  xdotdot - 3*(1 - x**2)*xdot + x = 0, ',
38     3  '   x(0) = 2, xdot(0) = 0'/
39     4  ' neq =',i2/
40     5  ' itol =',i3,'   rtol =',d10.1,'   atol =',d10.1//)
41c
42      do 195 meth = 1,2
43      do 190 miter = 0,3
44      mf = 10*meth + miter
45      write (lout,120) mf
46 120  format(///' Solution with mf =',i3//
47     1     5x,'t               x               xdot       nq      h'//)
48      t = 0.0d0
49      y(1) = 2.0d0
50      y(2) = 0.0d0
51      itask = 1
52      istate = 1
53      tout = tout1
54      ero = 0.0d0
55      do 170 iout = 1,nout
56        call dlsode(f1,neq,y,t,tout,itol,rtol,atol,itask,istate,
57     1     iopt,rwork,lrw,iwork,liw,jac1,mf)
58        hu = rwork(11)
59        nqu = iwork(14)
60        write (lout,140) t,y(1),y(2),nqu,hu
61 140    format(d15.5,d16.5,d14.3,i5,d14.3)
62        if (istate .lt. 0) go to 175
63        iopar = iout - 2*(iout/2)
64        if (iopar .ne. 0) go to 170
65        er = abs(y(1))/atol
66        ero = max(ero,er)
67        if (er .gt. 1000.0d0) then
68          write (lout,150)
69 150      format(//' Warning: error exceeds 1000 * tolerance'//)
70          nerr = nerr + 1
71        endif
72 170    tout = tout + dtout
73 175  continue
74      if (istate .lt. 0) nerr = nerr + 1
75      nst = iwork(11)
76      nfe = iwork(12)
77      nje = iwork(13)
78      lenrw = iwork(17)
79      leniw = iwork(18)
80      nfea = nfe
81      if (miter .eq. 2) nfea = nfe - neq*nje
82      if (miter .eq. 3) nfea = nfe - nje
83      write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
84 180  format(//' Final statistics for this run:'/
85     1  ' rwork size =',i4,'   iwork size =',i4/
86     2  ' number of steps =',i5/
87     3  ' number of f-s   =',i5/
88     4  ' (excluding J-s) =',i5/
89     5  ' number of J-s   =',i5/
90     6  ' error overrun =',d10.2)
91 190  continue
92 195  continue
93
94      write (lout,300) nerr
95 300  format(////' Number of errors encountered =',i3)
96      stop
97      end
98
99c$$$      subroutine f1 (neq, t, y, ydot)
100c$$$      integer neq
101c$$$      double precision t, y, ydot
102c$$$      dimension y(neq), ydot(neq)
103c$$$      ydot(1) = y(2)
104c$$$      ydot(2) = 3.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1)
105c$$$      return
106c$$$      end
107c$$$
108c$$$      subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd)
109c$$$      integer neq, ml, mu, nrowpd
110c$$$      double precision t, y, pd
111c$$$      dimension y(neq), pd(nrowpd,neq)
112c$$$      pd(1,1) = 0.0d0
113c$$$      pd(1,2) = 1.0d0
114c$$$      pd(2,1) = -6.0d0*y(1)*y(2) - 1.0d0
115c$$$      pd(2,2) = 3.0d0*(1.0d0 - y(1)*y(1))
116c$$$      return
117c$$$      end
118