1c-----------------------------------------------------------------------
2c Demonstration program for the DLSODIS package.
3c This is the version of 14 June 2001.
4c
5c This version is in double precision.
6c
7c This program solves a semi-discretized form of the Burgers equation,
8c
9c     u  = -(u*u/2)  + eta * u
10c      t           x          xx
11c
12c for  -1 .le. x .le. 1, t .ge. 0.
13c Here eta = 0.05.
14c Boundary conditions: u(-1,t) = u(1,t) and du/dx(-1,t) = du/dx(1,t).
15c Initial profile: square wave
16c     u(0,x) = 0    for 1/2 .lt. abs(x) .le. 1
17c     u(0,x) = 1/2  for abs(x) = 1/2
18c     u(0,x) = 1    for 0 .le. abs(x) .lt. 1/2
19c
20c An ODE system is generated by a simplified Galerkin treatment
21c of the spatial variable x.
22c
23c Reference:
24c R. C. Y. Chin, G. W. Hedstrom, and K. E. Karlsson,
25c A Simplified Galerkin Method for Hyperbolic Equations,
26c Math. Comp., vol. 33, no. 146 (April 1979), pp. 647-658.
27c
28c The problem is run with the DLSODIS package with a 12-node mesh,
29c for various appropriate values of the method flag mf.
30c Output is on unit lout, set to 6 in a data statement below.
31c
32c Problem specific data:
33c npts  = number of unknowns (npts = 0 mod 4)
34c nnz   = number of non-zeros in Jacobian before fill in
35c nnza  = number of non-zeros in Jacobian after fill in
36c lrwk  = length of real work array (taking into account fill in)
37c liwk  = length of integer work array
38c ipia  = pointer to ia in iw (ia(j) = iw(ipia+j-1)
39c ipja  = pointer to ja in iw (ja(j) = iw(ipja+j-1)
40c ipic  = pointer to ic in iw array (ic(j) = iw(ipic+j-1))
41c ipjc  = pointer to jc in iw array (jc(j) = iw(ipjc+j-1))
42c-----------------------------------------------------------------------
43      parameter (npts = 12, nnz = 3*npts, nnza = 5*npts)
44      parameter (lrwk = 20+3*nnza+28*npts)
45      parameter (liwk = 32+2*nnza+2*npts)
46      parameter (ipia = 31, ipja = 31+npts+1)
47      parameter (ll = 32+npts+nnz, ipic = ll,ipjc = ll+npts+1)
48c
49      external res, addasp, jacsp
50      integer i, io, istate, itol, iw, j, lout, lrw, liw, lyh,
51     1   meth, miter, mf, moss, n, nout, nerr, n14, n34
52      double precision eta, delta, r4d, eodsq, a, b, t, tout, tlast,
53     1    tinit, errfac
54      double precision atol, rtol, rw, y, ydoti, elkup
55      double precision zero, fourth, half, one, hun
56      dimension y(npts), ydoti(npts), tout(4), atol(2), rtol(2)
57      dimension rw(lrwk), iw(liwk)
58c Pass problem parameters in the Common block test1.
59      common /test1/ r4d, eodsq
60c
61c Set problem parameters and run parameters
62      data eta/0.05d0/, a/-1.0d0/, b/1.0d0/
63      data zero/0.0d0/, fourth/0.25d0/, half/0.5d0/, one/1.0d0/,
64     1   hun/100.0d0/
65      data tinit/0.0d0/, tlast/0.4d0/
66      data tout/0.10d0, 0.20d0, 0.30d0, 0.40d0/
67      data lout/6/, nout/4/
68      data itol/1/, rtol/1.0d-3, 1.0d-6/, atol/1.0d-3, 1.0d-6/
69c
70      nerr = 0
71      lrw = lrwk
72      liw = liwk
73c
74c Compute the mesh width delta and other parameters.
75      delta = (b - a)/npts
76      r4d = fourth/delta
77      eodsq = eta/delta**2
78      n14 = npts/4 + 1
79      n34 = 3 * (npts/4) + 1
80      n = npts
81c
82c Set the initial profile (for output purposes only).
83      do 10 i = 1,n
84   10   y(i) = zero
85      y(n14) = half
86      do 20 i = n14+1,n34-1
87   20   y(i) = one
88      y(n34) = half
89c
90      write(lout,1000)
91      write(lout,1100) eta,a,b,tinit,tlast,n
92      write(lout,1200) (y(i), i=1,n)
93c
94c Set the initial sparse data structures for coefficient matrix A
95c and the Jacobian matrix C
96      call struct(iw(ipia), iw(ipja), iw(ipic), iw(ipjc), n)
97c
98c The j loop is over error tolerances.
99      do 200 j = 1,2
100c
101c This method flag loop is for demonstration only.
102      do 200 moss = 0,4
103        do 100 meth = 1,2
104          do 100 miter = 1,2
105   35       mf = 100*moss + 10*meth + miter
106c
107c Set the initial profile.
108            do 40 i = 1,n
109   40         y(i) = zero
110            y(n14) = half
111            do 50 i = n14+1,n34-1
112   50         y(i) = one
113            y(n34) = half
114c
115            t = tinit
116            istate = 0
117c
118            write(lout,1500) itol, rtol(j), atol(j), mf
119c
120c output loop for each case
121            do 80 io = 1,nout
122c
123c Call DLSODIS and print results.
124              call dlsodis (res, addasp, jacsp, n, y, ydoti, t,
125     1             tout(io), itol, rtol(j), atol(j), 1, istate, 0,
126     2             rw, lrw, iw, liw, mf)
127              write(lout,2000) t, rw(11), iw(14), (y(i), i=1,n)
128c
129c If istate is not 2 on return, print message and go to next case.
130              if (istate .ne. 2) then
131                write(lout,4000) mf, t, istate
132                nerr = nerr + 1
133                go to 100
134                endif
135   80       continue
136            write(lout,3000) mf, iw(11), iw(12), iw(13),
137     1                       iw(17), iw(18), iw(20), iw(21)
138c
139c Estimate final error and print result.
140            lyh = iw(22)
141            errfac = elkup(n, y, rw(lyh), itol, rtol(j), atol(j), lout)
142            if (errfac .lt. hun) then
143              write(lout,5000) errfac
144            else
145              write(lout,5001) errfac
146              nerr = nerr + 1
147              endif
148  100     continue
149  200   continue
150c
151      write(lout,6000) nerr
152      stop
153c
154 1000 format(20x,' Demonstration Program for DLSODIS' )
155 1100 format(//10x,'-- Simplified Galerkin solution of ',
156     1       'Burgers equation --'///
157     1       13x,'Diffusion coefficient is eta =',d10.2/
158     1       13x,'Uniform mesh on interval',d12.3,' to ',d12.3/
159     2       13x,'Periodic boundary conditions'/
160     2       13x,'Initial data are as follows:'//20x,'t0 = ',d12.5/
161     2       20x,'tlast = ',d12.5/20x,'n  = ',i3//)
162c
163 1200 format(/'Initial profile:',/20(6d12.4/))
164c
165 1500 format(///85('*')///'Run with itol =',i2,'  rtol =',d12.2,
166     1       '  atol =',d12.2,'   mf = ',i3//)
167c
168 2000 format(' Output for time t =',d12.5,'  current h =',d12.5,
169     1       '  current order =',i2/20(6d12.4/))
170c
171 3000 format(/'Final statistics for mf = ',i3,': ',
172     1       i5,' steps,',i6,' res,',i6,' Jacobians,'/
173     2       20x,' rw size =',i6,',    iw size =',i6/
174     3       20x,i4,' extra res for each jac,',i4,' decomps')
175c
176 4000 format(/'Final time reached for mf = ',i3,
177     1       ' was t = ',d12.5/25x,'at which istate = ',i2//)
178 5000 format('Final output is correct to within ',d9.2,
179     1       '  times local error tolerance.'/)
180 5001 format('Final output is wrong by ',d9.2,
181     1       '  times local error tolerance.'/)
182 6000 format(///85('*')//
183     1       'Run completed: number of errors encountered =',i3)
184c
185c end of main program for the DLSODIS demonstration program
186      end
187
188      subroutine struct(ia, ja, ic, jc, n)
189c This subroutine computes the initial sparse data structure of
190c the mass (ia,ja) and Jacobian (ic,jc) matrices.
191c
192      integer ia(*), ja(*), ic(*), jc(*), n,  jj, k, l, m
193c
194      write(6,1200)
195      k = 0
196      do 33 l = 1,n
197         ia(l) = (l-1)*3+1
198         ic(l) = (l-1)*3+1
199         do 32 m = l,l+2
200            k = k + 1
201            ja(k) = m - 1
202            jc(k) = m - 1
203   32    continue
204   33 continue
205      ia(n+1) = 3*n + 1
206      ic(n+1) = 3*n+1
207      ja(1) = n
208      jc(1) = n
209      ja(k) = 1
210      jc(k) = 1
211c
212      write(6,1300) (ia(jj),jj=1,n+1)
213      write(6,1350) (ja(jj),jj=1,k)
214      write(6,1400) (ic(jj),jj=1,n+1)
215      write(6,1450) (jc(jj),jj=1,k)
216      return
2171200  format('Initial sparse data structures'/)
2181300  format(' ia  ',15i4/10(5x,15i4/))
2191350  format(' ja  ',15i4/10(5x,15i4/))
2201400  format(' ic  ',15i4/10(5x,15i4/))
2211450  format(' jc  ',15i4/10(5x,15i4/))
222      end
223
224      subroutine res (n, t, y, v, r, ires)
225c This subroutine computes the residual vector
226c   r = g(t,y) - A(t,y)*v .
227c If ires = -1, only g(t,y) is returned in r, since A(t,y) does
228c not depend on y.
229c No changes need to be made to this routine if n is changed.
230c
231      integer n, ires,  i
232      double precision t, y(n), v(n), r(n), r4d, eodsq, one, four, six,
233     1   fact1, fact4
234      common /test1/ r4d, eodsq
235      data one /1.0d0/, four /4.0d0/, six /6.0d0/
236c
237      call gfun (n, t, y, r)
238      if (ires .eq. -1) return
239c
240      fact1 = one/six
241      fact4 = four/six
242c
243      r(1) = r(1) - (fact4*v(1) + fact1*(v(2) + v(n)))
244      do 10 i = 2,n-1
245         r(i) = r(i) - (fact4*v(i) + fact1*(v(i-1) + v(i+1)))
246   10 continue
247      r(n) = r(n) - (fact4*v(n) + fact1*(v(1) + v(n-1)))
248      return
249c end of subroutine res for the DLSODIS demonstration program
250      end
251
252      subroutine gfun (n, t, y, g)
253c This subroutine computes the right-hand side function g(y,t).
254c It uses r4d = 1/(4*delta) and eodsq = eta/delta**2
255c from the Common block test1.
256c
257      integer n, i
258      double precision t, y(n), g(n), r4d, eodsq, two
259      common /test1/ r4d, eodsq
260      data two/2.0d0/
261c
262      g(1) = r4d*(y(n)**2 - y(2)**2) + eodsq*(y(2) - two*y(1) + y(n))
263c
264      do 20 i = 2,n-1
265        g(i) = r4d*(y(i-1)**2 - y(i+1)**2)
266     1        + eodsq*(y(i+1) - two*y(i) + y(i-1))
267   20   continue
268c
269      g(n) = r4d*(y(n-1)**2 - y(1)**2) + eodsq*(y(1)-two*y(n)+y(n-1))
270c
271      return
272c end of subroutine gfun for the DLSODIS demonstration program
273      end
274
275      subroutine addasp (n, t, y, j, ip, jp, pa)
276c This subroutine computes the sparse matrix A by columns, adds it to
277c pa, and returns the sum in pa.
278c The matrix A is periodic tridiagonal, of order n, with nonzero elements
279c (reading across) of  1/6, 4/6, 1/6, with 1/6 in the lower left and
280c upper right corners.
281c
282      integer n, j, ip(*), jp(*),  jm1, jp1
283      double precision t, y(n), pa(n), fact1, fact4, one, four, six
284      data one/1.0d0/, four/4.0d0/, six/6.0d0/
285c
286c Compute the elements of A.
287      fact1 = one/six
288      fact4 = four/six
289      jm1 = j - 1
290      jp1 = j + 1
291      if (j .eq. n) jp1 = 1
292      if (j .eq. 1) jm1 = n
293c
294c Add the matrix A to the matrix pa (sparse).
295      pa(j) = pa(j) + fact4
296      pa(jp1) = pa(jp1) + fact1
297      pa(jm1) = pa(jm1) + fact1
298      return
299c end of subroutine addasp for the DLSODIS demonstration program
300      end
301
302      subroutine jacsp (n, t, y, s, j, ip, jp, pdj)
303c This subroutine computes the Jacobian dg/dy = d(g-A*s)/dy by
304c columns in sparse matrix format.  Only nonzeros are loaded.
305c It uses r4d = 1/(4*delta) and eodsq = eta/delta**2 from the Common
306c block test1.
307c
308      integer n, j, ip(*), jp(*),  jm1, jp1
309      double precision t, y(n), s(n), pdj(n), r4d, eodsq, two, diag, r2d
310      common /test1/ r4d, eodsq
311      data two/2.0d0/
312c
313      diag = -two*eodsq
314      r2d = two*r4d
315      jm1 = j - 1
316      jp1 = j + 1
317      if (j .eq. 1) jm1 = n
318      if (j .eq. n) jp1 = 1
319c
320      pdj(jm1) = -r2d*y(j) + eodsq
321      pdj(j) = diag
322      pdj(jp1) = r2d*y(j) + eodsq
323      return
324c end of subroutine jacsp for the DLSODIS demonstration program
325      end
326
327      double precision function elkup (n, y, ewt, itol, rtol,atol, lout)
328c This routine looks up approximately correct values of y at t = 0.4,
329c ytrue = y12 or y120 depending on whether n = 12 or 120.
330c These were obtained by running with very tight tolerances.
331c The returned value is
332c   elkup = norm of [ (y - ytrue) / (rtol*abs(ytrue) + atol) ].
333c
334      integer n, itol, lout, i
335      double precision y(n), ewt(n), rtol, atol, y12(12), y120(120),
336     1    y120a(16), y120b(16), y120c(16), y120d(16), y120e(16),
337     2    y120f(16), y120g(16), y120h(8), dvnorm
338      equivalence (y120a(1),y120(1)), (y120b(1),y120(17)),
339     1      (y120c(1),y120(33)), (y120d(1),y120(49)),
340     1      (y120e(1),y120(65)),
341     1      (y120f(1),y120(81)), (y120g(1),y120(97)),
342     1      (y120h(1),y120(113))
343      data y12/
344     1 1.60581860d-02, 3.23063251d-02, 1.21903380d-01, 2.70943828d-01,
345     1 4.60951522d-01, 6.57571216d-01, 8.25154453d-01, 9.35644796d-01,
346     1 9.90167557d-01, 9.22421221d-01, 5.85764902d-01, 1.81112615d-01/
347      data y120a /
348     1 1.89009068d-02, 1.63261891d-02, 1.47080563d-02, 1.39263623d-02,
349     1 1.38901341d-02, 1.45336989d-02, 1.58129308d-02, 1.77017162d-02,
350     1 2.01886844d-02, 2.32742221d-02, 2.69677715d-02, 3.12854037d-02,
351     1 3.62476563d-02, 4.18776225d-02, 4.81992825d-02, 5.52360652d-02/
352      data y120b /
353     1 6.30096338d-02, 7.15388849d-02, 8.08391507d-02, 9.09215944d-02,
354     1 1.01792784d-01, 1.13454431d-01, 1.25903273d-01, 1.39131085d-01,
355     1 1.53124799d-01, 1.67866712d-01, 1.83334757d-01, 1.99502830d-01,
356     1 2.16341144d-01, 2.33816600d-01, 2.51893167d-01, 2.70532241d-01/
357      data y120c /
358     1 2.89693007d-01, 3.09332757d-01, 3.29407198d-01, 3.49870723d-01,
359     1 3.70676646d-01, 3.91777421d-01, 4.13124817d-01, 4.34670077d-01,
360     1 4.56364053d-01, 4.78157319d-01, 5.00000270d-01, 5.21843218d-01,
361     1 5.43636473d-01, 5.65330432d-01, 5.86875670d-01, 6.08223037d-01/
362      data y120d /
363     1 6.29323777d-01, 6.50129662d-01, 6.70593142d-01, 6.90667536d-01,
364     1 7.10307235d-01, 7.29467947d-01, 7.48106966d-01, 7.66183477d-01,
365     1 7.83658878d-01, 8.00497138d-01, 8.16665158d-01, 8.32133153d-01,
366     1 8.46875019d-01, 8.60868691d-01, 8.74096465d-01, 8.86545273d-01/
367      data y120e /
368     1 8.98206892d-01, 9.09078060d-01, 9.19160487d-01, 9.28460742d-01,
369     1 9.36989986d-01, 9.44763554d-01, 9.51800339d-01, 9.58122004d-01,
370     1 9.63751979d-01, 9.68714242d-01, 9.73031887d-01, 9.76725449d-01,
371     1 9.79811001d-01, 9.82297985d-01, 9.84186787d-01, 9.85466039d-01/
372      data y120f /
373     1 9.86109629d-01, 9.86073433d-01, 9.85291781d-01, 9.83673704d-01,
374     1 9.81099057d-01, 9.77414704d-01, 9.72431015d-01, 9.65919133d-01,
375     1 9.57609585d-01, 9.47193093d-01, 9.34324619d-01, 9.18631922d-01,
376     1 8.99729965d-01, 8.77242371d-01, 8.50830623d-01, 8.20230644d-01/
377      data y120g /
378     1 7.85294781d-01, 7.46035145d-01, 7.02662039d-01, 6.55609682d-01,
379     1 6.05541326d-01, 5.53327950d-01, 4.99999118d-01, 4.46670394d-01,
380     1 3.94457322d-01, 3.44389410d-01, 2.97337561d-01, 2.53964948d-01,
381     1 2.14705729d-01, 1.79770169d-01, 1.49170367d-01, 1.22758681d-01/
382      data y120h /
383     1 1.00271052d-01, 8.13689920d-02, 6.56761515d-02, 5.28075160d-02,
384     1 4.23908624d-02, 3.40811650d-02, 2.75691506d-02, 2.25853507d-02/
385c
386      if ((n-12)*(n-120) .ne. 0) go to 300
387      if (n .eq. 120) go to 100
388c
389c Compute local error tolerance using correct y (n = 12).
390      call dewset(n, itol, rtol, atol, y12, ewt)
391c
392c Invert ewt and replace y by the error, y - ytrue.
393      do 20  i = 1, 12
394        ewt(i) = 1.0d0/ewt(i)
395 20     y(i) = y(i) - y12(i)
396      go to 200
397c
398c Compute local error tolerance using correct y (n = 120).
399 100  call dewset( n, itol, rtol, atol, y120, ewt )
400c
401c Invert ewt and replace y by the error, y - ytrue.
402      do 120  i = 1, 120
403        ewt(i) = 1.0d0/ewt(i)
404 120    y(i) = y(i) - y120(i)
405c
406c Find weighted norm of the error.
407 200  elkup = dvnorm (n, y, ewt)
408      return
409c
410c error return
411 300  write(lout,400) n
412      elkup = 1.0d3
413 400  format(/5x,'Illegal use of elkup for n =',i4)
414      return
415c end of function elkup for the DLSODIS demonstration program
416      end
417