1c-----------------------------------------------------------------------
2c Demonstration program for the DLSODES package.
3c This is the version of 14 June 2001.
4c
5c This version is in double precision.
6c
7c The package is used for each of the relevant values of mf to solve
8c the problem ydot = A * y, where A is the 9 by 9 sparse matrix
9c
10c               -4  1     1
11c                1 -4  1     1
12c                   1 -4        1
13c                        -4  1     1
14c       A =               1 -4  1     1
15c                            1 -4        1
16c                                 -4  1
17c                                  1 -4  1
18c                                     1 -4
19c
20c The initial conditions are  y(0) = (1, 2, 3, ..., 9).
21c Output is printed at t = 1, 2, and 3.
22c Each case is solved first with nominal (large) values of lrw and liw,
23c and then with values given by lenrw and leniw (optional outputs)
24c on the first run, as a check on these computed work array lengths.
25c If the errors are too large, or other difficulty occurs,
26c a warning message is printed.
27c All output is on unit lout, which is data-loaded to 6 below.
28c-----------------------------------------------------------------------
29      external fdem, jdem
30      integer i, ia, igrid, iopt, iout, irun, istate, itask, itol,
31     1  iwork, j, ja, k, l, leniw, lenrw, liw, lout, lrw,
32     2  m, meth, mf, miter, moss, neq, nerr, nfe, nfea,
33     3  ngp, nje, nlu, nnz, nout, nqu, nst, nzl, nzu
34      double precision atol, erm, ero, hu, rtol, rwork, t, tout, y
35      dimension y(9), ia(10), ja(50), iwork(90), rwork(1000)
36      equivalence (ia(1),iwork(31)), (ja(1),iwork(41))
37      data lout/6/
38c
39c Write heading and set fixed parameters.
40      write(lout,10)
41 10   format(/'Demonstration problem for the DLSODES package'//)
42      nerr = 0
43      igrid = 3
44      neq = igrid**2
45      t = 0.0d0
46      itol = 1
47      rtol = 0.0d0
48      atol = 1.0d-5
49      itask = 1
50      iopt = 0
51      do 20 i = 1,neq
52 20     y(i) = i
53      ia(1) = 1
54      k = 1
55      do 60 m = 1,igrid
56        do 50 l = 1,igrid
57          j = l + (m - 1)*igrid
58          if (m .gt. 1) then
59            ja(k) = j - igrid
60            k = k + 1
61          endif
62 30       if (l .gt. 1) then
63            ja(k) = j - 1
64            k = k + 1
65          endif
66 35       ja(k) = j
67          k = k + 1
68          if (l .lt. igrid) then
69            ja(k) = j + 1
70            k = k + 1
71          endif
72 40       ia(j+1) = k
73 50       continue
74 60     continue
75      write (lout,80)neq,t,rtol,atol,(y(i),i=1,neq)
76 80   format(' neq =',i4,5x,'t0 =',f4.1,5x,'rtol =',d12.3,5x,
77     1   'atol =',d12.3//' Initial y vector =  ',9f5.1)
78c
79c Loop over all relevant values of mf.
80      do 193 moss = 0,2
81      do 192 meth = 1,2
82      do 191 miter = 0,3
83      if ( (miter.eq.0 .or. miter.eq.3) .and. moss.ne.0) go to 191
84      mf = 100*moss + 10*meth + miter
85      write (lout,100)
86 100  format(//80('*'))
87c First run: nominal work array lengths, 3 output points.
88      irun = 1
89      lrw = 1000
90      liw = 90
91      nout = 3
92 110  continue
93      write (lout,120)mf,lrw,liw
94 120  format(//'Run with mf =',i4,'.',5x,
95     1       'Input work lengths lrw, liw =',2i6/)
96      do 125 i = 1,neq
97 125    y(i) = i
98      t = 0.0d0
99      tout = 1.0d0
100      istate = 1
101      ero = 0.0d0
102c Loop over output points.  Do output and accuracy check at each.
103      do 170 iout = 1,nout
104        call dlsodes (fdem, neq, y, t, tout, itol, rtol, atol, itask,
105     1                istate, iopt, rwork, lrw, iwork, liw, jdem, mf)
106        nst = iwork(11)
107        hu = rwork(11)
108        nqu = iwork(14)
109        call edit (y, iout, erm)
110        write(lout,140)t,nst,hu,nqu,erm,(y(i),i=1,neq)
111 140    format('At t =',f5.1,3x,'nst =',i4,3x,'hu =',d12.3,3x,
112     1    'nqu =',i3,3x,' max. err. =',d11.3/
113     2    '  y array =    ',4d15.6/5d15.6)
114        if (istate .lt. 0) go to 175
115        erm = erm/atol
116        ero = max(ero,erm)
117        if (erm .gt. 100.0d0) then
118          write (lout,150)
119 150      format(//' Warning: error exceeds 100 * tolerance'//)
120          nerr = nerr + 1
121        endif
122        tout = tout + 1.0d0
123 170    continue
124 175  continue
125      if (istate .lt. 0) nerr = nerr + 1
126      if (irun .eq. 2) go to 191
127c Print final statistics (first run only)
128      nst = iwork(11)
129      nfe = iwork(12)
130      nje = iwork(13)
131      lenrw = iwork(17)
132      leniw = iwork(18)
133      nnz = iwork(19)
134      ngp = iwork(20)
135      nlu = iwork(21)
136      nzl = iwork(25)
137      nzu = iwork(26)
138      nfea = nfe
139      if (miter .eq. 2) nfea = nfe - ngp*nje
140      if (miter .eq. 3) nfea = nfe - nje
141      write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
142 180  format(/'Final statistics for this run:'/
143     1  ' rwork size =',i4,'   iwork size =',i4/
144     2  ' number of steps =',i5/
145     3  ' number of f-s   =',i5/
146     4  ' (excluding J-s) =',i5/
147     5  ' number of J-s   =',i5/
148     6  ' error overrun =',d10.2)
149      if (miter .eq. 1 .or. miter .eq. 2)
150     1   write (lout,185)nnz,ngp,nlu,nzl,nzu
151 185  format(' number of nonzeros in J = ',i5/
152     1  ' number of J index groups =',i5/
153     2  ' number of LU decomp-s    =',i5/
154     3  ' nonzeros in strict lower factor =',i5/
155     4  ' nonzeros in strict upper factor =',i5)
156      if (istate .lt. 0) go to 191
157      if (miter .eq. 1 .or. miter .eq. 2)
158     1   call ssout (neq, rwork(21), iwork, lout)
159c Return for second run: minimal work array lengths, 1 output point.
160      irun = irun + 1
161      lrw = lenrw
162      liw = leniw
163      nout = 1
164      go to 110
165 191  continue
166 192  continue
167 193  continue
168c
169      write (lout,100)
170      write (lout,200) nerr
171 200  format(//'Number of errors encountered =',i3)
172      stop
173      end
174
175      subroutine fdem (neq, t, y, ydot)
176      integer neq,  i, igrid, j, l, m
177      double precision t, y, ydot
178      dimension y(neq), ydot(neq)
179      data igrid/3/
180      do 5 i = 1,neq
181 5      ydot(i) = 0.0d0
182      do 20 m = 1,igrid
183        do 10 l = 1,igrid
184          j = l + (m - 1)*igrid
185          if (m .ne. 1) ydot(j-igrid) = ydot(j-igrid) + y(j)
186          if (l .ne. 1) ydot(j-1) = ydot(j-1) + y(j)
187          ydot(j) = ydot(j) - 4.0d0*y(j)
188          if (l .ne. igrid) ydot(j+1) = ydot(j+1) + y(j)
189 10       continue
190 20     continue
191      return
192      end
193
194      subroutine jdem (neq, t, y, j, ia, ja, pdj)
195      integer neq, j, ia, ja,  igrid, l, m
196      double precision t, y, pdj
197      dimension y(neq), ia(*), ja(*), pdj(neq)
198      data igrid/3/
199      m = (j - 1)/igrid + 1
200      l = j - (m - 1)*igrid
201      pdj(j) = -4.0d0
202      if (m .ne. 1) pdj(j-igrid) = 1.0d0
203      if (l .ne. 1) pdj(j-1) = 1.0d0
204      if (l .ne. igrid) pdj(j+1) = 1.0d0
205      return
206      end
207
208      subroutine edit (y, iout, erm)
209      integer iout,  i, neq
210      double precision y, erm,   er, yex
211      dimension y(*),yex(9,3)
212      data neq /9/
213      data yex /6.687279d-01, 9.901910d-01, 7.603061d-01,
214     1   8.077979d-01, 1.170226e+00, 8.810605d-01, 5.013331d-01,
215     2   7.201389d-01, 5.379644d-01, 1.340488d-01, 1.917157d-01,
216     3   1.374034d-01, 1.007882d-01, 1.437868d-01, 1.028010d-01,
217     4   3.844343d-02, 5.477593d-02, 3.911435d-02, 1.929166d-02,
218     5   2.735444d-02, 1.939611d-02, 1.055981d-02, 1.496753d-02,
219     6   1.060897d-02, 2.913689d-03, 4.128975d-03, 2.925977d-03/
220      erm = 0.0d0
221      do 10 i = 1,neq
222        er = abs(y(i) - yex(i,iout))
223 10     erm = max(erm,er)
224      return
225      end
226
227      subroutine ssout (neq, iwk, iwork, lout)
228      integer neq, iwk, iwork, lout
229      integer i, i1, i2, ipian, ipjan, nnz
230      dimension iwk(*), iwork(*)
231      ipian = iwork(23)
232      ipjan = iwork(24)
233      nnz = iwork(19)
234      i1 = ipian
235      i2 = i1 + neq
236      write (lout,10)(iwk(i),i=i1,i2)
237 10   format(/' structure descriptor array ian ='/(20i4))
238      i1 = ipjan
239      i2 = i1 + nnz - 1
240      write (lout,20)(iwk(i),i=i1,i2)
241 20   format(/' structure descriptor array jan ='/(20i4))
242      return
243      end
244