==========Source and Output for LSODE Demonstration Program===================== c----------------------------------------------------------------------- c Demonstration program for the DLSODE package. c This is the version of 14 June 2001. c c This version is in double precision. c c The package is used to solve two simple problems, c one with a full Jacobian, the other with a banded Jacobian, c with all 8 of the appropriate values of mf in each case. c If the errors are too large, or other difficulty occurs, c a warning message is printed. All output is on unit lout = 6. c----------------------------------------------------------------------- external f1, jac1, f2, jac2 integer i, iopar, iopt, iout, istate, itask, itol, iwork, 1 leniw, lenrw, liw, lout, lrw, mband, meth, mf, miter, 2 ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst double precision atol, dtout, er, erm, ero, hu, rtol, rwork, t, 1 tout, tout1, y dimension y(25), rwork(697), iwork(45) data lout/6/, tout1/1.39283880203d0/, dtout/2.214773875d0/ c nerr = 0 itol = 1 rtol = 0.0d0 atol = 1.0d-6 lrw = 697 liw = 45 iopt = 0 c c First problem c neq = 2 nout = 4 write (lout,110) neq,itol,rtol,atol 110 format(/' Demonstration program for DLSODE package'/// 1 ' Problem 1: Van der Pol oscillator:'/ 2 ' xdotdot - 3*(1 - x**2)*xdot + x = 0, ', 3 ' x(0) = 2, xdot(0) = 0'/ 4 ' neq =',i2/ 5 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1//) c do 195 meth = 1,2 do 190 miter = 0,3 mf = 10*meth + miter write (lout,120) mf 120 format(///' Solution with mf =',i3// 1 5x,'t x xdot nq h'//) t = 0.0d0 y(1) = 2.0d0 y(2) = 0.0d0 itask = 1 istate = 1 tout = tout1 ero = 0.0d0 do 170 iout = 1,nout call dlsode(f1,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac1,mf) hu = rwork(11) nqu = iwork(14) write (lout,140) t,y(1),y(2),nqu,hu 140 format(d15.5,d16.5,d14.3,i5,d14.3) if (istate .lt. 0) go to 175 iopar = iout - 2*(iout/2) if (iopar .ne. 0) go to 170 er = abs(y(1))/atol ero = max(ero,er) if (er .gt. 1000.0d0) then write (lout,150) 150 format(//' Warning: error exceeds 1000 * tolerance'//) nerr = nerr + 1 endif 170 tout = tout + dtout 175 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (miter .eq. 2) nfea = nfe - neq*nje if (miter .eq. 3) nfea = nfe - nje write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 180 format(//' Final statistics for this run:'/ 1 ' rwork size =',i4,' iwork size =',i4/ 2 ' number of steps =',i5/ 3 ' number of f-s =',i5/ 4 ' (excluding J-s) =',i5/ 5 ' number of J-s =',i5/ 6 ' error overrun =',d10.2) 190 continue 195 continue c c Second problem c neq = 25 ml = 5 mu = 0 iwork(1) = ml iwork(2) = mu mband = ml + mu + 1 nout = 5 write (lout,210) neq,ml,mu,itol,rtol,atol 210 format(///70('-')/// 1 ' Problem 2: ydot = A * y , where', 2 ' A is a banded lower triangular matrix'/ 3 12x, 'derived from 2-D advection PDE'/ 4 ' neq =',i3,' ml =',i2,' mu =',i2/ 5 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1//) do 295 meth = 1,2 do 290 miter = 0,5 if (miter .eq. 1 .or. miter .eq. 2) go to 290 mf = 10*meth + miter write (lout,220) mf 220 format(///' Solution with mf =',i3// 1 5x,'t max.err. nq h'//) t = 0.0d0 do 230 i = 2,neq 230 y(i) = 0.0d0 y(1) = 1.0d0 itask = 1 istate = 1 tout = 0.01d0 ero = 0.0d0 do 270 iout = 1,nout call dlsode(f2,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac2,mf) call edit2(y,t,erm) hu = rwork(11) nqu = iwork(14) write (lout,240) t,erm,nqu,hu 240 format(d15.5,d14.3,i5,d14.3) if (istate .lt. 0) go to 275 er = erm/atol ero = max(ero,er) if (er .gt. 1000.0d0) then write (lout,150) nerr = nerr + 1 endif 270 tout = tout*10.0d0 275 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (miter .eq. 5) nfea = nfe - mband*nje if (miter .eq. 3) nfea = nfe - nje write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 290 continue 295 continue write (lout,300) nerr 300 format(////' Number of errors encountered =',i3) stop end subroutine f1 (neq, t, y, ydot) integer neq double precision t, y, ydot dimension y(neq), ydot(neq) ydot(1) = y(2) ydot(2) = 3.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1) return end subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd double precision t, y, pd dimension y(neq), pd(nrowpd,neq) pd(1,1) = 0.0d0 pd(1,2) = 1.0d0 pd(2,1) = -6.0d0*y(1)*y(2) - 1.0d0 pd(2,2) = 3.0d0*(1.0d0 - y(1)*y(1)) return end subroutine f2 (neq, t, y, ydot) integer neq, i, j, k, ng double precision t, y, ydot, alph1, alph2, d dimension y(neq), ydot(neq) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ do 10 j = 1,ng do 10 i = 1,ng k = i + (j - 1)*ng d = -2.0d0*y(k) if (i .ne. 1) d = d + y(k-1)*alph1 if (j .ne. 1) d = d + y(k-ng)*alph2 10 ydot(k) = d return end subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng double precision t, y, pd, alph1, alph2 dimension y(neq), pd(nrowpd,neq) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ mband = ml + mu + 1 mu1 = mu + 1 mu2 = mu + 2 do 10 j = 1,neq pd(mu1,j) = -2.0d0 pd(mu2,j) = alph1 10 pd(mband,j) = alph2 do 20 j = ng,neq,ng 20 pd(mu2,j) = 0.0d0 return end subroutine edit2 (y, t, erm) integer i, j, k, ng double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt dimension y(25) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ erm = 0.0d0 if (t .eq. 0.0d0) return ex = 0.0d0 if (t .le. 30.0d0) ex = exp(-2.0d0*t) a2 = 1.0d0 do 60 j = 1,ng a1 = 1.0d0 do 50 i = 1,ng k = i + (j - 1)*ng yt = t**(i+j-2)*ex*a1*a2 er = abs(y(k)-yt) erm = max(erm,er) a1 = a1*alph1/i 50 continue a2 = a2*alph2/j 60 continue return end ................................................................................ Demonstration program for DLSODE package Problem 1: Van der Pol oscillator: xdotdot - 3*(1 - x**2)*xdot + x = 0, x(0) = 2, xdot(0) = 0 neq = 2 itol = 1 rtol = 0.0E+00 atol = 0.1E-05 Solution with mf = 10 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 3 0.123E+00 0.36076E+01 -0.77986E-04 -0.317E+01 5 0.217E-01 0.58224E+01 -0.16801E+01 0.291E+00 3 0.475E-01 0.80372E+01 0.11669E-03 0.317E+01 5 0.234E-01 Final statistics for this run: rwork size = 52 iwork size = 20 number of steps = 297 number of f-s = 352 (excluding J-s) = 352 number of J-s = 0 error overrun = 0.12E+03 Solution with mf = 11 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.121E+00 0.36076E+01 -0.17732E-04 -0.317E+01 5 0.187E-01 0.58224E+01 -0.16801E+01 0.291E+00 6 0.963E-01 0.80372E+01 0.25894E-04 0.317E+01 5 0.190E-01 Final statistics for this run: rwork size = 58 iwork size = 22 number of steps = 203 number of f-s = 281 (excluding J-s) = 281 number of J-s = 29 error overrun = 0.26E+02 Solution with mf = 12 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.121E+00 0.36076E+01 -0.17732E-04 -0.317E+01 5 0.187E-01 0.58224E+01 -0.16801E+01 0.291E+00 6 0.963E-01 0.80372E+01 0.25894E-04 0.317E+01 5 0.190E-01 Final statistics for this run: rwork size = 58 iwork size = 22 number of steps = 203 number of f-s = 339 (excluding J-s) = 281 number of J-s = 29 error overrun = 0.26E+02 Solution with mf = 13 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.739E-01 0.36076E+01 0.34401E-04 -0.317E+01 6 0.260E-01 0.58224E+01 -0.16801E+01 0.291E+00 4 0.133E+00 0.80372E+01 -0.59053E-04 0.317E+01 5 0.205E-01 Final statistics for this run: rwork size = 56 iwork size = 20 number of steps = 198 number of f-s = 315 (excluding J-s) = 289 number of J-s = 26 error overrun = 0.59E+02 Solution with mf = 20 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.549E-01 0.36076E+01 -0.56579E-04 -0.317E+01 5 0.143E-01 0.58224E+01 -0.16801E+01 0.291E+00 4 0.583E-01 0.80372E+01 0.10387E-03 0.317E+01 5 0.149E-01 Final statistics for this run: rwork size = 38 iwork size = 20 number of steps = 289 number of f-s = 321 (excluding J-s) = 321 number of J-s = 0 error overrun = 0.10E+03 Solution with mf = 21 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.676E-01 0.36076E+01 -0.48977E-04 -0.317E+01 5 0.141E-01 0.58224E+01 -0.16801E+01 0.291E+00 5 0.126E+00 0.80372E+01 0.96867E-04 0.317E+01 5 0.142E-01 Final statistics for this run: rwork size = 44 iwork size = 22 number of steps = 262 number of f-s = 345 (excluding J-s) = 345 number of J-s = 30 error overrun = 0.97E+02 Solution with mf = 22 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.676E-01 0.36076E+01 -0.48977E-04 -0.317E+01 5 0.141E-01 0.58224E+01 -0.16801E+01 0.291E+00 5 0.126E+00 0.80372E+01 0.96867E-04 0.317E+01 5 0.142E-01 Final statistics for this run: rwork size = 44 iwork size = 22 number of steps = 262 number of f-s = 405 (excluding J-s) = 345 number of J-s = 30 error overrun = 0.97E+02 Solution with mf = 23 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.709E-01 0.36076E+01 -0.46705E-04 -0.317E+01 5 0.139E-01 0.58224E+01 -0.16801E+01 0.291E+00 3 0.719E-01 0.80372E+01 0.54700E-04 0.317E+01 5 0.154E-01 Final statistics for this run: rwork size = 42 iwork size = 20 number of steps = 271 number of f-s = 414 (excluding J-s) = 383 number of J-s = 31 error overrun = 0.55E+02 ---------------------------------------------------------------------- Problem 2: ydot = A * y , where A is a banded lower triangular matrix derived from 2-D advection PDE neq = 25 ml = 5 mu = 0 itol = 1 rtol = 0.0E+00 atol = 0.1E-05 Solution with mf = 10 t max.err. nq h 0.10000E-01 0.556E-06 2 0.766E-02 0.10000E+00 0.655E-05 3 0.249E-01 0.10000E+01 0.274E-05 4 0.520E-01 0.10000E+02 0.114E-05 3 0.117E+00 0.10000E+03 0.221E-05 2 0.262E+00 Final statistics for this run: rwork size = 420 iwork size = 20 number of steps = 524 number of f-s = 552 (excluding J-s) = 552 number of J-s = 0 error overrun = 0.65E+01 Solution with mf = 13 t max.err. nq h 0.10000E-01 0.839E-06 2 0.949E-02 0.10000E+00 0.208E-05 3 0.250E-01 0.10000E+01 0.127E-03 3 0.168E-01 0.10000E+02 0.113E-04 3 0.385E+00 0.10000E+03 0.145E-05 2 0.149E+02 Final statistics for this run: rwork size = 447 iwork size = 20 number of steps = 129 number of f-s = 235 (excluding J-s) = 201 number of J-s = 34 error overrun = 0.13E+03 Solution with mf = 14 t max.err. nq h 0.10000E-01 0.877E-06 2 0.965E-02 0.10000E+00 0.206E-05 3 0.250E-01 0.10000E+01 0.126E-05 5 0.935E-01 0.10000E+02 0.311E-06 6 0.442E+00 0.10000E+03 0.159E-07 2 0.291E+02 Final statistics for this run: rwork size = 697 iwork size = 45 number of steps = 92 number of f-s = 113 (excluding J-s) = 113 number of J-s = 18 error overrun = 0.21E+01 Solution with mf = 15 t max.err. nq h 0.10000E-01 0.877E-06 2 0.965E-02 0.10000E+00 0.206E-05 3 0.250E-01 0.10000E+01 0.126E-05 5 0.935E-01 0.10000E+02 0.311E-06 6 0.442E+00 0.10000E+03 0.160E-07 2 0.291E+02 Final statistics for this run: rwork size = 697 iwork size = 45 number of steps = 92 number of f-s = 221 (excluding J-s) = 113 number of J-s = 18 error overrun = 0.21E+01 Solution with mf = 20 t max.err. nq h 0.10000E-01 0.465E-06 2 0.483E-02 0.10000E+00 0.131E-05 3 0.148E-01 0.10000E+01 0.427E-05 5 0.635E-01 0.10000E+02 0.192E-05 4 0.351E+00 0.10000E+03 0.929E-07 1 0.455E+00 Final statistics for this run: rwork size = 245 iwork size = 20 number of steps = 330 number of f-s = 530 (excluding J-s) = 530 number of J-s = 0 error overrun = 0.43E+01 Solution with mf = 23 t max.err. nq h 0.10000E-01 0.101E-05 2 0.598E-02 0.10000E+00 0.446E-06 3 0.146E-01 0.10000E+01 0.153E-05 5 0.738E-01 0.10000E+02 0.578E-06 4 0.324E+00 0.10000E+03 0.908E-08 1 0.992E+02 Final statistics for this run: rwork size = 272 iwork size = 20 number of steps = 180 number of f-s = 325 (excluding J-s) = 274 number of J-s = 51 error overrun = 0.15E+01 Solution with mf = 24 t max.err. nq h 0.10000E-01 0.104E-05 2 0.608E-02 0.10000E+00 0.463E-06 3 0.146E-01 0.10000E+01 0.247E-05 5 0.666E-01 0.10000E+02 0.828E-06 5 0.391E+00 0.10000E+03 0.384E-09 1 0.108E+03 Final statistics for this run: rwork size = 522 iwork size = 45 number of steps = 118 number of f-s = 136 (excluding J-s) = 136 number of J-s = 18 error overrun = 0.25E+01 Solution with mf = 25 t max.err. nq h 0.10000E-01 0.104E-05 2 0.608E-02 0.10000E+00 0.463E-06 3 0.146E-01 0.10000E+01 0.247E-05 5 0.666E-01 0.10000E+02 0.828E-06 5 0.391E+00 0.10000E+03 0.384E-09 1 0.108E+03 Final statistics for this run: rwork size = 522 iwork size = 45 number of steps = 118 number of f-s = 244 (excluding J-s) = 136 number of J-s = 18 error overrun = 0.25E+01 Number of errors encountered = 0 ==========Source and Output for LSODES Demonstration Program==================== c----------------------------------------------------------------------- c Demonstration program for the DLSODES package. c This is the version of 14 June 2001. c c This version is in double precision. c c The package is used for each of the relevant values of mf to solve c the problem ydot = A * y, where A is the 9 by 9 sparse matrix c c -4 1 1 c 1 -4 1 1 c 1 -4 1 c -4 1 1 c A = 1 -4 1 1 c 1 -4 1 c -4 1 c 1 -4 1 c 1 -4 c c The initial conditions are y(0) = (1, 2, 3, ..., 9). c Output is printed at t = 1, 2, and 3. c Each case is solved first with nominal (large) values of lrw and liw, c and then with values given by lenrw and leniw (optional outputs) c on the first run, as a check on these computed work array lengths. c If the errors are too large, or other difficulty occurs, c a warning message is printed. c All output is on unit lout, which is data-loaded to 6 below. c----------------------------------------------------------------------- external fdem, jdem integer i, ia, igrid, iopt, iout, irun, istate, itask, itol, 1 iwork, j, ja, k, l, leniw, lenrw, liw, lout, lrw, 2 m, meth, mf, miter, moss, neq, nerr, nfe, nfea, 3 ngp, nje, nlu, nnz, nout, nqu, nst, nzl, nzu double precision atol, erm, ero, hu, rtol, rwork, t, tout, y dimension y(9), ia(10), ja(50), iwork(90), rwork(1000) equivalence (ia(1),iwork(31)), (ja(1),iwork(41)) data lout/6/ c c Write heading and set fixed parameters. write(lout,10) 10 format(/'Demonstration problem for the DLSODES package'//) nerr = 0 igrid = 3 neq = igrid**2 t = 0.0d0 itol = 1 rtol = 0.0d0 atol = 1.0d-5 itask = 1 iopt = 0 do 20 i = 1,neq 20 y(i) = i ia(1) = 1 k = 1 do 60 m = 1,igrid do 50 l = 1,igrid j = l + (m - 1)*igrid if (m .gt. 1) then ja(k) = j - igrid k = k + 1 endif 30 if (l .gt. 1) then ja(k) = j - 1 k = k + 1 endif 35 ja(k) = j k = k + 1 if (l .lt. igrid) then ja(k) = j + 1 k = k + 1 endif 40 ia(j+1) = k 50 continue 60 continue write (lout,80)neq,t,rtol,atol,(y(i),i=1,neq) 80 format(' neq =',i4,5x,'t0 =',f4.1,5x,'rtol =',d12.3,5x, 1 'atol =',d12.3//' Initial y vector = ',9f5.1) c c Loop over all relevant values of mf. do 193 moss = 0,2 do 192 meth = 1,2 do 191 miter = 0,3 if ( (miter.eq.0 .or. miter.eq.3) .and. moss.ne.0) go to 191 mf = 100*moss + 10*meth + miter write (lout,100) 100 format(//80('*')) c First run: nominal work array lengths, 3 output points. irun = 1 lrw = 1000 liw = 90 nout = 3 110 continue write (lout,120)mf,lrw,liw 120 format(//'Run with mf =',i4,'.',5x, 1 'Input work lengths lrw, liw =',2i6/) do 125 i = 1,neq 125 y(i) = i t = 0.0d0 tout = 1.0d0 istate = 1 ero = 0.0d0 c Loop over output points. Do output and accuracy check at each. do 170 iout = 1,nout call dlsodes (fdem, neq, y, t, tout, itol, rtol, atol, itask, 1 istate, iopt, rwork, lrw, iwork, liw, jdem, mf) nst = iwork(11) hu = rwork(11) nqu = iwork(14) call edit (y, iout, erm) write(lout,140)t,nst,hu,nqu,erm,(y(i),i=1,neq) 140 format('At t =',f5.1,3x,'nst =',i4,3x,'hu =',d12.3,3x, 1 'nqu =',i3,3x,' max. err. =',d11.3/ 2 ' y array = ',4d15.6/5d15.6) if (istate .lt. 0) go to 175 erm = erm/atol ero = max(ero,erm) if (erm .gt. 100.0d0) then write (lout,150) 150 format(//' Warning: error exceeds 100 * tolerance'//) nerr = nerr + 1 endif tout = tout + 1.0d0 170 continue 175 continue if (istate .lt. 0) nerr = nerr + 1 if (irun .eq. 2) go to 191 c Print final statistics (first run only) nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nnz = iwork(19) ngp = iwork(20) nlu = iwork(21) nzl = iwork(25) nzu = iwork(26) nfea = nfe if (miter .eq. 2) nfea = nfe - ngp*nje if (miter .eq. 3) nfea = nfe - nje write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 180 format(/'Final statistics for this run:'/ 1 ' rwork size =',i4,' iwork size =',i4/ 2 ' number of steps =',i5/ 3 ' number of f-s =',i5/ 4 ' (excluding J-s) =',i5/ 5 ' number of J-s =',i5/ 6 ' error overrun =',d10.2) if (miter .eq. 1 .or. miter .eq. 2) 1 write (lout,185)nnz,ngp,nlu,nzl,nzu 185 format(' number of nonzeros in J = ',i5/ 1 ' number of J index groups =',i5/ 2 ' number of LU decomp-s =',i5/ 3 ' nonzeros in strict lower factor =',i5/ 4 ' nonzeros in strict upper factor =',i5) if (istate .lt. 0) go to 191 if (miter .eq. 1 .or. miter .eq. 2) 1 call ssout (neq, rwork(21), iwork, lout) c Return for second run: minimal work array lengths, 1 output point. irun = irun + 1 lrw = lenrw liw = leniw nout = 1 go to 110 191 continue 192 continue 193 continue c write (lout,100) write (lout,200) nerr 200 format(//'Number of errors encountered =',i3) stop end subroutine fdem (neq, t, y, ydot) integer neq, i, igrid, j, l, m double precision t, y, ydot dimension y(neq), ydot(neq) data igrid/3/ do 5 i = 1,neq 5 ydot(i) = 0.0d0 do 20 m = 1,igrid do 10 l = 1,igrid j = l + (m - 1)*igrid if (m .ne. 1) ydot(j-igrid) = ydot(j-igrid) + y(j) if (l .ne. 1) ydot(j-1) = ydot(j-1) + y(j) ydot(j) = ydot(j) - 4.0d0*y(j) if (l .ne. igrid) ydot(j+1) = ydot(j+1) + y(j) 10 continue 20 continue return end subroutine jdem (neq, t, y, j, ia, ja, pdj) integer neq, j, ia, ja, igrid, l, m double precision t, y, pdj dimension y(neq), ia(*), ja(*), pdj(neq) data igrid/3/ m = (j - 1)/igrid + 1 l = j - (m - 1)*igrid pdj(j) = -4.0d0 if (m .ne. 1) pdj(j-igrid) = 1.0d0 if (l .ne. 1) pdj(j-1) = 1.0d0 if (l .ne. igrid) pdj(j+1) = 1.0d0 return end subroutine edit (y, iout, erm) integer iout, i, neq double precision y, erm, er, yex dimension y(*),yex(9,3) data neq /9/ data yex /6.687279d-01, 9.901910d-01, 7.603061d-01, 1 8.077979d-01, 1.170226e+00, 8.810605d-01, 5.013331d-01, 2 7.201389d-01, 5.379644d-01, 1.340488d-01, 1.917157d-01, 3 1.374034d-01, 1.007882d-01, 1.437868d-01, 1.028010d-01, 4 3.844343d-02, 5.477593d-02, 3.911435d-02, 1.929166d-02, 5 2.735444d-02, 1.939611d-02, 1.055981d-02, 1.496753d-02, 6 1.060897d-02, 2.913689d-03, 4.128975d-03, 2.925977d-03/ erm = 0.0d0 do 10 i = 1,neq er = abs(y(i) - yex(i,iout)) 10 erm = max(erm,er) return end subroutine ssout (neq, iwk, iwork, lout) integer neq, iwk, iwork, lout integer i, i1, i2, ipian, ipjan, nnz dimension iwk(*), iwork(*) ipian = iwork(23) ipjan = iwork(24) nnz = iwork(19) i1 = ipian i2 = i1 + neq write (lout,10)(iwk(i),i=i1,i2) 10 format(/' structure descriptor array ian ='/(20i4)) i1 = ipjan i2 = i1 + nnz - 1 write (lout,20)(iwk(i),i=i1,i2) 20 format(/' structure descriptor array jan ='/(20i4)) return end ................................................................................ Demonstration problem for the DLSODES package neq = 9 t0 = 0.0 rtol = 0.000E+00 atol = 0.100E-04 Initial y vector = 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 ******************************************************************************** Run with mf = 10. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 42 hu = 0.241E-01 nqu = 4 max. err. = 0.394E-05 y array = 0.668730E+00 0.990188E+00 0.760308E+00 0.807799E+00 0.117023E+01 0.881061E+00 0.501332E+00 0.720142E+00 0.537962E+00 At t = 2.0 nst = 71 hu = 0.680E-01 nqu = 3 max. err. = 0.273E-04 y array = 0.134047E+00 0.191717E+00 0.137407E+00 0.100802E+00 0.143808E+00 0.102820E+00 0.384623E-01 0.548033E-01 0.391361E-01 At t = 3.0 nst = 90 hu = 0.455E-01 nqu = 3 max. err. = 0.121E-04 y array = 0.193008E-01 0.273568E-01 0.194059E-01 0.105663E-01 0.149796E-01 0.106158E-01 0.291803E-02 0.413489E-02 0.293048E-02 Final statistics for this run: rwork size = 164 iwork size = 30 number of steps = 90 number of f-s = 98 (excluding J-s) = 98 number of J-s = 0 error overrun = 0.27E+01 Run with mf = 10. Input work lengths lrw, liw = 164 30 At t = 1.0 nst = 42 hu = 0.241E-01 nqu = 4 max. err. = 0.394E-05 y array = 0.668730E+00 0.990188E+00 0.760308E+00 0.807799E+00 0.117023E+01 0.881061E+00 0.501332E+00 0.720142E+00 0.537962E+00 ******************************************************************************** Run with mf = 11. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 At t = 2.0 nst = 37 hu = 0.115E+00 nqu = 5 max. err. = 0.593E-05 y array = 0.134047E+00 0.191713E+00 0.137403E+00 0.100789E+00 0.143787E+00 0.102803E+00 0.384477E-01 0.547819E-01 0.391199E-01 At t = 3.0 nst = 44 hu = 0.155E+00 nqu = 5 max. err. = 0.386E-05 y array = 0.192920E-01 0.273552E-01 0.193970E-01 0.105624E-01 0.149714E-01 0.106120E-01 0.291609E-02 0.413247E-02 0.292857E-02 Final statistics for this run: rwork size = 308 iwork size = 67 number of steps = 44 number of f-s = 56 (excluding J-s) = 56 number of J-s = 1 error overrun = 0.99E+00 number of nonzeros in J = 27 number of J index groups = 0 number of LU decomp-s = 9 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 11. Input work lengths lrw, liw = 308 67 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 ******************************************************************************** Run with mf = 12. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 At t = 2.0 nst = 37 hu = 0.115E+00 nqu = 5 max. err. = 0.593E-05 y array = 0.134047E+00 0.191713E+00 0.137403E+00 0.100789E+00 0.143787E+00 0.102803E+00 0.384477E-01 0.547819E-01 0.391199E-01 At t = 3.0 nst = 44 hu = 0.155E+00 nqu = 5 max. err. = 0.386E-05 y array = 0.192920E-01 0.273552E-01 0.193970E-01 0.105624E-01 0.149714E-01 0.106120E-01 0.291609E-02 0.413247E-02 0.292857E-02 Final statistics for this run: rwork size = 315 iwork size = 67 number of steps = 44 number of f-s = 60 (excluding J-s) = 56 number of J-s = 1 error overrun = 0.99E+00 number of nonzeros in J = 27 number of J index groups = 4 number of LU decomp-s = 9 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 12. Input work lengths lrw, liw = 315 67 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 ******************************************************************************** Run with mf = 13. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 28 hu = 0.825E-01 nqu = 5 max. err. = 0.591E-05 y array = 0.668729E+00 0.990192E+00 0.760304E+00 0.807797E+00 0.117022E+01 0.881060E+00 0.501334E+00 0.720136E+00 0.537970E+00 At t = 2.0 nst = 41 hu = 0.825E-01 nqu = 5 max. err. = 0.468E-05 y array = 0.134048E+00 0.191713E+00 0.137405E+00 0.100784E+00 0.143782E+00 0.102798E+00 0.384467E-01 0.547752E-01 0.391129E-01 At t = 3.0 nst = 58 hu = 0.536E-01 nqu = 4 max. err. = 0.801E-04 y array = 0.193021E-01 0.273573E-01 0.194762E-01 0.105617E-01 0.149669E-01 0.106567E-01 0.291302E-02 0.412818E-02 0.292528E-02 Final statistics for this run: rwork size = 175 iwork size = 30 number of steps = 58 number of f-s = 90 (excluding J-s) = 80 number of J-s = 10 error overrun = 0.80E+01 Run with mf = 13. Input work lengths lrw, liw = 175 30 At t = 1.0 nst = 28 hu = 0.825E-01 nqu = 5 max. err. = 0.591E-05 y array = 0.668729E+00 0.990192E+00 0.760304E+00 0.807797E+00 0.117022E+01 0.881060E+00 0.501334E+00 0.720136E+00 0.537970E+00 ******************************************************************************** Run with mf = 20. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 39 hu = 0.549E-01 nqu = 5 max. err. = 0.378E-04 y array = 0.668726E+00 0.990193E+00 0.760309E+00 0.807791E+00 0.117020E+01 0.881056E+00 0.501356E+00 0.720124E+00 0.538002E+00 At t = 2.0 nst = 53 hu = 0.677E-01 nqu = 5 max. err. = 0.113E-04 y array = 0.134039E+00 0.191719E+00 0.137397E+00 0.100792E+00 0.143779E+00 0.102808E+00 0.384485E-01 0.547872E-01 0.391222E-01 At t = 3.0 nst = 64 hu = 0.123E+00 nqu = 5 max. err. = 0.869E-05 y array = 0.192944E-01 0.273518E-01 0.193999E-01 0.105634E-01 0.149762E-01 0.106132E-01 0.291807E-02 0.413507E-02 0.293054E-02 Final statistics for this run: rwork size = 101 iwork size = 30 number of steps = 64 number of f-s = 77 (excluding J-s) = 77 number of J-s = 0 error overrun = 0.38E+01 Run with mf = 20. Input work lengths lrw, liw = 101 30 At t = 1.0 nst = 39 hu = 0.549E-01 nqu = 5 max. err. = 0.378E-04 y array = 0.668726E+00 0.990193E+00 0.760309E+00 0.807791E+00 0.117020E+01 0.881056E+00 0.501356E+00 0.720124E+00 0.538002E+00 ******************************************************************************** Run with mf = 21. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 At t = 2.0 nst = 52 hu = 0.105E+00 nqu = 5 max. err. = 0.132E-04 y array = 0.134044E+00 0.191709E+00 0.137402E+00 0.100788E+00 0.143786E+00 0.102805E+00 0.384531E-01 0.547890E-01 0.391276E-01 At t = 3.0 nst = 61 hu = 0.132E+00 nqu = 5 max. err. = 0.134E-04 y array = 0.192907E-01 0.273543E-01 0.193977E-01 0.105672E-01 0.149788E-01 0.106186E-01 0.292280E-02 0.414233E-02 0.293619E-02 Final statistics for this run: rwork size = 245 iwork size = 67 number of steps = 61 number of f-s = 71 (excluding J-s) = 71 number of J-s = 2 error overrun = 0.25E+01 number of nonzeros in J = 27 number of J index groups = 0 number of LU decomp-s = 8 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 21. Input work lengths lrw, liw = 245 67 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 ******************************************************************************** Run with mf = 22. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 At t = 2.0 nst = 52 hu = 0.105E+00 nqu = 5 max. err. = 0.132E-04 y array = 0.134044E+00 0.191709E+00 0.137402E+00 0.100788E+00 0.143786E+00 0.102805E+00 0.384531E-01 0.547890E-01 0.391276E-01 At t = 3.0 nst = 61 hu = 0.132E+00 nqu = 5 max. err. = 0.134E-04 y array = 0.192907E-01 0.273543E-01 0.193977E-01 0.105672E-01 0.149788E-01 0.106186E-01 0.292280E-02 0.414233E-02 0.293619E-02 Final statistics for this run: rwork size = 252 iwork size = 67 number of steps = 61 number of f-s = 79 (excluding J-s) = 71 number of J-s = 2 error overrun = 0.25E+01 number of nonzeros in J = 27 number of J index groups = 4 number of LU decomp-s = 8 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 22. Input work lengths lrw, liw = 252 67 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 ******************************************************************************** Run with mf = 23. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 39 hu = 0.564E-01 nqu = 5 max. err. = 0.334E-04 y array = 0.668727E+00 0.990191E+00 0.760307E+00 0.807793E+00 0.117019E+01 0.881055E+00 0.501339E+00 0.720129E+00 0.537977E+00 At t = 2.0 nst = 53 hu = 0.720E-01 nqu = 5 max. err. = 0.268E-04 y array = 0.134056E+00 0.191713E+00 0.137402E+00 0.100785E+00 0.143781E+00 0.102798E+00 0.384437E-01 0.548027E-01 0.391161E-01 At t = 3.0 nst = 80 hu = 0.200E-01 nqu = 4 max. err. = 0.490E-04 y array = 0.192918E-01 0.273431E-01 0.194010E-01 0.105628E-01 0.149652E-01 0.106029E-01 0.290678E-02 0.412271E-02 0.297494E-02 Final statistics for this run: rwork size = 112 iwork size = 30 number of steps = 80 number of f-s = 137 (excluding J-s) = 122 number of J-s = 15 error overrun = 0.49E+01 Run with mf = 23. Input work lengths lrw, liw = 112 30 At t = 1.0 nst = 39 hu = 0.564E-01 nqu = 5 max. err. = 0.334E-04 y array = 0.668727E+00 0.990191E+00 0.760307E+00 0.807793E+00 0.117019E+01 0.881055E+00 0.501339E+00 0.720129E+00 0.537977E+00 ******************************************************************************** Run with mf = 111. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 At t = 2.0 nst = 37 hu = 0.115E+00 nqu = 5 max. err. = 0.593E-05 y array = 0.134047E+00 0.191713E+00 0.137403E+00 0.100789E+00 0.143787E+00 0.102803E+00 0.384477E-01 0.547819E-01 0.391199E-01 At t = 3.0 nst = 44 hu = 0.155E+00 nqu = 5 max. err. = 0.386E-05 y array = 0.192920E-01 0.273552E-01 0.193970E-01 0.105624E-01 0.149714E-01 0.106120E-01 0.291609E-02 0.413247E-02 0.292857E-02 Final statistics for this run: rwork size = 308 iwork size = 30 number of steps = 44 number of f-s = 56 (excluding J-s) = 56 number of J-s = 1 error overrun = 0.99E+00 number of nonzeros in J = 27 number of J index groups = 0 number of LU decomp-s = 9 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 111. Input work lengths lrw, liw = 308 30 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 ******************************************************************************** Run with mf = 112. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 At t = 2.0 nst = 37 hu = 0.115E+00 nqu = 5 max. err. = 0.593E-05 y array = 0.134047E+00 0.191713E+00 0.137403E+00 0.100789E+00 0.143787E+00 0.102803E+00 0.384477E-01 0.547819E-01 0.391199E-01 At t = 3.0 nst = 44 hu = 0.155E+00 nqu = 5 max. err. = 0.386E-05 y array = 0.192920E-01 0.273552E-01 0.193970E-01 0.105624E-01 0.149714E-01 0.106120E-01 0.291609E-02 0.413247E-02 0.292857E-02 Final statistics for this run: rwork size = 315 iwork size = 30 number of steps = 44 number of f-s = 60 (excluding J-s) = 56 number of J-s = 1 error overrun = 0.99E+00 number of nonzeros in J = 27 number of J index groups = 4 number of LU decomp-s = 9 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 112. Input work lengths lrw, liw = 315 30 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 ******************************************************************************** Run with mf = 121. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 At t = 2.0 nst = 52 hu = 0.105E+00 nqu = 5 max. err. = 0.132E-04 y array = 0.134044E+00 0.191709E+00 0.137402E+00 0.100788E+00 0.143786E+00 0.102805E+00 0.384531E-01 0.547890E-01 0.391276E-01 At t = 3.0 nst = 61 hu = 0.132E+00 nqu = 5 max. err. = 0.134E-04 y array = 0.192907E-01 0.273543E-01 0.193977E-01 0.105672E-01 0.149788E-01 0.106186E-01 0.292280E-02 0.414233E-02 0.293619E-02 Final statistics for this run: rwork size = 245 iwork size = 30 number of steps = 61 number of f-s = 71 (excluding J-s) = 71 number of J-s = 2 error overrun = 0.25E+01 number of nonzeros in J = 27 number of J index groups = 0 number of LU decomp-s = 8 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 121. Input work lengths lrw, liw = 245 30 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 ******************************************************************************** Run with mf = 122. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 At t = 2.0 nst = 52 hu = 0.105E+00 nqu = 5 max. err. = 0.132E-04 y array = 0.134044E+00 0.191709E+00 0.137402E+00 0.100788E+00 0.143786E+00 0.102805E+00 0.384531E-01 0.547890E-01 0.391276E-01 At t = 3.0 nst = 61 hu = 0.132E+00 nqu = 5 max. err. = 0.134E-04 y array = 0.192907E-01 0.273543E-01 0.193977E-01 0.105672E-01 0.149788E-01 0.106186E-01 0.292280E-02 0.414233E-02 0.293619E-02 Final statistics for this run: rwork size = 252 iwork size = 30 number of steps = 61 number of f-s = 79 (excluding J-s) = 71 number of J-s = 2 error overrun = 0.25E+01 number of nonzeros in J = 27 number of J index groups = 4 number of LU decomp-s = 8 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 122. Input work lengths lrw, liw = 252 30 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 ******************************************************************************** Run with mf = 211. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 At t = 2.0 nst = 37 hu = 0.115E+00 nqu = 5 max. err. = 0.593E-05 y array = 0.134047E+00 0.191713E+00 0.137403E+00 0.100789E+00 0.143787E+00 0.102803E+00 0.384477E-01 0.547819E-01 0.391199E-01 At t = 3.0 nst = 44 hu = 0.155E+00 nqu = 5 max. err. = 0.386E-05 y array = 0.192920E-01 0.273552E-01 0.193970E-01 0.105624E-01 0.149714E-01 0.106120E-01 0.291609E-02 0.413247E-02 0.292857E-02 Final statistics for this run: rwork size = 308 iwork size = 30 number of steps = 44 number of f-s = 56 (excluding J-s) = 56 number of J-s = 1 error overrun = 0.99E+00 number of nonzeros in J = 27 number of J index groups = 0 number of LU decomp-s = 9 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 211. Input work lengths lrw, liw = 308 30 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 ******************************************************************************** Run with mf = 212. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 At t = 2.0 nst = 37 hu = 0.115E+00 nqu = 5 max. err. = 0.593E-05 y array = 0.134047E+00 0.191713E+00 0.137403E+00 0.100789E+00 0.143787E+00 0.102803E+00 0.384477E-01 0.547819E-01 0.391199E-01 At t = 3.0 nst = 44 hu = 0.155E+00 nqu = 5 max. err. = 0.386E-05 y array = 0.192920E-01 0.273552E-01 0.193970E-01 0.105624E-01 0.149714E-01 0.106120E-01 0.291609E-02 0.413247E-02 0.292857E-02 Final statistics for this run: rwork size = 315 iwork size = 30 number of steps = 44 number of f-s = 60 (excluding J-s) = 56 number of J-s = 1 error overrun = 0.99E+00 number of nonzeros in J = 27 number of J index groups = 4 number of LU decomp-s = 9 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 212. Input work lengths lrw, liw = 315 30 At t = 1.0 nst = 27 hu = 0.862E-01 nqu = 5 max. err. = 0.988E-05 y array = 0.668727E+00 0.990190E+00 0.760307E+00 0.807796E+00 0.117022E+01 0.881060E+00 0.501339E+00 0.720139E+00 0.537974E+00 ******************************************************************************** Run with mf = 221. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 At t = 2.0 nst = 52 hu = 0.105E+00 nqu = 5 max. err. = 0.132E-04 y array = 0.134044E+00 0.191709E+00 0.137402E+00 0.100788E+00 0.143786E+00 0.102805E+00 0.384531E-01 0.547890E-01 0.391276E-01 At t = 3.0 nst = 61 hu = 0.132E+00 nqu = 5 max. err. = 0.134E-04 y array = 0.192907E-01 0.273543E-01 0.193977E-01 0.105672E-01 0.149788E-01 0.106186E-01 0.292280E-02 0.414233E-02 0.293619E-02 Final statistics for this run: rwork size = 245 iwork size = 30 number of steps = 61 number of f-s = 71 (excluding J-s) = 71 number of J-s = 2 error overrun = 0.25E+01 number of nonzeros in J = 27 number of J index groups = 0 number of LU decomp-s = 8 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 221. Input work lengths lrw, liw = 245 30 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 ******************************************************************************** Run with mf = 222. Input work lengths lrw, liw = 1000 90 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 At t = 2.0 nst = 52 hu = 0.105E+00 nqu = 5 max. err. = 0.132E-04 y array = 0.134044E+00 0.191709E+00 0.137402E+00 0.100788E+00 0.143786E+00 0.102805E+00 0.384531E-01 0.547890E-01 0.391276E-01 At t = 3.0 nst = 61 hu = 0.132E+00 nqu = 5 max. err. = 0.134E-04 y array = 0.192907E-01 0.273543E-01 0.193977E-01 0.105672E-01 0.149788E-01 0.106186E-01 0.292280E-02 0.414233E-02 0.293619E-02 Final statistics for this run: rwork size = 252 iwork size = 30 number of steps = 61 number of f-s = 79 (excluding J-s) = 71 number of J-s = 2 error overrun = 0.25E+01 number of nonzeros in J = 27 number of J index groups = 4 number of LU decomp-s = 8 nonzeros in strict lower factor = 8 nonzeros in strict upper factor = 14 structure descriptor array ian = 1 3 6 8 11 15 18 21 25 28 structure descriptor array jan = 1 2 3 1 2 3 2 1 5 4 2 6 5 4 3 6 5 7 4 8 9 7 5 8 9 6 8 Run with mf = 222. Input work lengths lrw, liw = 252 30 At t = 1.0 nst = 38 hu = 0.573E-01 nqu = 5 max. err. = 0.254E-04 y array = 0.668726E+00 0.990191E+00 0.760308E+00 0.807793E+00 0.117021E+01 0.881059E+00 0.501348E+00 0.720133E+00 0.537990E+00 ******************************************************************************** Number of errors encountered = 0 ==========Source and Output for LSODA Demonstration Program===================== c----------------------------------------------------------------------- c Demonstration program for the DLSODA package. c This is the version of 14 June 2001. c c This version is in double precision. c c The package is used to solve two simple problems, c one with a full Jacobian, the other with a banded Jacobian, c with the 2 appropriate values of jt in each case. c If the errors are too large, or other difficulty occurs, c a warning message is printed. All output is on unit lout = 6. c----------------------------------------------------------------------- external f1, jac1, f2, jac2 integer i, iopar, iopt, iout, istate, itask, itol, iwork, 1 jt, leniw, lenrw, liw, lout, lrw, mband, mused, 2 ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst double precision atol, dtout, dtout0, dtout1, er, erm, ero, hu, 1 rtol, rwork, t, tout, tout1, tsw, y dimension y(25), rwork(522), iwork(45) data lout/6/, tout1/16.921743d0/, dtout/17.341162d0/ c nerr = 0 itol = 1 rtol = 0.0d0 atol = 1.0d-8 lrw = 522 liw = 45 iopt = 0 c c First problem c neq = 2 nout = 4 write (lout,110) neq,itol,rtol,atol 110 format(/'Demonstration program for DLSODA package'//// 1 ' Problem 1: Van der Pol oscillator:'/ 2 ' xdotdot - 20*(1 - x**2)*xdot + x = 0, ', 3 ' x(0) = 2, xdot(0) = 0'/' neq =',i2/ 4 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1//) c do 190 jt = 1,2 write (lout,120) jt 120 format(//' Solution with jt =',i3// 1 ' t x xdot meth', 2 ' nq h tsw'//) t = 0.0d0 y(1) = 2.0d0 y(2) = 0.0d0 itask = 1 istate = 1 dtout0 = 0.5d0*tout1 dtout1 = 0.5d0*dtout tout = dtout0 ero = 0.0d0 do 170 iout = 1,nout call dlsoda(f1,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac1,jt) hu = rwork(11) tsw = rwork(15) nqu = iwork(14) mused = iwork(19) write (lout,140) t,y(1),y(2),mused,nqu,hu,tsw 140 format(d12.5,d16.5,d14.3,2i6,2d13.3) if (istate .lt. 0) go to 175 iopar = iout - 2*(iout/2) if (iopar .ne. 0) go to 160 er = abs(y(1)) ero = max(ero,er) if (er .gt. 1.0d-2) then write (lout,150) 150 format(//' Warning: value at root exceeds 1.0d-2'//) nerr = nerr + 1 endif 160 if (iout .eq. 1) tout = tout + dtout0 if (iout .gt. 1) tout = tout + dtout1 170 continue 175 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (jt .eq. 2) nfea = nfe - neq*nje write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 180 format(//' Final statistics for this run:'/ 1 ' rwork size =',i4,' iwork size =',i4/ 2 ' number of steps =',i5/ 3 ' number of f-s =',i5/ 4 ' (excluding J-s) =',i5/ 5 ' number of J-s =',i5/ 6 ' max. error at root =',d10.2) 190 continue c c Second problem c neq = 25 ml = 5 mu = 0 iwork(1) = ml iwork(2) = mu mband = ml + mu + 1 atol = 1.0d-6 nout = 5 write (lout,210) neq,ml,mu,itol,rtol,atol 210 format(///80('-')/// 1 ' Problem 2: ydot = A * y , where', 2 ' A is a banded lower triangular matrix'/ 2 ' derived from 2-D advection PDE'/ 3 ' neq =',i3,' ml =',i2,' mu =',i2/ 4 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1//) do 290 jt = 4,5 write (lout,220) jt 220 format(//' Solution with jt =',i3// 1 ' t max.err. meth ', 2 'nq h tsw'//) t = 0.0d0 do 230 i = 2,neq 230 y(i) = 0.0d0 y(1) = 1.0d0 itask = 1 istate = 1 tout = 0.01d0 ero = 0.0d0 do 270 iout = 1,nout call dlsoda(f2,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac2,jt) call edit2(y,t,erm) hu = rwork(11) tsw = rwork(15) nqu = iwork(14) mused = iwork(19) write (lout,240) t,erm,mused,nqu,hu,tsw 240 format(d15.5,d14.3,2i6,2d14.3) if (istate .lt. 0) go to 275 er = erm/atol ero = max(ero,er) if (er .gt. 1000.0d0) then write (lout,150) nerr = nerr + 1 endif 270 tout = tout*10.0d0 275 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (jt .eq. 5) nfea = nfe - mband*nje write (lout,280) lenrw,leniw,nst,nfe,nfea,nje,ero 280 format(//' Final statistics for this run:'/ 1 ' rwork size =',i4,' iwork size =',i4/ 2 ' number of steps =',i5/ 3 ' number of f-s =',i5/ 4 ' (excluding J-s) =',i5/ 5 ' number of J-s =',i5/ 6 ' error overrun =',d10.2) 290 continue write (lout,300) nerr 300 format(///' Number of errors encountered =',i3) stop end subroutine f1 (neq, t, y, ydot) integer neq double precision t, y, ydot dimension y(neq), ydot(neq) ydot(1) = y(2) ydot(2) = 20.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1) return end subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd double precision t, y, pd dimension y(neq), pd(nrowpd,neq) pd(1,1) = 0.0d0 pd(1,2) = 1.0d0 pd(2,1) = -40.0d0*y(1)*y(2) - 1.0d0 pd(2,2) = 20.0d0*(1.0d0 - y(1)*y(1)) return end subroutine f2 (neq, t, y, ydot) integer neq, i, j, k, ng double precision t, y, ydot, alph1, alph2, d dimension y(neq), ydot(neq) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ do 10 j = 1,ng do 10 i = 1,ng k = i + (j - 1)*ng d = -2.0d0*y(k) if (i .ne. 1) d = d + y(k-1)*alph1 if (j .ne. 1) d = d + y(k-ng)*alph2 10 ydot(k) = d return end subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng double precision t, y, pd, alph1, alph2 dimension y(neq), pd(nrowpd,neq) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ mband = ml + mu + 1 mu1 = mu + 1 mu2 = mu + 2 do 10 j = 1,neq pd(mu1,j) = -2.0d0 pd(mu2,j) = alph1 10 pd(mband,j) = alph2 do 20 j = ng,neq,ng 20 pd(mu2,j) = 0.0d0 return end subroutine edit2 (y, t, erm) integer i, j, k, ng double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt dimension y(*) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ erm = 0.0d0 if (t .eq. 0.0d0) return ex = 0.0d0 if (t .le. 30.0d0) ex = exp(-2.0d0*t) a2 = 1.0d0 do 60 j = 1,ng a1 = 1.0d0 do 50 i = 1,ng k = i + (j - 1)*ng yt = t**(i+j-2)*ex*a1*a2 er = abs(y(k)-yt) erm = max(erm,er) a1 = a1*alph1/i 50 continue a2 = a2*alph2/j 60 continue return end ................................................................................ Demonstration program for DLSODA package Problem 1: Van der Pol oscillator: xdotdot - 20*(1 - x**2)*xdot + x = 0, x(0) = 2, xdot(0) = 0 neq = 2 itol = 1 rtol = 0.0E+00 atol = 0.1E-07 Solution with jt = 1 t x xdot meth nq h tsw 0.84609E+01 0.16731E+01 -0.464E-01 2 4 0.209E+00 0.311E+00 0.16922E+02 -0.11574E-03 -0.141E+02 1 7 0.206E-02 0.158E+02 0.25592E+02 -0.16828E+01 0.459E-01 2 4 0.240E+00 0.174E+02 0.34263E+02 0.21448E-03 0.141E+02 1 8 0.293E-02 0.332E+02 Final statistics for this run: rwork size = 52 iwork size = 22 number of steps = 695 number of f-s = 1305 (excluding J-s) = 1305 number of J-s = 30 max. error at root = 0.21E-03 Solution with jt = 2 t x xdot meth nq h tsw 0.84609E+01 0.16731E+01 -0.464E-01 2 4 0.209E+00 0.311E+00 0.16922E+02 -0.11574E-03 -0.141E+02 1 7 0.206E-02 0.158E+02 0.25592E+02 -0.16828E+01 0.459E-01 2 4 0.240E+00 0.174E+02 0.34263E+02 0.21448E-03 0.141E+02 1 8 0.293E-02 0.332E+02 Final statistics for this run: rwork size = 52 iwork size = 22 number of steps = 695 number of f-s = 1365 (excluding J-s) = 1305 number of J-s = 30 max. error at root = 0.21E-03 -------------------------------------------------------------------------------- Problem 2: ydot = A * y , where A is a banded lower triangular matrix derived from 2-D advection PDE neq = 25 ml = 5 mu = 0 itol = 1 rtol = 0.0E+00 atol = 0.1E-05 Solution with jt = 4 t max.err. meth nq h tsw 0.10000E-01 0.476E-06 1 2 0.714E-02 0.000E+00 0.10000E+00 0.988E-06 1 4 0.343E-01 0.000E+00 0.10000E+01 0.431E-06 1 5 0.724E-01 0.000E+00 0.10000E+02 0.558E-07 1 3 0.323E+00 0.000E+00 0.10000E+03 0.127E-11 2 1 0.239E+03 0.170E+02 Final statistics for this run: rwork size = 522 iwork size = 45 number of steps = 105 number of f-s = 207 (excluding J-s) = 207 number of J-s = 3 error overrun = 0.99E+00 Solution with jt = 5 t max.err. meth nq h tsw 0.10000E-01 0.476E-06 1 2 0.714E-02 0.000E+00 0.10000E+00 0.988E-06 1 4 0.343E-01 0.000E+00 0.10000E+01 0.431E-06 1 5 0.724E-01 0.000E+00 0.10000E+02 0.558E-07 1 3 0.323E+00 0.000E+00 0.10000E+03 0.127E-11 2 1 0.239E+03 0.170E+02 Final statistics for this run: rwork size = 522 iwork size = 45 number of steps = 105 number of f-s = 225 (excluding J-s) = 207 number of J-s = 3 error overrun = 0.99E+00 Number of errors encountered = 0 ==========Source and Output for LSODAR Demonstration Program==================== c----------------------------------------------------------------------- c Demonstration program for the DLSODAR package. c This is the version of 14 June 2001. c c This version is in double precision. c c The DLSODAR package is used to solve two simple problems, c one nonstiff and one intermittently stiff. c If the errors are too large, or other difficulty occurs, c a warning message is printed. All output is on unit lout = 6. c----------------------------------------------------------------------- external f1, gr1, f2, jac2, gr2 integer iopt, iout, istate, itask, itol, iwork, jroot, jt, 1 kroot, leniw, lenrw, liw, lrw, lout, neq, nerr, ng, 2 nfe, nfea, nge, nje, nst double precision atol, er, ero, errt, rtol, rwork, 1 t, tout, tzero, y, yt dimension y(2), atol(2), rwork(57), iwork(22), jroot(2) data lout/6/ c nerr = 0 c----------------------------------------------------------------------- c First problem. c The initial value problem is: c dy/dt = ((2*log(y) + 8)/t - 5)*y, y(1) = 1, 1 .le. t .le. 6 c The solution is y(t) = exp(-t**2 + 5*t - 4) c The two root functions are: c g1 = ((2*log(y)+8)/t - 5)*y (= dy/dt) (with root at t = 2.5), c g2 = log(y) - 2.2491 (with roots at t = 2.47 and 2.53) c----------------------------------------------------------------------- c Set all input parameters and print heading. neq = 1 y(1) = 1.0d0 t = 1.0d0 tout = 2.0d0 itol = 1 rtol = 1.0d-6 atol(1) = 1.0d-6 itask = 1 istate = 1 iopt = 0 lrw = 44 liw = 21 jt = 2 ng = 2 write (lout,110) itol,rtol,atol(1),jt 110 format(/' Demonstration program for DLSODAR package'//// 1 ' First problem'/// 2 ' Problem is dy/dt = ((2*log(y)+8)/t - 5)*y, y(1) = 1'// 3 ' Solution is y(t) = exp(-t**2 + 5*t - 4)'// 4 ' Root functions are:'/ 5 10x,' g1 = dy/dt (root at t = 2.5)'/ 6 10x,' g2 = log(y) - 2.2491 (roots at t = 2.47 and t = 2.53)'// 7 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1// 8 ' jt =',i3///) c c Call DLSODAR in loop over tout values 2,3,4,5,6. ero = 0.0d0 do 180 iout = 1,5 120 continue call dlsodar(f1,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jdum,jt,gr1,ng,jroot) c c Print y and error in y, and print warning if error too large. yt = exp(-t*t + 5.0d0*t - 4.0d0) er = y(1) - yt write (lout,130) t,y(1),er 130 format(' At t =',d15.7,5x,'y =',d15.7,5x,'error =',d12.4) if (istate .lt. 0) go to 185 er = abs(er)/(rtol*abs(y(1)) + atol(1)) ero = max(ero,er) if (er .gt. 1000.0d0) then write (lout,140) 140 format(//' Warning: error exceeds 1000 * tolerance'//) nerr = nerr + 1 endif if (istate .ne. 3) go to 175 c c If a root was found, write results and check root location. c Then reset istate to 2 and return to DLSODAR call. write (lout,150) t,jroot(1),jroot(2) 150 format(/' Root found at t =',d15.7,5x,'jroot =',2i5) if (jroot(1) .eq. 1) errt = t - 2.5d0 if (jroot(2) .eq. 1 .and. t .le. 2.5d0) errt = t - 2.47d0 if (jroot(2) .eq. 1 .and. t .gt. 2.5d0) errt = t - 2.53d0 write (lout,160) errt 160 format(' Error in t location of root is',d12.4/) if (abs(errt) .gt. 1.0d-3) then write (lout,170) 170 format(//' Warning: root error exceeds 1.0d-3'//) nerr = nerr + 1 endif istate = 2 go to 120 c c If no root found, increment tout and loop back. 175 tout = tout + 1.0d0 180 continue c c Problem complete. Print final statistics. 185 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) nge = iwork(10) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (jt .eq. 2) nfea = nfe - neq*nje write (lout,190) lenrw,leniw,nst,nfe,nfea,nje,nge,ero 190 format(//' Final statistics for this run:'/ 1 ' rwork size =',i4,' iwork size =',i4/ 2 ' number of steps =',i5/ 3 ' number of f-s =',i5/ 4 ' (excluding j-s) =',i5/ 5 ' number of j-s =',i5/ 6 ' number of g-s =',i5/ 7 ' error overrun =',d10.2) c c----------------------------------------------------------------------- c Second problem (Van der Pol oscillator). c The initial value problem is (after reduction of 2nd order ODE): c dy1/dt = y2, dy2/dt = 100*(1 - y1**2)*y2 - y1, c y1(0) = 2, y2(0) = 0, 0 .le. t .le. 200 c The root function is g = y1. c An analytic solution is not known, but the zeros of y1 are known c to 15 figures for purposes of checking the accuracy. c----------------------------------------------------------------------- c Set tolerance parameters and print heading. itol = 2 rtol = 1.0d-6 atol(1) = 1.0d-6 atol(2) = 1.0d-4 write (lout,200) itol,rtol,atol(1),atol(2) 200 format(////80('*')//' Second problem (Van der Pol oscillator)'// 1 ' Problem is dy1/dt = y2, dy2/dt = 100*(1-y1**2)*y2 - y1'/ 2 ' y1(0) = 2, y2(0) = 0'// 3 ' Root function is g = y1'// 4 ' itol =',i3,' rtol =',d10.1,' atol =',2d10.1) c c Loop over jt = 1, 2. Set remaining parameters and print jt. do 290 jt = 1,2 neq = 2 y(1) = 2.0d0 y(2) = 0.0d0 t = 0.0d0 tout = 20.0d0 itask = 1 istate = 1 iopt = 0 lrw = 57 liw = 22 ng = 1 write (lout,210) jt 210 format(///' Solution with jt =',i2//) c c Call DLSODAR in loop over tout values 20,40,...,200. do 270 iout = 1,10 220 continue call dlsodar(f2,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac2,jt,gr2,ng,jroot) c c Print y1 and y2. write (lout,230) t,y(1),y(2) 230 format(' At t =',d15.7,5x,'y1 =',d15.7,5x,'y2 =',d15.7) if (istate .lt. 0) go to 275 if (istate .ne. 3) go to 265 c c If a root was found, write results and check root location. c Then reset istate to 2 and return to DLSODAR call. write (lout,240) t 240 format(/' Root found at t =',d15.7) kroot = int(t/81.2d0 + 0.5d0) tzero = 81.17237787055d0 + (kroot-1)*81.41853556212d0 errt = t - tzero write (lout,250) errt 250 format(' Error in t location of root is',d12.4//) if (abs(errt) .gt. 1.0d-1) then write (lout,260) 260 format(//' Warning: root error exceeds 1.0d-1'//) nerr = nerr + 1 endif istate = 2 go to 220 c c If no root found, increment tout and loop back. 265 tout = tout + 20.0d0 270 continue c c Problem complete. Print final statistics. 275 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) nge = iwork(10) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (jt .eq. 2) nfea = nfe - neq*nje write (lout,280) lenrw,leniw,nst,nfe,nfea,nje,nge 280 format(//' Final statistics for this run:'/ 1 ' rwork size =',i4,' iwork size =',i4/ 2 ' number of steps =',i5/ 3 ' number of f-s =',i5/ 4 ' (excluding j-s) =',i5/ 5 ' number of j-s =',i5/ 6 ' number of g-s =',i5) 290 continue c write (lout,300) nerr 300 format(///' Total number of errors encountered =',i3) stop end subroutine f1 (neq, t, y, ydot) integer neq double precision t, y, ydot dimension y(1), ydot(1) ydot(1) = ((2.0d0*log(y(1)) + 8.0d0)/t - 5.0d0)*y(1) return end subroutine gr1 (neq, t, y, ng, groot) integer neq, ng double precision t, y, groot dimension y(1), groot(2) groot(1) = ((2.0d0*log(y(1)) + 8.0d0)/t - 5.0d0)*y(1) groot(2) = log(y(1)) - 2.2491d0 return end subroutine f2 (neq, t, y, ydot) integer neq double precision t, y, ydot dimension y(2), ydot(2) ydot(1) = y(2) ydot(2) = 100.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1) return end subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd double precision t, y, pd dimension y(2), pd(nrowpd,2) pd(1,1) = 0.0d0 pd(1,2) = 1.0d0 pd(2,1) = -200.0d0*y(1)*y(2) - 1.0d0 pd(2,2) = 100.0d0*(1.0d0 - y(1)*y(1)) return end subroutine gr2 (neq, t, y, ng, groot) integer neq, ng double precision t, y, groot dimension y(2), groot(1) groot(1) = y(1) return end ................................................................................ Demonstration program for DLSODAR package First problem Problem is dy/dt = ((2*log(y)+8)/t - 5)*y, y(1) = 1 Solution is y(t) = exp(-t**2 + 5*t - 4) Root functions are: g1 = dy/dt (root at t = 2.5) g2 = log(y) - 2.2491 (roots at t = 2.47 and t = 2.53) itol = 1 rtol = 0.1E-05 atol = 0.1E-05 jt = 2 At t = 0.2000000E+01 y = 0.7389071E+01 error = 0.1534E-04 At t = 0.2469971E+01 y = 0.9479201E+01 error = 0.1631E-04 Root found at t = 0.2469971E+01 jroot = 0 1 Error in t location of root is -0.2865E-04 At t = 0.2500001E+01 y = 0.9487752E+01 error = 0.1613E-04 Root found at t = 0.2500001E+01 jroot = 1 0 Error in t location of root is 0.6800E-06 At t = 0.2530028E+01 y = 0.9479201E+01 error = 0.1601E-04 Root found at t = 0.2530028E+01 jroot = 0 1 Error in t location of root is 0.2813E-04 At t = 0.3000000E+01 y = 0.7389083E+01 error = 0.2667E-04 At t = 0.4000000E+01 y = 0.1000007E+01 error = 0.6703E-05 At t = 0.5000000E+01 y = 0.1831582E-01 error = 0.1814E-06 At t = 0.6000000E+01 y = 0.4536004E-04 error = -0.3989E-07 Final statistics for this run: rwork size = 42 iwork size = 21 number of steps = 71 number of f-s = 147 (excluding j-s) = 147 number of j-s = 0 number of g-s = 114 error overrun = 0.34E+01 ******************************************************************************** Second problem (Van der Pol oscillator) Problem is dy1/dt = y2, dy2/dt = 100*(1-y1**2)*y2 - y1 y1(0) = 2, y2(0) = 0 Root function is g = y1 itol = 2 rtol = 0.1E-05 atol = 0.1E-05 0.1E-03 Solution with jt = 1 At t = 0.2000000E+02 y1 = 0.1858228E+01 y2 = -0.7575094E-02 At t = 0.4000000E+02 y1 = 0.1693230E+01 y2 = -0.9068584E-02 At t = 0.6000000E+02 y1 = 0.1484608E+01 y2 = -0.1232742E-01 At t = 0.8000000E+02 y1 = 0.1086291E+01 y2 = -0.5840716E-01 At t = 0.8116520E+02 y1 = -0.1308482E-12 y2 = -0.6713980E+02 Root found at t = 0.8116520E+02 Error in t location of root is -0.7180E-02 At t = 0.1000000E+03 y1 = -0.1868862E+01 y2 = 0.7497304E-02 At t = 0.1200000E+03 y1 = -0.1705927E+01 y2 = 0.8930077E-02 At t = 0.1400000E+03 y1 = -0.1501740E+01 y2 = 0.1196163E-01 At t = 0.1600000E+03 y1 = -0.1148800E+01 y2 = 0.3568399E-01 At t = 0.1625761E+03 y1 = 0.1153602E-11 y2 = 0.6713972E+02 Root found at t = 0.1625761E+03 Error in t location of root is -0.1485E-01 At t = 0.1800000E+03 y1 = 0.1879384E+01 y2 = -0.7422067E-02 At t = 0.2000000E+03 y1 = 0.1718431E+01 y2 = -0.8798201E-02 Final statistics for this run: rwork size = 55 iwork size = 22 number of steps = 478 number of f-s = 931 (excluding j-s) = 931 number of j-s = 42 number of g-s = 513 Solution with jt = 2 At t = 0.2000000E+02 y1 = 0.1858228E+01 y2 = -0.7575094E-02 At t = 0.4000000E+02 y1 = 0.1693230E+01 y2 = -0.9068584E-02 At t = 0.6000000E+02 y1 = 0.1484608E+01 y2 = -0.1232742E-01 At t = 0.8000000E+02 y1 = 0.1086291E+01 y2 = -0.5840716E-01 At t = 0.8116520E+02 y1 = -0.8500767E-12 y2 = -0.6713980E+02 Root found at t = 0.8116520E+02 Error in t location of root is -0.7180E-02 At t = 0.1000000E+03 y1 = -0.1868862E+01 y2 = 0.7497304E-02 At t = 0.1200000E+03 y1 = -0.1705927E+01 y2 = 0.8930077E-02 At t = 0.1400000E+03 y1 = -0.1501740E+01 y2 = 0.1196163E-01 At t = 0.1600000E+03 y1 = -0.1148800E+01 y2 = 0.3568399E-01 At t = 0.1625761E+03 y1 = 0.1057454E-11 y2 = 0.6713972E+02 Root found at t = 0.1625761E+03 Error in t location of root is -0.1485E-01 At t = 0.1800000E+03 y1 = 0.1879384E+01 y2 = -0.7422067E-02 At t = 0.2000000E+03 y1 = 0.1718431E+01 y2 = -0.8798201E-02 Final statistics for this run: rwork size = 55 iwork size = 22 number of steps = 478 number of f-s = 1015 (excluding j-s) = 931 number of j-s = 42 number of g-s = 516 Total number of errors encountered = 0 ==========Source and Output for LSODPK Demonstration Program==================== c----------------------------------------------------------------------- c Demonstration program for DLSODPK. c ODE system from ns-species interaction pde in 2 dimensions. c This is the version of 14 June 2001. c c This version is in double precision. c----------------------------------------------------------------------- c This program solves a stiff ODE system that arises from a system c of partial differential equations. The PDE system is a food web c population model, with predator-prey interaction and diffusion on c the unit square in two dimensions. The dependent variable vector is c c 1 2 ns c c = (c , c , ..., c ) c c and the PDEs are as follows: c c i i i c dc /dt = d(i)*(c + c ) + f (x,y,c) (i=1,...,ns) c xx yy i c c where c i ns j c f (x,y,c) = c *(b(i) + sum a(i,j)*c ) c i j=1 c c The number of species is ns = 2*np, with the first np being prey and c the last np being predators. The coefficients a(i,j), b(i), d(i) are: c c a(i,i) = -a (all i) c a(i,j) = -g (i .le. np, j .gt. np) c a(i,j) = e (i .gt. np, j .le. np) c b(i) = b*(1 + alpha*x*y) (i .le. np) c b(i) = -b*(1 + alpha*x*y) (i .gt. np) c d(i) = dprey (i .le. np) c d(i) = dpred (i .gt. np) c c The various scalar parameters are set in subroutine setpar. c c The boundary conditions are: normal derivative = 0. c A polynomial in x and y is used to set the initial conditions. c c The PDEs are discretized by central differencing on a mx by my mesh. c c The ODE system is solved by DLSODPK using method flag values c mf = 10, 21, 22, 23, 24, 29. The final time is tmax = 10, except c that for mf = 10 it is tmax = 1.0d-3 because the problem is stiff, c and for mf = 23 and 24 it is tmax = 2 because the lack of symmetry c in the problem makes these methods more costly. c c Two preconditioner matrices are used. One uses a fixed number of c Gauss-Seidel iterations based on the diffusion terms only. c The other preconditioner is a block-diagonal matrix based on c the partial derivatives of the interaction terms f only, using c block-grouping (computing only a subset of the ns by ns blocks). c For mf = 21 and 22, these two preconditioners are applied on c the left and right, respectively, and for mf = 23 and 24 the product c of the two is used as the one preconditioner matrix. c For mf = 29, the inverse of the product is applied. c c Two output files are written: one with the problem description and c and performance statistics on unit 6, and one with solution profiles c at selected output times (for mf = 22 only) on unit 8. c----------------------------------------------------------------------- c Note: In addition to the main program and 10 subroutines c given below, this program requires the LINPACK subroutines c DGEFA and DGESL, and the BLAS routine DAXPY. c----------------------------------------------------------------------- c Reference: c Peter N. Brown and Alan C. Hindmarsh, c Reduced Storage Matrix Methods in Stiff ODE Systems, c J. Appl. Math. & Comp., 31 (1989), pp. 40-91; c Also LLNL Report UCRL-95088, Rev. 1, June 1987. c----------------------------------------------------------------------- external fweb, jacbg, solsbg integer ns, mx, my, mxns, 1 mp, mq, mpsq, itmax, 2 meshx,meshy,ngx,ngy,ngrp,mxmp,jgx,jgy,jigx,jigy,jxr,jyr integer i, imod3, iopt, iout, istate, itask, itol, iwork, 1 jacflg, jpre, leniw, lenrw, liw, lrw, mf, 2 ncfl, ncfn, neq, nfe, nfldif, nfndif, nli, nlidif, nni, nnidif, 3 nout, npe, nps, nqu, nsdif, nst double precision aa, ee, gg, bb, dprey, dpred, 1 ax, ay, acoef, bcoef, dx, dy, alph, diff, cox, coy, 2 uround, srur double precision avdim, atol, cc, hu, rcfl, rcfn, dumach, 1 rtol, rwork, t, tout c c The problem Common blocks below allow for up to 20 species, c up to a 50x50 mesh, and up to a 20x20 group structure. common /pcom0/ aa, ee, gg, bb, dprey, dpred common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph, 1 diff(20), cox(20), coy(20), ns, mx, my, mxns common /pcom2/ uround, srur, mp, mq, mpsq, itmax common /pcom3/ meshx, meshy, ngx, ngy, ngrp, mxmp, 2 jgx(21), jgy(21), jigx(50), jigy(50), jxr(20), jyr(20) c c The dimension of cc below must be .ge. 2*neq, where neq = ns*mx*my. c The dimension lrw of rwork must be .ge. 17*neq + ns*ns*ngrp + 61, c and the dimension liw of iwork must be .ge. 35 + ns*ngrp. dimension cc(576), rwork(5213), iwork(67) data lrw/5213/, liw/67/ c open (unit=6, file='demout', status='new') open (unit=8, file='ccout', status='new') c ax = 1.0d0 ay = 1.0d0 c c Call setpar to set problem parameters. call setpar c c Set remaining problem parameters. neq = ns*mx*my mxns = mx*ns dx = ax/(mx-1) dy = ay/(my-1) do 10 i = 1,ns cox(i) = diff(i)/dx**2 10 coy(i) = diff(i)/dy**2 c c Write heading. write(6,20)ns, mx,my,neq 20 format(' Demonstration program for DLSODPK package'// 1 ' Food web problem with ns species, ns =',i4/ 2 ' Predator-prey interaction and diffusion on a 2-d square'// 3 ' Mesh dimensions (mx,my) =',2i4/ 4 ' Total system size is neq =',i7//) write(6,25) aa,ee,gg,bb,dprey,dpred,alph 25 format(' Matrix parameters: a =',d12.4,' e =',d12.4, 1 ' g =',d12.4/20x,' b =',d12.4// 2 ' Diffusion coefficients: dprey =',d12.4,' dpred =',d12.4/ 3 ' Rate parameter alpha =',d12.4//) c c Set remaining method parameters. jpre = 3 jacflg = 1 iwork(3) = jpre iwork(4) = jacflg iopt = 0 mp = ns mq = mx*my mpsq = ns*ns uround = dumach() srur = sqrt(uround) meshx = mx meshy = my mxmp = meshx*mp ngx = 2 ngy = 2 ngrp = ngx*ngy call gset (meshx, ngx, jgx, jigx, jxr) call gset (meshy, ngy, jgy, jigy, jyr) iwork(1) = mpsq*ngrp iwork(2) = mp*ngrp itmax = 5 itol = 1 rtol = 1.0d-5 atol = rtol itask = 1 write(6,30)ngrp,ngx,ngy,itmax,rtol,atol 30 format(' Preconditioning uses interaction-only block-diagonal', 1 ' matrix'/' with block-grouping, and Gauss-Seidel iterations'// 2 ' Number of diagonal block groups = ngrp =',i4, 3 ' (ngx by ngy, ngx =',i2,' ngy =',i2,' )'// 4 ' G-S preconditioner uses itmax iterations, itmax =',i3// 5 ' Tolerance parameters: rtol =',d10.2,' atol =',d10.2) c c c Loop over mf values 10, 21, ..., 24, 29. c do 90 mf = 10,29 if (mf .gt. 10 .and. mf .lt. 21) go to 90 if (mf .gt. 24 .and. mf .lt. 29) go to 90 write(6,40)mf 40 format(//80('-')//' Solution with mf =',i3// 1 ' t nstep nfe nni nli npe nq', 2 4x,'h avdim ncf rate lcf rate') c t = 0.0d0 tout = 1.0d-8 nout = 18 if (mf .eq. 10) nout = 6 if (mf .eq. 23 .or. mf .eq. 24) nout = 10 call cinit (cc) if (mf .eq. 22) call outweb (t, cc, ns, mx, my, 8) istate = 1 nli = 0 nni = 0 ncfn = 0 ncfl = 0 nst = 0 c c Loop over output times and call DLSODPK. c do 70 iout = 1,nout call dlsodpk (fweb, neq, cc, t, tout, itol, rtol, atol, itask, 1 istate, iopt, rwork, lrw, iwork, liw, jacbg, solsbg, mf) nsdif = iwork(11) - nst nst = iwork(11) nfe = iwork(12) npe = iwork(13) nnidif = iwork(19) - nni nni = iwork(19) nlidif = iwork(20) - nli nli = iwork(20) nfndif = iwork(22) - ncfn ncfn = iwork(22) nfldif = iwork(23) - ncfl ncfl = iwork(23) nqu = iwork(14) hu = rwork(11) avdim = 0.0d0 rcfn = 0.0d0 rcfl = 0.0d0 if (nnidif .gt. 0) avdim = real(nlidif)/real(nnidif) if (nsdif .gt. 0) rcfn = real(nfndif)/real(nsdif) if (nnidif .gt. 0) rcfl = real(nfldif)/real(nnidif) write(6,50)t,nst,nfe,nni,nli,npe,nqu,hu,avdim,rcfn,rcfl 50 format(d10.2,i5,i6,3i5,i4,2d11.2,d10.2,d12.2) imod3 = iout - 3*(iout/3) if (mf .eq. 22 .and. imod3 .eq. 0) call outweb (t,cc,ns,mx,my,8) if (istate .eq. 2) go to 65 write(6,60)t 60 format(//' final time reached =',d12.4//) go to 75 65 continue if (tout .gt. 0.9d0) tout = tout + 1.0d0 if (tout .lt. 0.9d0) tout = tout*10.0d0 70 continue c 75 continue nst = iwork(11) nfe = iwork(12) npe = iwork(13) lenrw = iwork(17) leniw = iwork(18) nni = iwork(19) nli = iwork(20) nps = iwork(21) if (nni .gt. 0) avdim = real(nli)/real(nni) ncfn = iwork(22) ncfl = iwork(23) write (6,80) lenrw,leniw,nst,nfe,npe,nps,nni,nli,avdim, 1 ncfn,ncfl 80 format(//' Final statistics for this run:'/ 1 ' rwork size =',i8,' iwork size =',i6/ 2 ' number of time steps =',i5/ 3 ' number of f evaluations =',i5/ 4 ' number of preconditioner evals. =',i5/ 4 ' number of preconditioner solves =',i5/ 5 ' number of nonlinear iterations =',i5/ 5 ' number of linear iterations =',i5/ 6 ' average subspace dimension =',f8.4/ 7 i5,' nonlinear conv. failures,',i5,' linear conv. failures') c 90 continue stop c------ end of main program for DLSODPK demonstration program ---------- end subroutine setpar c----------------------------------------------------------------------- c This routine sets the problem parameters. c It set ns, mx, my, and problem coefficients acoef, bcoef, diff, alph, c using parameters np, aa, ee, gg, bb, dprey, dpred. c----------------------------------------------------------------------- integer ns, mx, my, mxns integer i, j, np double precision aa, ee, gg, bb, dprey, dpred, 1 ax, ay, acoef, bcoef, dx, dy, alph, diff, cox, coy common /pcom0/ aa, ee, gg, bb, dprey, dpred common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph, 1 diff(20), cox(20), coy(20), ns, mx, my, mxns c np = 3 mx = 6 my = 6 aa = 1.0d0 ee = 1.0d4 gg = 0.5d-6 bb = 1.0d0 dprey = 1.0d0 dpred = 0.5d0 alph = 1.0d0 ns = 2*np do 70 j = 1,np do 60 i = 1,np acoef(np+i,j) = ee acoef(i,np+j) = -gg 60 continue acoef(j,j) = -aa acoef(np+j,np+j) = -aa bcoef(j) = bb bcoef(np+j) = -bb diff(j) = dprey diff(np+j) = dpred 70 continue c return c------------ end of subroutine setpar ------------------------------- end subroutine gset (m, ng, jg, jig, jr) c----------------------------------------------------------------------- c This routine sets arrays jg, jig, and jr describing c a uniform partition of (1,2,...,m) into ng groups. c----------------------------------------------------------------------- integer m, ng, jg, jig, jr dimension jg(*), jig(*), jr(*) integer ig, j, len1, mper, ngm1 c mper = m/ng do 10 ig = 1,ng 10 jg(ig) = 1 + (ig - 1)*mper jg(ng+1) = m + 1 c ngm1 = ng - 1 len1 = ngm1*mper do 20 j = 1,len1 20 jig(j) = 1 + (j-1)/mper len1 = len1 + 1 do 25 j = len1,m 25 jig(j) = ng c do 30 ig = 1,ngm1 30 jr(ig) = 0.5d0 + (ig - 0.5d0)*mper jr(ng) = 0.5d0*(1 + ngm1*mper + m) c return c------------ end of subroutine gset --------------------------------- end subroutine cinit (cc) c----------------------------------------------------------------------- c This routine computes and loads the vector of initial values. c----------------------------------------------------------------------- double precision cc dimension cc(*) integer ns, mx, my, mxns integer i, ici, ioff, iyoff, jx, jy double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy double precision argx, argy, x, y common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph, 1 diff(20), cox(20), coy(20), ns, mx, my, mxns c do 20 jy = 1,my y = (jy-1)*dy argy = 16.0d0*y*y*(ay-y)*(ay-y) iyoff = mxns*(jy-1) do 10 jx = 1,mx x = (jx-1)*dx argx = 16.0d0*x*x*(ax-x)*(ax-x) ioff = iyoff + ns*(jx-1) do 5 i = 1,ns ici = ioff + i cc(ici) = 10.0d0 + i*argx*argy 5 continue 10 continue 20 continue return c------------ end of subroutine cinit -------------------------------- end subroutine outweb (t, c, ns, mx, my, lun) c----------------------------------------------------------------------- c This routine prints the values of the individual species densities c at the current time t. The write statements use unit lun. c----------------------------------------------------------------------- integer ns, mx, my, lun double precision t, c dimension c(ns,mx,my) integer i, jx, jy c write(lun,10) t 10 format(/80('-')/30x,'At time t = ',d16.8/80('-') ) c do 40 i = 1,ns write(lun,20) i 20 format(' the species c(',i2,') values are:') do 30 jy = my,1,-1 write(lun,25) (c(i,jx,jy),jx=1,mx) 25 format(6(1x,g12.6)) 30 continue write(lun,35) 35 format(80('-'),/) 40 continue c return c------------ end of subroutine outweb ------------------------------- end subroutine fweb (neq, t, cc, cdot) c----------------------------------------------------------------------- c This routine computes the derivative of cc and returns it in cdot. c The interaction rates are computed by calls to webr, and these are c saved in cc(neq+1),...,cc(2*neq) for use in preconditioning. c----------------------------------------------------------------------- integer neq double precision t, cc, cdot dimension cc(neq), cdot(neq) integer ns, mx, my, mxns integer i, ic, ici, idxl, idxu, idyl, idyu, iyoff, jx, jy double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy double precision dcxli, dcxui, dcyli, dcyui, x, y common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph, 1 diff(20), cox(20), coy(20), ns, mx, my, mxns c do 100 jy = 1,my y = (jy-1)*dy iyoff = mxns*(jy-1) idyu = mxns if (jy .eq. my) idyu = -mxns idyl = mxns if (jy .eq. 1) idyl = -mxns do 90 jx = 1,mx x = (jx-1)*dx ic = iyoff + ns*(jx-1) + 1 c Get interaction rates at one point (x,y). call webr (x, y, t, cc(ic), cc(neq+ic)) idxu = ns if (jx .eq. mx) idxu = -ns idxl = ns if (jx .eq. 1) idxl = -ns do 80 i = 1,ns ici = ic + i - 1 c Do differencing in y. dcyli = cc(ici) - cc(ici-idyl) dcyui = cc(ici+idyu) - cc(ici) c Do differencing in x. dcxli = cc(ici) - cc(ici-idxl) dcxui = cc(ici+idxu) - cc(ici) c Collect terms and load cdot elements. cdot(ici) = coy(i)*(dcyui - dcyli) + cox(i)*(dcxui - dcxli) 1 + cc(neq+ici) 80 continue 90 continue 100 continue return c------------ end of subroutine fweb --------------------------------- end subroutine webr (x, y, t, c, rate) c----------------------------------------------------------------------- c This routine computes the interaction rates for the species c c(1),...,c(ns), at one spatial point and at time t. c----------------------------------------------------------------------- double precision x, y, t, c, rate dimension c(*), rate(*) integer ns, mx, my, mxns integer i double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy double precision fac common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph, 1 diff(20), cox(20), coy(20), ns, mx, my, mxns c do 10 i = 1,ns 10 rate(i) = 0.0d0 do 15 i = 1,ns call daxpy (ns, c(i), acoef(1,i), 1, rate, 1) 15 continue fac = 1.0d0 + alph*x*y do 20 i = 1,ns 20 rate(i) = c(i)*(bcoef(i)*fac + rate(i)) return c------------ end of subroutine webr --------------------------------- end subroutine jacbg (f, neq, t, cc, ccsv, rewt, f0, f1, hl0, 1 bd, ipbd, ier) c----------------------------------------------------------------------- c This routine generates part of the block-diagonal part of the c Jacobian, multiplies by -hl0, adds the identity matrix, c and calls DGEFA to do LU decomposition of each diagonal block. c The computation of the diagonal blocks uses the block and grouping c information in /pcom1/ and /pcom2/. One block per group is computed. c The Jacobian elements are generated by difference quotients c using calls to the routine fbg. c----------------------------------------------------------------------- c The two Common blocks below are used for internal communication. c The variables used are: c mp = size of blocks in block-diagonal preconditioning matrix. c mq = number of blocks in each direction (neq = mp*mq). c mpsq = mp*mp. c uround = unit roundoff, generated by a call uround = dumach(). c srur = sqrt(uround). c meshx = x mesh size c meshy = y mesh size (mesh is meshx by meshy) c ngx = no. groups in x direction in block-grouping scheme. c ngy = no. groups in y direction in block-grouping scheme. c ngrp = total number of groups = ngx*ngy. c mxmp = meshx*mp. c jgx = length ngx+1 array of group boundaries in x direction. c group igx has x indices jx = jgx(igx),...,jgx(igx+1)-1. c jigx = length meshx array of x group indices vs x node index. c x node index jx is in x group jigx(jx). c jxr = length ngx array of x indices representing the x groups. c the index for x group igx is jx = jxr(igx). c jgy, jigy, jyr = analogous arrays for grouping in y direction. c----------------------------------------------------------------------- external f integer neq, ipbd, ier double precision t, cc, ccsv, rewt, f0, f1, hl0, bd dimension cc(neq), ccsv(neq), rewt(neq), f0(neq), f1(neq), 1 bd(*), ipbd(*) integer mp, mq, mpsq, itmax, 2 meshx,meshy,ngx,ngy,ngrp,mxmp,jgx,jgy,jigx,jigy,jxr,jyr integer i, ibd, idiag, if0, if00, ig, igx, igy, iip, 1 j, jj, jx, jy, n double precision uround, srur double precision fac, r, r0, dvnorm c common /pcom2/ uround, srur, mp, mq, mpsq, itmax common /pcom3/ meshx, meshy, ngx, ngy, ngrp, mxmp, 1 jgx(21), jgy(21), jigx(50), jigy(50), jxr(20), jyr(20) c n = neq c c----------------------------------------------------------------------- c Make mp calls to fbg to approximate each diagonal block of Jacobian. c Here cc(neq+1),...,cc(2*neq) contains the base fb value. c r0 is a minimum increment factor for the difference quotient. c----------------------------------------------------------------------- 200 fac = dvnorm (n, f0, rewt) r0 = 1000.0d0*abs(hl0)*uround*n*fac if (r0 .eq. 0.0d0) r0 = 1.0d0 ibd = 0 do 240 igy = 1,ngy jy = jyr(igy) if00 = (jy - 1)*mxmp do 230 igx = 1,ngx jx = jxr(igx) if0 = if00 + (jx - 1)*mp do 220 j = 1,mp jj = if0 + j r = max(srur*abs(cc(jj)),r0/rewt(jj)) cc(jj) = cc(jj) + r fac = -hl0/r call fbg (neq, t, cc, jx, jy, f1) do 210 i = 1,mp 210 bd(ibd+i) = (f1(i) - cc(neq+if0+i))*fac cc(jj) = ccsv(jj) ibd = ibd + mp 220 continue 230 continue 240 continue c c Add identity matrix and do LU decompositions on blocks. -------------- 260 continue ibd = 1 iip = 1 do 280 ig = 1,ngrp idiag = ibd do 270 i = 1,mp bd(idiag) = bd(idiag) + 1.0d0 270 idiag = idiag + (mp + 1) call dgefa (bd(ibd), mp, mp, ipbd(iip), ier) if (ier .ne. 0) go to 290 ibd = ibd + mpsq iip = iip + mp 280 continue 290 return c------------ end of subroutine jacbg -------------------------------- end subroutine fbg (neq, t, cc, jx, jy, cdot) c----------------------------------------------------------------------- c This routine computes one block of the interaction terms of the c system, namely block (jx,jy), for use in preconditioning. c----------------------------------------------------------------------- integer neq, jx, jy double precision t, cc, cdot dimension cc(neq), cdot(neq) integer ns, mx, my, mxns integer iblok, ic double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy double precision x, y c common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph, 1 diff(20), cox(20), coy(20), ns, mx, my, mxns c iblok = jx + (jy-1)*mx y = (jy-1)*dy x = (jx-1)*dx ic = ns*(iblok-1) + 1 call webr (x, y, t, cc(ic), cdot) return c------------ end of subroutine fbg ---------------------------------- end subroutine solsbg (n, t, cc, f0, wk, hl0, bd, ipbd, v, lr, ier) c----------------------------------------------------------------------- c This routine applies one or two inverse preconditioner matrices c to the array v, using the interaction-only block-diagonal Jacobian c with block-grouping, and Gauss-Seidel applied to the diffusion terms. c When lr = 1 or 3, it calls gs for a Gauss-Seidel approximation c to ((I-hl0*Jd)-inverse)*v, and stores the result in v. c When lr = 2 or 3, it computes ((I-hl0*dg/dc)-inverse)*v, using LU c factors of the blocks in bd, and pivot information in ipbd. c In both cases, the array v is overwritten with the solution. c----------------------------------------------------------------------- integer n, ipbd, lr, ier double precision t, cc, f0, wk, hl0, bd, v dimension cc(n), f0(n), wk(n), bd(*), ipbd(*), v(n) integer mp, mq, mpsq, itmax, 2 meshx,meshy,ngx,ngy,ngrp,mxmp,jgx,jgy,jigx,jigy,jxr,jyr integer ibd, ig0, igm1, igx, igy, iip, iv, jx, jy double precision uround, srur c common /pcom2/ uround, srur, mp, mq, mpsq, itmax common /pcom3/ meshx, meshy, ngx, ngy, ngrp, mxmp, 1 jgx(21), jgy(21), jigx(50), jigy(50), jxr(20), jyr(20) c ier = 0 c if (lr.eq.0 .or. lr.eq.1 .or. lr.eq.3) call gs (n, hl0, v, wk) if (lr.eq.0 .or. lr.eq.2 .or. lr.eq.3) then iv = 1 do 20 jy = 1,meshy igy = jigy(jy) ig0 = (igy - 1)*ngx do 10 jx = 1,meshx igx = jigx(jx) igm1 = igx - 1 + ig0 ibd = 1 + igm1*mpsq iip = 1 + igm1*mp call dgesl (bd(ibd), mp, mp, ipbd(iip), v(iv), 0) iv = iv + mp 10 continue 20 continue endif c return c------------ end of subroutine solsbg ------------------------------- end subroutine gs (n, hl0, z, x) c----------------------------------------------------------------------- c This routine performs itmax Gauss-Seidel iterations c to compute an approximation to P-inverse*z, where P = I - hl0*Jd, and c Jd represents the diffusion contributions to the Jacobian. c z contains the answer on return. c The dimensions below assume ns .le. 20. c----------------------------------------------------------------------- integer n double precision hl0, z, x dimension z(n), x(n) integer ns, mx, my, mxns, 1 mp, mq, mpsq, itmax integer i, ic, ici, iter, iyoff, jx, jy double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy, 2 uround, srur double precision beta,beta2,cof1,elamda,gamma,gamma2 dimension beta(20), gamma(20), beta2(20), gamma2(20), cof1(20) common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph, 1 diff(20), cox(20), coy(20), ns, mx, my, mxns common /pcom2/ uround, srur, mp, mq, mpsq, itmax c c----------------------------------------------------------------------- c Write matrix as P = D - L - U. c Load local arrays beta, beta2, gamma, gamma2, and cof1. c----------------------------------------------------------------------- do 10 i = 1,ns elamda = 1.d0/(1.d0 + 2.d0*hl0*(cox(i) + coy(i))) beta(i) = hl0*cox(i)*elamda beta2(i) = 2.d0*beta(i) gamma(i) = hl0*coy(i)*elamda gamma2(i) = 2.d0*gamma(i) cof1(i) = elamda 10 continue c----------------------------------------------------------------------- c Begin iteration loop. c Load array x with (D-inverse)*z for first iteration. c----------------------------------------------------------------------- iter = 1 c do 50 jy = 1,my iyoff = mxns*(jy-1) do 40 jx = 1,mx ic = iyoff + ns*(jx-1) do 30 i = 1,ns ici = ic + i x(ici) = cof1(i)*z(ici) z(ici) = 0.d0 30 continue 40 continue 50 continue go to 160 c----------------------------------------------------------------------- c Calculate (D-inverse)*U*x. c----------------------------------------------------------------------- 70 continue iter = iter + 1 jy = 1 jx = 1 ic = ns*(jx-1) do 75 i = 1,ns ici = ic + i 75 x(ici) = beta2(i)*x(ici+ns) + gamma2(i)*x(ici+mxns) do 85 jx = 2,mx-1 ic = ns*(jx-1) do 80 i = 1,ns ici = ic + i 80 x(ici) = beta(i)*x(ici+ns) + gamma2(i)*x(ici+mxns) 85 continue jx = mx ic = ns*(jx-1) do 90 i = 1,ns ici = ic + i 90 x(ici) = gamma2(i)*x(ici+mxns) do 115 jy = 2,my-1 iyoff = mxns*(jy-1) jx = 1 ic = iyoff do 95 i = 1,ns ici = ic + i 95 x(ici) = beta2(i)*x(ici+ns) + gamma(i)*x(ici+mxns) do 105 jx = 2,mx-1 ic = iyoff + ns*(jx-1) do 100 i = 1,ns ici = ic + i 100 x(ici) = beta(i)*x(ici+ns) + gamma(i)*x(ici+mxns) 105 continue jx = mx ic = iyoff + ns*(jx-1) do 110 i = 1,ns ici = ic + i 110 x(ici) = gamma(i)*x(ici+mxns) 115 continue jy = my iyoff = mxns*(jy-1) jx = 1 ic = iyoff do 120 i = 1,ns ici = ic + i 120 x(ici) = beta2(i)*x(ici+ns) do 130 jx = 2,mx-1 ic = iyoff + ns*(jx-1) do 125 i = 1,ns ici = ic + i 125 x(ici) = beta(i)*x(ici+ns) 130 continue jx = mx ic = iyoff + ns*(jx-1) do 135 i = 1,ns ici = ic + i 135 x(ici) = 0.0d0 c----------------------------------------------------------------------- c Calculate (I - (D-inverse)*L)-inverse * x. c----------------------------------------------------------------------- 160 continue jy = 1 do 175 jx = 2,mx-1 ic = ns*(jx-1) do 170 i = 1,ns ici = ic + i 170 x(ici) = x(ici) + beta(i)*x(ici-ns) 175 continue jx = mx ic = ns*(jx-1) do 180 i = 1,ns ici = ic + i 180 x(ici) = x(ici) + beta2(i)*x(ici-ns) do 210 jy = 2,my-1 iyoff = mxns*(jy-1) jx = 1 ic = iyoff do 185 i = 1,ns ici = ic + i 185 x(ici) = x(ici) + gamma(i)*x(ici-mxns) do 200 jx = 2,mx-1 ic = iyoff + ns*(jx-1) do 195 i = 1,ns ici = ic + i x(ici) = (x(ici) + beta(i)*x(ici-ns)) 1 + gamma(i)*x(ici-mxns) 195 continue 200 continue jx = mx ic = iyoff + ns*(jx-1) do 205 i = 1,ns ici = ic + i x(ici) = (x(ici) + beta2(i)*x(ici-ns)) 1 + gamma(i)*x(ici-mxns) 205 continue 210 continue jy = my iyoff = mxns*(jy-1) jx = 1 ic = iyoff do 215 i = 1,ns ici = ic + i 215 x(ici) = x(ici) + gamma2(i)*x(ici-mxns) do 225 jx = 2,mx-1 ic = iyoff + ns*(jx-1) do 220 i = 1,ns ici = ic + i x(ici) = (x(ici) + beta(i)*x(ici-ns)) 1 + gamma2(i)*x(ici-mxns) 220 continue 225 continue jx = mx ic = iyoff + ns*(jx-1) do 230 i = 1,ns ici = ic + i x(ici) = (x(ici) + beta2(i)*x(ici-ns)) 1 + gamma2(i)*x(ici-mxns) 230 continue c----------------------------------------------------------------------- c Add increment x to z. c----------------------------------------------------------------------- do 300 i = 1,n 300 z(i) = z(i) + x(i) c if (iter .lt. itmax) go to 70 return c------------ end of subroutine gs ----------------------------------- end ................................................................................ Demonstration program for DLSODPK package Food web problem with ns species, ns = 6 Predator-prey interaction and diffusion on a 2-d square Mesh dimensions (mx,my) = 6 6 Total system size is neq = 216 Matrix parameters: a = 0.1000E+01 e = 0.1000E+05 g = 0.5000E-06 b = 0.1000E+01 Diffusion coefficients: dprey = 0.1000E+01 dpred = 0.5000E+00 Rate parameter alpha = 0.1000E+01 Preconditioning uses interaction-only block-diagonal matrix with block-grouping, and Gauss-Seidel iterations Number of diagonal block groups = ngrp = 4 (ngx by ngy, ngx = 2 ngy = 2 ) G-S preconditioner uses itmax iterations, itmax = 5 Tolerance parameters: rtol = 0.10E-04 atol = 0.10E-04 -------------------------------------------------------------------------------- Solution with mf = 10 t nstep nfe nni nli npe nq h avdim ncf rate lcf rate 0.10E-07 3 4 0 0 0 2 0.10E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-06 3 4 0 0 0 2 0.10E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-05 7 13 0 0 0 3 0.35E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-04 22 41 0 0 0 6 0.10E-05 0.00E+00 0.00E+00 0.00E+00 0.10E-03 81 155 0 0 0 2 0.24E-05 0.00E+00 0.34E-01 0.00E+00 0.10E-02 397 899 0 0 0 2 0.43E-05 0.00E+00 0.22E+00 0.00E+00 Final statistics for this run: rwork size = 3476 iwork size = 30 number of time steps = 397 number of f evaluations = 899 number of preconditioner evals. = 0 number of preconditioner solves = 0 number of nonlinear iterations = 0 number of linear iterations = 0 average subspace dimension = 0.0000 73 nonlinear conv. failures, 0 linear conv. failures -------------------------------------------------------------------------------- Solution with mf = 21 t nstep nfe nni nli npe nq h avdim ncf rate lcf rate 0.10E-07 3 5 3 1 2 2 0.65E-07 0.33E+00 0.00E+00 0.00E+00 0.10E-06 4 8 5 2 2 2 0.65E-07 0.50E+00 0.00E+00 0.00E+00 0.10E-05 10 22 13 8 4 3 0.37E-06 0.75E+00 0.00E+00 0.00E+00 0.10E-04 33 78 43 34 6 5 0.51E-06 0.87E+00 0.00E+00 0.00E+00 0.10E-03 112 270 134 135 14 5 0.60E-05 0.11E+01 0.00E+00 0.00E+00 0.10E-02 131 321 156 164 18 2 0.35E-03 0.13E+01 0.00E+00 0.00E+00 0.10E-01 138 343 167 175 20 3 0.18E-02 0.10E+01 0.00E+00 0.00E+00 0.10E+00 162 402 194 207 23 4 0.71E-02 0.12E+01 0.00E+00 0.00E+00 0.10E+01 206 540 242 297 27 4 0.47E-01 0.19E+01 0.00E+00 0.00E+00 0.20E+01 221 608 259 348 29 4 0.96E-01 0.30E+01 0.00E+00 0.00E+00 0.30E+01 230 656 269 386 30 4 0.16E+00 0.38E+01 0.00E+00 0.00E+00 0.40E+01 236 694 276 417 31 4 0.21E+00 0.44E+01 0.00E+00 0.29E+00 0.50E+01 240 718 280 437 31 4 0.27E+00 0.50E+01 0.00E+00 0.75E+00 0.60E+01 244 742 284 457 31 4 0.27E+00 0.50E+01 0.00E+00 0.75E+00 0.70E+01 247 766 288 477 32 4 0.36E+00 0.50E+01 0.00E+00 0.10E+01 0.80E+01 250 790 292 497 33 4 0.48E+00 0.50E+01 0.00E+00 0.10E+01 0.90E+01 252 802 294 507 33 4 0.48E+00 0.50E+01 0.00E+00 0.10E+01 0.10E+02 254 814 296 517 33 4 0.48E+00 0.50E+01 0.00E+00 0.10E+01 Final statistics for this run: rwork size = 3861 iwork size = 59 number of time steps = 254 number of f evaluations = 814 number of preconditioner evals. = 33 number of preconditioner solves = 1554 number of nonlinear iterations = 296 number of linear iterations = 517 average subspace dimension = 1.7466 0 nonlinear conv. failures, 20 linear conv. failures -------------------------------------------------------------------------------- Solution with mf = 22 t nstep nfe nni nli npe nq h avdim ncf rate lcf rate 0.10E-07 3 5 3 1 2 2 0.65E-07 0.33E+00 0.00E+00 0.00E+00 0.10E-06 4 8 5 2 2 2 0.65E-07 0.50E+00 0.00E+00 0.00E+00 0.10E-05 10 22 13 8 4 3 0.37E-06 0.75E+00 0.00E+00 0.00E+00 0.10E-04 33 78 43 34 6 5 0.51E-06 0.87E+00 0.00E+00 0.00E+00 0.10E-03 112 270 134 135 14 5 0.60E-05 0.11E+01 0.00E+00 0.00E+00 0.10E-02 131 321 156 164 18 2 0.35E-03 0.13E+01 0.00E+00 0.00E+00 0.10E-01 138 343 167 175 20 3 0.18E-02 0.10E+01 0.00E+00 0.00E+00 0.10E+00 162 402 194 207 23 4 0.71E-02 0.12E+01 0.00E+00 0.00E+00 0.10E+01 206 543 242 300 27 4 0.47E-01 0.19E+01 0.00E+00 0.00E+00 0.20E+01 221 612 259 352 29 4 0.95E-01 0.31E+01 0.00E+00 0.00E+00 0.30E+01 230 661 269 391 30 4 0.16E+00 0.39E+01 0.00E+00 0.10E+00 0.40E+01 236 695 275 419 30 4 0.20E+00 0.47E+01 0.00E+00 0.67E+00 0.50E+01 241 731 281 449 31 4 0.27E+00 0.50E+01 0.00E+00 0.10E+01 0.60E+01 244 749 284 464 31 4 0.27E+00 0.50E+01 0.00E+00 0.10E+01 0.70E+01 247 773 288 484 32 4 0.36E+00 0.50E+01 0.00E+00 0.10E+01 0.80E+01 250 797 292 504 33 3 0.51E+00 0.50E+01 0.00E+00 0.10E+01 0.90E+01 252 809 294 514 33 3 0.51E+00 0.50E+01 0.00E+00 0.10E+01 0.10E+02 254 827 297 529 34 2 0.82E+00 0.50E+01 0.00E+00 0.10E+01 Final statistics for this run: rwork size = 3877 iwork size = 54 number of time steps = 254 number of f evaluations = 827 number of preconditioner evals. = 34 number of preconditioner solves = 1580 number of nonlinear iterations = 297 number of linear iterations = 529 average subspace dimension = 1.7811 0 nonlinear conv. failures, 27 linear conv. failures -------------------------------------------------------------------------------- Solution with mf = 23 t nstep nfe nni nli npe nq h avdim ncf rate lcf rate 0.10E-07 3 5 3 1 2 2 0.65E-07 0.33E+00 0.00E+00 0.00E+00 0.10E-06 4 8 5 2 2 2 0.65E-07 0.50E+00 0.00E+00 0.00E+00 0.10E-05 10 22 13 8 4 3 0.37E-06 0.75E+00 0.00E+00 0.00E+00 0.10E-04 33 78 43 34 6 5 0.51E-06 0.87E+00 0.00E+00 0.00E+00 0.10E-03 112 272 134 137 14 5 0.60E-05 0.11E+01 0.00E+00 0.00E+00 0.10E-02 131 324 156 167 18 2 0.35E-03 0.14E+01 0.00E+00 0.00E+00 0.10E-01 139 360 166 193 20 3 0.15E-02 0.26E+01 0.00E+00 0.40E+00 0.10E+00 172 503 209 293 28 4 0.53E-02 0.23E+01 0.91E-01 0.26E+00 0.10E+01 227 845 278 566 42 4 0.11E-01 0.40E+01 0.91E-01 0.52E+00 0.20E+01 280 1287 360 926 70 2 0.17E-01 0.44E+01 0.25E+00 0.76E+00 Final statistics for this run: rwork size = 3404 iwork size = 54 number of time steps = 280 number of f evaluations = 1287 number of preconditioner evals. = 70 number of preconditioner solves = 926 number of nonlinear iterations = 360 number of linear iterations = 926 average subspace dimension = 2.5722 21 nonlinear conv. failures, 113 linear conv. failures -------------------------------------------------------------------------------- Solution with mf = 24 t nstep nfe nni nli npe nq h avdim ncf rate lcf rate 0.10E-07 3 5 3 1 2 2 0.65E-07 0.33E+00 0.00E+00 0.00E+00 0.10E-06 4 8 5 2 2 2 0.65E-07 0.50E+00 0.00E+00 0.00E+00 0.10E-05 10 22 13 8 4 3 0.37E-06 0.75E+00 0.00E+00 0.00E+00 0.10E-04 33 78 43 34 6 5 0.51E-06 0.87E+00 0.00E+00 0.00E+00 0.10E-03 112 270 134 135 14 5 0.60E-05 0.11E+01 0.00E+00 0.00E+00 0.10E-02 131 321 156 164 18 2 0.35E-03 0.13E+01 0.00E+00 0.00E+00 0.10E-01 139 356 167 188 20 3 0.16E-02 0.22E+01 0.00E+00 0.27E+00 0.10E+00 162 412 193 218 23 4 0.71E-02 0.12E+01 0.00E+00 0.00E+00 0.10E+01 223 780 274 505 38 4 0.23E-01 0.35E+01 0.82E-01 0.43E+00 0.20E+01 263 1085 335 749 59 3 0.17E-01 0.40E+01 0.23E+00 0.56E+00 Final statistics for this run: rwork size = 3404 iwork size = 54 number of time steps = 263 number of f evaluations = 1085 number of preconditioner evals. = 59 number of preconditioner solves = 749 number of nonlinear iterations = 335 number of linear iterations = 749 average subspace dimension = 2.2358 14 nonlinear conv. failures, 72 linear conv. failures -------------------------------------------------------------------------------- Solution with mf = 29 t nstep nfe nni nli npe nq h avdim ncf rate lcf rate 0.10E-07 3 4 3 0 2 2 0.65E-07 0.00E+00 0.00E+00 0.00E+00 0.10E-06 4 6 5 0 2 2 0.65E-07 0.00E+00 0.00E+00 0.00E+00 0.10E-05 10 14 13 0 4 3 0.37E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-04 32 42 41 0 6 5 0.56E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-03 114 135 134 0 13 5 0.47E-05 0.00E+00 0.00E+00 0.00E+00 0.10E-02 136 162 161 0 18 2 0.30E-03 0.00E+00 0.00E+00 0.00E+00 0.10E-01 144 173 172 0 20 3 0.16E-02 0.00E+00 0.00E+00 0.00E+00 0.10E+00 168 200 199 0 23 4 0.79E-02 0.00E+00 0.00E+00 0.00E+00 0.10E+01 245 353 352 0 39 2 0.17E-01 0.00E+00 0.13E-01 0.00E+00 0.20E+01 293 467 466 0 57 2 0.28E-01 0.00E+00 0.10E+00 0.00E+00 0.30E+01 330 566 565 0 76 2 0.37E-01 0.00E+00 0.22E+00 0.00E+00 0.40E+01 356 631 630 0 87 2 0.31E-01 0.00E+00 0.15E+00 0.00E+00 0.50E+01 384 697 696 0 98 1 0.72E-01 0.00E+00 0.14E+00 0.00E+00 0.60E+01 399 742 741 0 109 2 0.21E+00 0.00E+00 0.27E+00 0.00E+00 0.70E+01 411 783 782 0 117 1 0.20E+00 0.00E+00 0.33E+00 0.00E+00 0.80E+01 414 788 787 0 118 2 0.41E+00 0.00E+00 0.00E+00 0.00E+00 0.90E+01 416 791 790 0 118 2 0.41E+00 0.00E+00 0.00E+00 0.00E+00 0.10E+02 418 793 792 0 119 3 0.74E+00 0.00E+00 0.00E+00 0.00E+00 Final statistics for this run: rwork size = 2756 iwork size = 54 number of time steps = 418 number of f evaluations = 793 number of preconditioner evals. = 119 number of preconditioner solves = 777 number of nonlinear iterations = 792 number of linear iterations = 0 average subspace dimension = 0.0000 30 nonlinear conv. failures, 0 linear conv. failures ................................................................................ -------------------------------------------------------------------------------- At time t = 0.00000000E+00 -------------------------------------------------------------------------------- the species c( 1) values are: 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.1678 10.3775 10.3775 10.1678 10.0000 10.0000 10.3775 10.8493 10.8493 10.3775 10.0000 10.0000 10.3775 10.8493 10.8493 10.3775 10.0000 10.0000 10.1678 10.3775 10.3775 10.1678 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- the species c( 2) values are: 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.3355 10.7550 10.7550 10.3355 10.0000 10.0000 10.7550 11.6987 11.6987 10.7550 10.0000 10.0000 10.7550 11.6987 11.6987 10.7550 10.0000 10.0000 10.3355 10.7550 10.7550 10.3355 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- the species c( 3) values are: 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.5033 11.1325 11.1325 10.5033 10.0000 10.0000 11.1325 12.5480 12.5480 11.1325 10.0000 10.0000 11.1325 12.5480 12.5480 11.1325 10.0000 10.0000 10.5033 11.1325 11.1325 10.5033 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- the species c( 4) values are: 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.6711 11.5099 11.5099 10.6711 10.0000 10.0000 11.5099 13.3974 13.3974 11.5099 10.0000 10.0000 11.5099 13.3974 13.3974 11.5099 10.0000 10.0000 10.6711 11.5099 11.5099 10.6711 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- the species c( 5) values are: 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.8389 11.8874 11.8874 10.8389 10.0000 10.0000 11.8874 14.2467 14.2467 11.8874 10.0000 10.0000 11.8874 14.2467 14.2467 11.8874 10.0000 10.0000 10.8389 11.8874 11.8874 10.8389 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- the species c( 6) values are: 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 11.0066 12.2649 12.2649 11.0066 10.0000 10.0000 12.2649 15.0961 15.0961 12.2649 10.0000 10.0000 12.2649 15.0961 15.0961 12.2649 10.0000 10.0000 11.0066 12.2649 12.2649 11.0066 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.10000000E-05 -------------------------------------------------------------------------------- the species c( 1) values are: 9.99991 9.99992 9.99993 9.99993 9.99993 9.99992 9.99992 10.1677 10.3774 10.3774 10.1677 9.99993 9.99993 10.3774 10.8492 10.8492 10.3774 9.99993 9.99993 10.3774 10.8492 10.8492 10.3774 9.99993 9.99992 10.1677 10.3774 10.3774 10.1677 9.99992 9.99991 9.99992 9.99993 9.99993 9.99992 9.99991 -------------------------------------------------------------------------------- the species c( 2) values are: 9.99991 9.99993 9.99995 9.99995 9.99993 9.99992 9.99993 10.3355 10.7549 10.7549 10.3355 9.99993 9.99995 10.7549 11.6985 11.6985 10.7549 9.99995 9.99995 10.7549 11.6985 11.6985 10.7549 9.99995 9.99993 10.3355 10.7549 10.7549 10.3355 9.99993 9.99991 9.99993 9.99995 9.99995 9.99993 9.99991 -------------------------------------------------------------------------------- the species c( 3) values are: 9.99991 9.99994 9.99997 9.99997 9.99994 9.99992 9.99994 10.5032 11.1323 11.1323 10.5032 9.99994 9.99997 11.1323 12.5478 12.5478 11.1323 9.99997 9.99997 11.1323 12.5478 12.5478 11.1323 9.99997 9.99994 10.5032 11.1323 11.1323 10.5032 9.99994 9.99991 9.99994 9.99997 9.99997 9.99994 9.99991 -------------------------------------------------------------------------------- the species c( 4) values are: 13.4987 13.4987 13.4987 13.4987 13.4987 13.4987 13.4987 14.5503 15.8929 15.8929 14.5503 13.4987 13.4987 15.8929 19.0303 19.0303 15.8929 13.4987 13.4987 15.8929 19.0303 19.0303 15.8929 13.4987 13.4987 14.5503 15.8929 15.8929 14.5503 13.4987 13.4987 13.4987 13.4987 13.4987 13.4987 13.4987 -------------------------------------------------------------------------------- the species c( 5) values are: 13.4987 13.4987 13.4987 13.4987 13.4987 13.4987 13.4987 14.7791 16.4141 16.4141 14.7791 13.4987 13.4987 16.4141 20.2367 20.2367 16.4141 13.4987 13.4987 16.4141 20.2367 20.2367 16.4141 13.4987 13.4987 14.7791 16.4141 16.4141 14.7791 13.4987 13.4987 13.4987 13.4987 13.4987 13.4987 13.4987 -------------------------------------------------------------------------------- the species c( 6) values are: 13.4987 13.4987 13.4988 13.4987 13.4987 13.4987 13.4987 15.0078 16.9353 16.9353 15.0078 13.4987 13.4988 16.9353 21.4431 21.4431 16.9353 13.4987 13.4988 16.9353 21.4431 21.4431 16.9353 13.4988 13.4987 15.0078 16.9353 16.9353 15.0078 13.4987 13.4987 13.4987 13.4988 13.4988 13.4987 13.4987 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.10000000E-02 -------------------------------------------------------------------------------- the species c( 1) values are: 9.90702 9.91664 9.92836 9.93033 9.92253 9.91674 9.91472 10.0746 10.2769 10.2785 10.0795 9.92253 9.92446 10.2748 10.7181 10.7194 10.2785 9.93033 9.92445 10.2744 10.7173 10.7181 10.2769 9.92836 9.91469 10.0734 10.2744 10.2748 10.0746 9.91664 9.90697 9.91469 9.92445 9.92446 9.91472 9.90702 -------------------------------------------------------------------------------- the species c( 2) values are: 9.90741 9.92474 9.94623 9.94820 9.93064 9.91713 9.92282 10.2412 10.6440 10.6457 10.2461 9.93064 9.94232 10.6419 11.5267 11.5281 10.6457 9.94820 9.94231 10.6415 11.5258 11.5267 10.6440 9.94623 9.92279 10.2400 10.6415 10.6419 10.2412 9.92474 9.90736 9.92279 9.94231 9.94232 9.92282 9.90741 -------------------------------------------------------------------------------- the species c( 3) values are: 9.90780 9.93284 9.96409 9.96606 9.93874 9.91752 9.93092 10.4078 11.0109 11.0127 10.4127 9.93874 9.96018 11.0088 12.3339 12.3354 11.0127 9.96606 9.96017 11.0083 12.3329 12.3339 11.0109 9.96409 9.93089 10.4065 11.0083 11.0088 10.4078 9.93284 9.90776 9.93089 9.96017 9.96018 9.93092 9.90780 -------------------------------------------------------------------------------- the species c( 4) values are: 297231. 297749. 298393. 298451. 297925. 297520. 297692. 307244. 319327. 319378. 307390. 297925. 298276. 319264. 345799. 345840. 319378. 298451. 298276. 319252. 345771. 345799. 319327. 298393. 297691. 307208. 319252. 319264. 307244. 297749. 297229. 297691. 298276. 298276. 297692. 297231. -------------------------------------------------------------------------------- the species c( 5) values are: 297231. 297749. 298393. 298451. 297925. 297520. 297692. 307244. 319327. 319378. 307390. 297925. 298276. 319264. 345799. 345840. 319378. 298451. 298276. 319252. 345771. 345799. 319327. 298393. 297691. 307208. 319252. 319264. 307244. 297749. 297229. 297691. 298276. 298276. 297692. 297231. -------------------------------------------------------------------------------- the species c( 6) values are: 297231. 297749. 298393. 298451. 297925. 297520. 297692. 307244. 319327. 319378. 307390. 297925. 298276. 319264. 345799. 345840. 319378. 298451. 298276. 319252. 345771. 345799. 319327. 298393. 297691. 307208. 319252. 319264. 307244. 297749. 297229. 297691. 298276. 298276. 297692. 297231. -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.10000000E+01 -------------------------------------------------------------------------------- the species c( 1) values are: 1.58846 1.59918 1.62146 1.64759 1.67030 1.68143 1.58527 1.59498 1.61542 1.63946 1.66027 1.67030 1.57751 1.58542 1.60234 1.62229 1.63946 1.64759 1.56815 1.57406 1.58700 1.60234 1.61542 1.62146 1.56043 1.56457 1.57406 1.58542 1.59498 1.59918 1.55727 1.56043 1.56815 1.57751 1.58527 1.58846 -------------------------------------------------------------------------------- the species c( 2) values are: 1.59061 1.60135 1.62365 1.64981 1.67255 1.68369 1.58742 1.59714 1.61761 1.64167 1.66251 1.67255 1.57965 1.58757 1.60451 1.62449 1.64167 1.64981 1.57028 1.57620 1.58916 1.60451 1.61761 1.62365 1.56255 1.56670 1.57620 1.58757 1.59714 1.60135 1.55939 1.56255 1.57028 1.57965 1.58742 1.59061 -------------------------------------------------------------------------------- the species c( 3) values are: 1.59265 1.60340 1.62572 1.65191 1.67468 1.68583 1.58946 1.59918 1.61967 1.64377 1.66462 1.67468 1.58168 1.58960 1.60656 1.62656 1.64377 1.65191 1.57230 1.57823 1.59119 1.60656 1.61967 1.62572 1.56456 1.56872 1.57823 1.58960 1.59918 1.60340 1.56140 1.56456 1.57230 1.58168 1.58946 1.59265 -------------------------------------------------------------------------------- the species c( 4) values are: 47716.6 48038.6 48707.3 49491.7 50173.7 50507.6 47621.0 47912.3 48526.1 49247.9 49872.6 50173.7 47388.0 47625.3 48133.2 48732.4 49247.9 49491.7 47106.8 47284.3 47672.8 48133.2 48526.1 48707.3 46874.9 46999.4 47284.3 47625.3 47912.3 48038.6 46780.1 46874.9 47106.8 47388.0 47621.0 47716.6 -------------------------------------------------------------------------------- the species c( 5) values are: 47716.6 48038.6 48707.3 49491.7 50173.7 50507.6 47621.0 47912.3 48526.1 49247.9 49872.6 50173.7 47388.0 47625.3 48133.2 48732.4 49247.9 49491.7 47106.8 47284.3 47672.8 48133.2 48526.1 48707.3 46874.9 46999.4 47284.3 47625.3 47912.3 48038.6 46780.1 46874.9 47106.8 47388.0 47621.0 47716.6 -------------------------------------------------------------------------------- the species c( 6) values are: 47716.6 48038.6 48707.3 49491.7 50173.7 50507.6 47621.0 47912.3 48526.1 49247.9 49872.6 50173.7 47388.0 47625.3 48133.2 48732.4 49247.9 49491.7 47106.8 47284.3 47672.8 48133.2 48526.1 48707.3 46874.9 46999.4 47284.3 47625.3 47912.3 48038.6 46780.1 46874.9 47106.8 47388.0 47621.0 47716.6 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.40000000E+01 -------------------------------------------------------------------------------- the species c( 1) values are: 1.19534 1.20367 1.22108 1.24156 1.25934 1.26799 1.19279 1.20034 1.21635 1.23521 1.25152 1.25934 1.18656 1.19272 1.20601 1.22172 1.23521 1.24156