share/odepack/fortran/opkdemo2.f

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