c----------------------------------------------------------------------- c demonstration program for the lsodes package. c this is the version of may 3, 1983. c c this version is in double precision. c c for computer systems requiring a program card, the following (with c the c in column 1 removed) may be used.. c c program lssdem(lssout,tape6=lssout) 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, miter1, moss, moss1, 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(1x/45h demonstration problem for the lsodes 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) = dfloat(i) ia(1) = 1 k = 1 do 60 m = 1,igrid do 50 l = 1,igrid j = l + (m - 1)*igrid if (m .eq. 1) go to 30 ja(k) = j - igrid k = k + 1 30 if (l .eq. 1) go to 35 ja(k) = j - 1 k = k + 1 35 ja(k) = j k = k + 1 if (l .eq. igrid) go to 40 ja(k) = j + 1 k = k + 1 40 ia(j+1) = k 50 continue 60 continue write (lout,80)neq,t,rtol,atol,(y(i),i=1,neq) 80 format(1x,6h neq =,i4,5x,4ht0 =,f4.1,5x,6hrtol =,e12.3,5x, 1 6hatol =,e12.3//1x,21h initial y vector = ,9f5.1////) c c loop over all relevant values of mf. do 193 moss1 = 1,3 moss = moss1 - 1 do 192 meth = 1,2 do 191 miter1 = 1,4 miter = miter1 - 1 if ( (miter.eq.0 .or. miter.eq.3) .and. moss.ne.0) go to 191 mf = 100*moss + 10*meth + miter if (mf .ne. 10) write (lout,100) 100 format(1h1) 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(///1x,4hmf =,i4,5x,29hinput work lengths lrw, liw =,2i6) do 125 i = 1,neq 125 y(i) = dfloat(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 lsodes (fdem, neq, y, t, tout, itol, rtol, atol, 1 itask, 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(//1x,7h at t =,f5.1,5x,5hnst =,i4,5x,4hhu =,e12.3,5x, 1 5hnqu =,i3,5x,12h max. err. =,e12.4/ 2 1x,15h y array = ,4e15.6/1x,5e15.6) if (istate .lt. 0) go to 175 erm = erm/atol ero = dmax1(ero,erm) if (erm .lt. 100.0d0) go to 160 write (lout,150) 150 format(//1x,40h warning.. error exceeds 100 * tolerance//) nerr = nerr + 1 160 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(//1x,32h final statistics for this run../ 1 1x,13h rwork size =,i4,15h iwork size =,i4/ 2 1x,18h number of steps =,i5/ 3 1x,18h number of f-s =,i5/ 4 1x,18h (excluding j-s) =,i5/ 5 1x,18h number of j-s =,i5/ 6 1x,16h error overrun =,e10.2) if (miter .eq. 1 .or. miter .eq. 2) 1 write (lout,185)nnz,ngp,nlu,nzl,nzu 185 format(1x,27h number of nonzeros in j = ,i5/ 1 1x,27h number of j index groups =,i5/ 2 1x,27h number of lu decomp-s =,i5/ 3 1x,34h nonzeros in strict lower factor =,i5/ 4 1x,34h 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,200) nerr 200 format(////1x,31h 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(1), ia(1), ja(1), pdj(1) 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(1),yex(9,3) data neq /9/ data yex /6.687279d-01, 9.901910d-01, 7.603061d-01, 1 8.077979d-01, 1.170226d+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 = dabs(y(i) - yex(i,iout)) 10 erm = dmax1(erm,er) return end subroutine ssout (neq, iwk, iwork, lout) integer neq, iwk, iwork, lout integer i, i1, i2, ipian, ipjan, nnz dimension iwk(1), iwork(1) ipian = iwork(23) ipjan = iwork(24) nnz = iwork(19) i1 = ipian i2 = i1 + neq write (lout,10)(iwk(i),i=i1,i2) 10 format(//1x,33h structure descriptor array ian =/(20i4)) i1 = ipjan i2 = i1 + nnz - 1 write (lout,20)(iwk(i),i=i1,i2) 20 format(/1x,33h structure descriptor array jan =/(20i4)) return end