1c-----------------------------------------------------------------------
2c Demonstration program for the DLSODA 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 the 2 appropriate values of jt 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, f2, jac2
14      integer i, iopar, iopt, iout, istate, itask, itol, iwork,
15     1   jt, leniw, lenrw, liw, lout, lrw, mband, mused,
16     2   ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst
17      double precision atol, dtout, dtout0, dtout1, er, erm, ero, hu,
18     1     rtol, rwork, t, tout, tout1, tsw, y
19      dimension y(25), rwork(522), iwork(45)
20      data lout/6/, tout1/16.921743d0/, dtout/17.341162d0/
21c
22      nerr = 0
23      itol = 1
24      rtol = 0.0d0
25      atol = 1.0d-8
26      lrw = 522
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 DLSODA package'////
36     1  ' Problem 1:   Van der Pol oscillator:'/
37     2  '              xdotdot - 20*(1 - x**2)*xdot + x = 0, ',
38     3  '   x(0) = 2, xdot(0) = 0'/' neq =',i2/
39     4  ' itol =',i3,'   rtol =',d10.1,'   atol =',d10.1//)
40c
41      do 190 jt = 1,2
42      write (lout,120) jt
43 120  format(//' Solution with jt =',i3//
44     1       '  t               x               xdot       meth',
45     2       '   nq     h           tsw'//)
46      t = 0.0d0
47      y(1) = 2.0d0
48      y(2) = 0.0d0
49      itask = 1
50      istate = 1
51      dtout0 = 0.5d0*tout1
52      dtout1 = 0.5d0*dtout
53      tout = dtout0
54      ero = 0.0d0
55      do 170 iout = 1,nout
56        call dlsoda(f1,neq,y,t,tout,itol,rtol,atol,itask,istate,
57     1              iopt,rwork,lrw,iwork,liw,jac1,jt)
58        hu = rwork(11)
59        tsw = rwork(15)
60        nqu = iwork(14)
61        mused = iwork(19)
62        write (lout,140) t,y(1),y(2),mused,nqu,hu,tsw
63 140    format(d12.5,d16.5,d14.3,2i6,2d13.3)
64        if (istate .lt. 0) go to 175
65        iopar = iout - 2*(iout/2)
66        if (iopar .ne. 0) go to 160
67        er = abs(y(1))
68        ero = max(ero,er)
69        if (er .gt. 1.0d-2) then
70          write (lout,150)
71 150      format(//' Warning: value at root exceeds 1.0d-2'//)
72          nerr = nerr + 1
73        endif
74 160    if (iout .eq. 1) tout = tout + dtout0
75        if (iout .gt. 1) tout = tout + dtout1
76 170    continue
77 175  continue
78      if (istate .lt. 0) nerr = nerr + 1
79      nst = iwork(11)
80      nfe = iwork(12)
81      nje = iwork(13)
82      lenrw = iwork(17)
83      leniw = iwork(18)
84      nfea = nfe
85      if (jt .eq. 2) nfea = nfe - neq*nje
86      write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
87 180  format(//' Final statistics for this run:'/
88     1  ' rwork size =',i4,'   iwork size =',i4/
89     2  ' number of steps =',i5/
90     3  ' number of f-s   =',i5/
91     4  ' (excluding J-s) =',i5/
92     5  ' number of J-s   =',i5/
93     6  ' max. error at root =',d10.2)
94 190  continue
95c
96c Second problem
97c
98      neq = 25
99      ml = 5
100      mu = 0
101      iwork(1) = ml
102      iwork(2) = mu
103      mband = ml + mu + 1
104      atol = 1.0d-6
105      nout = 5
106      write (lout,210) neq,ml,mu,itol,rtol,atol
107 210  format(///80('-')///
108     1  ' Problem 2: ydot = A * y , where',
109     2  '  A is a banded lower triangular matrix'/
110     2  '            derived from 2-D advection PDE'/
111     3  ' neq =',i3,'   ml =',i2,'   mu =',i2/
112     4  ' itol =',i3,'   rtol =',d10.1,'   atol =',d10.1//)
113      do 290 jt = 4,5
114      write (lout,220) jt
115 220  format(//' Solution with jt =',i3//
116     1       '     t             max.err.     meth   ',
117     2       'nq      h            tsw'//)
118      t = 0.0d0
119      do 230 i = 2,neq
120 230    y(i) = 0.0d0
121      y(1) = 1.0d0
122      itask = 1
123      istate = 1
124      tout = 0.01d0
125      ero = 0.0d0
126      do 270 iout = 1,nout
127        call dlsoda(f2,neq,y,t,tout,itol,rtol,atol,itask,istate,
128     1              iopt,rwork,lrw,iwork,liw,jac2,jt)
129        call edit2(y,t,erm)
130        hu = rwork(11)
131        tsw = rwork(15)
132        nqu = iwork(14)
133        mused = iwork(19)
134        write (lout,240) t,erm,mused,nqu,hu,tsw
135 240    format(d15.5,d14.3,2i6,2d14.3)
136        if (istate .lt. 0) go to 275
137        er = erm/atol
138        ero = max(ero,er)
139        if (er .gt. 1000.0d0) then
140          write (lout,150)
141          nerr = nerr + 1
142        endif
143 270    tout = tout*10.0d0
144 275  continue
145      if (istate .lt. 0) nerr = nerr + 1
146      nst = iwork(11)
147      nfe = iwork(12)
148      nje = iwork(13)
149      lenrw = iwork(17)
150      leniw = iwork(18)
151      nfea = nfe
152      if (jt .eq. 5) nfea = nfe - mband*nje
153      write (lout,280) lenrw,leniw,nst,nfe,nfea,nje,ero
154 280  format(//' Final statistics for this run:'/
155     1  ' rwork size =',i4,'   iwork size =',i4/
156     2  ' number of steps =',i5/
157     3  ' number of f-s   =',i5/
158     4  ' (excluding J-s) =',i5/
159     5  ' number of J-s   =',i5/
160     6  ' error overrun =',d10.2)
161 290  continue
162      write (lout,300) nerr
163 300  format(///' Number of errors encountered =',i3)
164      stop
165      end
166
167      subroutine f1 (neq, t, y, ydot)
168      integer neq
169      double precision t, y, ydot
170      dimension y(neq), ydot(neq)
171      ydot(1) = y(2)
172      ydot(2) = 20.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1)
173      return
174      end
175
176      subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd)
177      integer neq, ml, mu, nrowpd
178      double precision t, y, pd
179      dimension y(neq), pd(nrowpd,neq)
180      pd(1,1) = 0.0d0
181      pd(1,2) = 1.0d0
182      pd(2,1) = -40.0d0*y(1)*y(2) - 1.0d0
183      pd(2,2) = 20.0d0*(1.0d0 - y(1)*y(1))
184      return
185      end
186
187      subroutine f2 (neq, t, y, ydot)
188      integer neq, i, j, k, ng
189      double precision t, y, ydot, alph1, alph2, d
190      dimension y(neq), ydot(neq)
191      data alph1/1.0d0/, alph2/1.0d0/, ng/5/
192      do 10 j = 1,ng
193      do 10 i = 1,ng
194        k = i + (j - 1)*ng
195        d = -2.0d0*y(k)
196        if (i .ne. 1) d = d + y(k-1)*alph1
197        if (j .ne. 1) d = d + y(k-ng)*alph2
198 10     ydot(k) = d
199      return
200      end
201
202      subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
203      integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng
204      double precision t, y, pd, alph1, alph2
205      dimension y(neq), pd(nrowpd,neq)
206      data alph1/1.0d0/, alph2/1.0d0/, ng/5/
207      mband = ml + mu + 1
208      mu1 = mu + 1
209      mu2 = mu + 2
210      do 10 j = 1,neq
211        pd(mu1,j) = -2.0d0
212        pd(mu2,j) = alph1
213 10     pd(mband,j) = alph2
214      do 20 j = ng,neq,ng
215 20     pd(mu2,j) = 0.0d0
216      return
217      end
218
219      subroutine edit2 (y, t, erm)
220      integer i, j, k, ng
221      double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt
222      dimension y(*)
223      data alph1/1.0d0/, alph2/1.0d0/, ng/5/
224      erm = 0.0d0
225      if (t .eq. 0.0d0) return
226      ex = 0.0d0
227      if (t .le. 30.0d0) ex = exp(-2.0d0*t)
228      a2 = 1.0d0
229      do 60 j = 1,ng
230        a1 = 1.0d0
231        do 50 i = 1,ng
232          k = i + (j - 1)*ng
233          yt = t**(i+j-2)*ex*a1*a2
234          er = abs(y(k)-yt)
235          erm = max(erm,er)
236          a1 = a1*alph1/i
237 50       continue
238        a2 = a2*alph2/j
239 60     continue
240      return
241      end
242