share/odepack/fortran/opkdemo9.f

   1 c-----------------------------------------------------------------------
   2 c Demonstration program for the DLSODIS package.
   3 c This is the version of 14 June 2001.
   4 c
   5 c This version is in double precision.
   6 c
   7 c This program solves a semi-discretized form of the Burgers equation,
   8 c
   9 c     u  = -(u*u/2)  + eta * u
  10 c      t           x          xx
  11 c
  12 c for  -1 .le. x .le. 1, t .ge. 0.
  13 c Here eta = 0.05.
  14 c Boundary conditions: u(-1,t) = u(1,t) and du/dx(-1,t) = du/dx(1,t).
  15 c Initial profile: square wave
  16 c     u(0,x) = 0    for 1/2 .lt. abs(x) .le. 1
  17 c     u(0,x) = 1/2  for abs(x) = 1/2
  18 c     u(0,x) = 1    for 0 .le. abs(x) .lt. 1/2
  19 c
  20 c An ODE system is generated by a simplified Galerkin treatment
  21 c of the spatial variable x.
  22 c
  23 c Reference:
  24 c R. C. Y. Chin, G. W. Hedstrom, and K. E. Karlsson,
  25 c A Simplified Galerkin Method for Hyperbolic Equations,
  26 c Math. Comp., vol. 33, no. 146 (April 1979), pp. 647-658.
  27 c
  28 c The problem is run with the DLSODIS package with a 12-node mesh,
  29 c for various appropriate values of the method flag mf.
  30 c Output is on unit lout, set to 6 in a data statement below.
  31 c
  32 c Problem specific data:
  33 c npts  = number of unknowns (npts = 0 mod 4)
  34 c nnz   = number of non-zeros in Jacobian before fill in
  35 c nnza  = number of non-zeros in Jacobian after fill in
  36 c lrwk  = length of real work array (taking into account fill in)
  37 c liwk  = length of integer work array
  38 c ipia  = pointer to ia in iw (ia(j) = iw(ipia+j-1)
  39 c ipja  = pointer to ja in iw (ja(j) = iw(ipja+j-1)
  40 c ipic  = pointer to ic in iw array (ic(j) = iw(ipic+j-1))
  41 c ipjc  = pointer to jc in iw array (jc(j) = iw(ipjc+j-1))
  42 c-----------------------------------------------------------------------
  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)
  48 c
  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)
  58 c Pass problem parameters in the Common block test1.
  59       common /test1/ r4d, eodsq
  60 c
  61 c 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/
  69 c
  70       nerr = 0
  71       lrw = lrwk
  72       liw = liwk
  73 c
  74 c 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
  81 c
  82 c 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
  89 c
  90       write(lout,1000)
  91       write(lout,1100) eta,a,b,tinit,tlast,n
  92       write(lout,1200) (y(i), i=1,n)
  93 c
  94 c Set the initial sparse data structures for coefficient matrix A
  95 c and the Jacobian matrix C
  96       call struct(iw(ipia), iw(ipja), iw(ipic), iw(ipjc), n)
  97 c
  98 c The j loop is over error tolerances.
  99       do 200 j = 1,2
 100 c
 101 c 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
 106 c
 107 c 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
 114 c
 115             t = tinit
 116             istate = 0
 117 c
 118             write(lout,1500) itol, rtol(j), atol(j), mf
 119 c
 120 c output loop for each case
 121             do 80 io = 1,nout
 122 c
 123 c 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)
 128 c
 129 c 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)
 138 c
 139 c 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
 150 c
 151       write(lout,6000) nerr
 152       stop
 153 c
 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//)
 162 c
 163  1200 format(/'Initial profile:',/20(6d12.4/))
 164 c
 165  1500 format(///85('*')///'Run with itol =',i2,'  rtol =',d12.2,
 166      1       '  atol =',d12.2,'   mf = ',i3//)
 167 c
 168  2000 format(' Output for time t =',d12.5,'  current h =',d12.5,
 169      1       '  current order =',i2/20(6d12.4/))
 170 c
 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')
 175 c
 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)
 184 c
 185 c end of main program for the DLSODIS demonstration program
 186       end
 187
 188       subroutine struct(ia, ja, ic, jc, n)
 189 c This subroutine computes the initial sparse data structure of
 190 c the mass (ia,ja) and Jacobian (ic,jc) matrices.
 191 c
 192       integer ia(*), ja(*), ic(*), jc(*), n,  jj, k, l, m
 193 c
 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
 211 c
 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
 217 1200  format('Initial sparse data structures'/)
 218 1300  format(' ia  ',15i4/10(5x,15i4/))
 219 1350  format(' ja  ',15i4/10(5x,15i4/))
 220 1400  format(' ic  ',15i4/10(5x,15i4/))
 221 1450  format(' jc  ',15i4/10(5x,15i4/))
 222       end
 223
 224       subroutine res (n, t, y, v, r, ires)
 225 c This subroutine computes the residual vector
 226 c   r = g(t,y) - A(t,y)*v .
 227 c If ires = -1, only g(t,y) is returned in r, since A(t,y) does
 228 c not depend on y.
 229 c No changes need to be made to this routine if n is changed.
 230 c
 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/
 236 c
 237       call gfun (n, t, y, r)
 238       if (ires .eq. -1) return
 239 c
 240       fact1 = one/six
 241       fact4 = four/six
 242 c
 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
 249 c end of subroutine res for the DLSODIS demonstration program
 250       end
 251
 252       subroutine gfun (n, t, y, g)
 253 c This subroutine computes the right-hand side function g(y,t).
 254 c It uses r4d = 1/(4*delta) and eodsq = eta/delta**2
 255 c from the Common block test1.
 256 c
 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/
 261 c
 262       g(1) = r4d*(y(n)**2 - y(2)**2) + eodsq*(y(2) - two*y(1) + y(n))
 263 c
 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
 268 c
 269       g(n) = r4d*(y(n-1)**2 - y(1)**2) + eodsq*(y(1)-two*y(n)+y(n-1))
 270 c
 271       return
 272 c end of subroutine gfun for the DLSODIS demonstration program
 273       end
 274
 275       subroutine addasp (n, t, y, j, ip, jp, pa)
 276 c This subroutine computes the sparse matrix A by columns, adds it to
 277 c pa, and returns the sum in pa.
 278 c The matrix A is periodic tridiagonal, of order n, with nonzero elements
 279 c (reading across) of  1/6, 4/6, 1/6, with 1/6 in the lower left and
 280 c upper right corners.
 281 c
 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/
 285 c
 286 c 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
 293 c
 294 c 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
 299 c end of subroutine addasp for the DLSODIS demonstration program
 300       end
 301
 302       subroutine jacsp (n, t, y, s, j, ip, jp, pdj)
 303 c This subroutine computes the Jacobian dg/dy = d(g-A*s)/dy by
 304 c columns in sparse matrix format.  Only nonzeros are loaded.
 305 c It uses r4d = 1/(4*delta) and eodsq = eta/delta**2 from the Common
 306 c block test1.
 307 c
 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/
 312 c
 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
 319 c
 320       pdj(jm1) = -r2d*y(j) + eodsq
 321       pdj(j) = diag
 322       pdj(jp1) = r2d*y(j) + eodsq
 323       return
 324 c end of subroutine jacsp for the DLSODIS demonstration program
 325       end
 326
 327       double precision function elkup (n, y, ewt, itol, rtol,atol, lout)
 328 c This routine looks up approximately correct values of y at t = 0.4,
 329 c ytrue = y12 or y120 depending on whether n = 12 or 120.
 330 c These were obtained by running with very tight tolerances.
 331 c The returned value is
 332 c   elkup = norm of [ (y - ytrue) / (rtol*abs(ytrue) + atol) ].
 333 c
 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/
 385 c
 386       if ((n-12)*(n-120) .ne. 0) go to 300
 387       if (n .eq. 120) go to 100
 388 c
 389 c Compute local error tolerance using correct y (n = 12).
 390       call dewset(n, itol, rtol, atol, y12, ewt)
 391 c
 392 c 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
 397 c
 398 c Compute local error tolerance using correct y (n = 120).
 399  100  call dewset( n, itol, rtol, atol, y120, ewt )
 400 c
 401 c 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)
 405 c
 406 c Find weighted norm of the error.
 407  200  elkup = dvnorm (n, y, ewt)
 408       return
 409 c
 410 c 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
 415 c end of function elkup for the DLSODIS demonstration program
 416       end