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