share/odepack/fortran/opkdemo5.f

   1       program opkdemo5
   2 c-----------------------------------------------------------------------
   3 c Demonstration program for DLSODPK.
   4 c ODE system from ns-species interaction pde in 2 dimensions.
   5 c This is the version of 14 June 2001.
   6 c
   7 c This version is in double precision.
   8 c-----------------------------------------------------------------------
   9 c This program solves a stiff ODE system that arises from a system
  10 c of partial differential equations.  The PDE system is a food web
  11 c population model, with predator-prey interaction and diffusion on
  12 c the unit square in two dimensions.  The dependent variable vector is
  13 c
  14 c         1   2        ns
  15 c   c = (c , c , ..., c  )
  16 c
  17 c and the PDEs are as follows:
  18 c
  19 c     i               i      i
  20 c   dc /dt  =  d(i)*(c    + c   )  +  f (x,y,c)  (i=1,...,ns)
  21 c                     xx     yy        i
  22 c
  23 c where
  24 c                  i          ns         j
  25 c   f (x,y,c)  =  c *(b(i) + sum a(i,j)*c )
  26 c    i                       j=1
  27 c
  28 c The number of species is ns = 2*np, with the first np being prey and
  29 c the last np being predators.  The coefficients a(i,j), b(i), d(i) are:
  30 c
  31 c   a(i,i) = -a  (all i)
  32 c   a(i,j) = -g  (i .le. np, j .gt. np)
  33 c   a(i,j) =  e  (i .gt. np, j .le. np)
  34 c   b(i) =  b*(1 + alpha*x*y)  (i .le. np)
  35 c   b(i) = -b*(1 + alpha*x*y)  (i .gt. np)
  36 c   d(i) = dprey  (i .le. np)
  37 c   d(i) = dpred  (i .gt. np)
  38 c
  39 c The various scalar parameters are set in subroutine setpar.
  40 c
  41 c The boundary conditions are: normal derivative = 0.
  42 c A polynomial in x and y is used to set the initial conditions.
  43 c
  44 c The PDEs are discretized by central differencing on a mx by my mesh.
  45 c
  46 c The ODE system is solved by DLSODPK using method flag values
  47 c mf = 10, 21, 22, 23, 24, 29.  The final time is tmax = 10, except
  48 c that for mf = 10 it is tmax = 1.0d-3 because the problem is stiff,
  49 c and for mf = 23 and 24 it is tmax = 2 because the lack of symmetry
  50 c in the problem makes these methods more costly.
  51 c
  52 c Two preconditioner matrices are used.  One uses a fixed number of
  53 c Gauss-Seidel iterations based on the diffusion terms only.
  54 c The other preconditioner is a block-diagonal matrix based on
  55 c the partial derivatives of the interaction terms f only, using
  56 c block-grouping (computing only a subset of the ns by ns blocks).
  57 c For mf = 21 and 22, these two preconditioners are applied on
  58 c the left and right, respectively, and for mf = 23 and 24 the product
  59 c of the two is used as the one preconditioner matrix.
  60 c For mf = 29, the inverse of the product is applied.
  61 c
  62 c Two output files are written: one with the problem description and
  63 c and performance statistics on unit 6, and one with solution profiles
  64 c at selected output times (for mf = 22 only) on unit 8.
  65 c-----------------------------------------------------------------------
  66 c Note: In addition to the main program and 10 subroutines
  67 c given below, this program requires the LINPACK subroutines
  68 c DGEFA and DGESL, and the BLAS routine DAXPY.
  69 c-----------------------------------------------------------------------
  70 c Reference:
  71 c     Peter N. Brown and Alan C. Hindmarsh,
  72 c     Reduced Storage Matrix Methods in Stiff ODE Systems,
  73 c     J. Appl. Math. & Comp., 31 (1989), pp. 40-91;
  74 c     Also LLNL Report UCRL-95088, Rev. 1, June 1987.
  75 c-----------------------------------------------------------------------
  76       external fweb, jacbg, solsbg
  77       integer ns, mx, my, mxns,
  78      1        mp, mq, mpsq, itmax,
  79      2        meshx,meshy,ngx,ngy,ngrp,mxmp,jgx,jgy,jigx,jigy,jxr,jyr
  80       integer i, imod3, iopt, iout, istate, itask, itol, iwork,
  81      1   jacflg, jpre, leniw, lenrw, liw, lrw, mf,
  82      2   ncfl, ncfn, neq, nfe, nfldif, nfndif, nli, nlidif, nni, nnidif,
  83      3   nout, npe, nps, nqu, nsdif, nst
  84       double precision aa, ee, gg, bb, dprey, dpred,
  85      1     ax, ay, acoef, bcoef, dx, dy, alph, diff, cox, coy,
  86      2     uround, srur
  87       double precision avdim, atol, cc, hu, rcfl, rcfn, dumach,
  88      1   rtol, rwork, t, tout
  89       dimension neq(1), atol(1), rtol(1)
  90 c
  91 c The problem Common blocks below allow for up to 20 species,
  92 c up to a 50x50 mesh, and up to a 20x20 group structure.
  93       common /pcom0/ aa, ee, gg, bb, dprey, dpred
  94       common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph,
  95      1               diff(20), cox(20), coy(20), ns, mx, my, mxns
  96       common /pcom2/ uround, srur, mp, mq, mpsq, itmax
  97       common /pcom3/ meshx, meshy, ngx, ngy, ngrp, mxmp,
  98      2   jgx(21), jgy(21), jigx(50), jigy(50), jxr(20), jyr(20)
  99 c
 100 c The dimension of cc below must be .ge. 2*neq, where neq = ns*mx*my.
 101 c The dimension lrw of rwork must be .ge. 17*neq + ns*ns*ngrp + 61,
 102 c and the dimension liw of iwork must be .ge. 35 + ns*ngrp.
 103       dimension cc(576), rwork(5213), iwork(67)
 104       data lrw/5213/, liw/67/
 105 c
 106       open (unit=6, file='demout', status='new')
 107       open (unit=8, file='ccout', status='new')
 108 c
 109       ax = 1.0d0
 110       ay = 1.0d0
 111 c
 112 c Call setpar to set problem parameters.
 113       call setpar
 114 c
 115 c Set remaining problem parameters.
 116       neq(1) = ns*mx*my
 117       mxns = mx*ns
 118       dx = ax/(mx-1)
 119       dy = ay/(my-1)
 120       do 10 i = 1,ns
 121         cox(i) = diff(i)/dx**2
 122  10     coy(i) = diff(i)/dy**2
 123 c
 124 c Write heading.
 125       write(6,20)ns, mx,my,neq(1)
 126  20   format(' Demonstration program for DLSODPK package'//
 127      1   ' Food web problem with ns species, ns =',i4/
 128      2   ' Predator-prey interaction and diffusion on a 2-d square'//
 129      3   ' Mesh dimensions (mx,my) =',2i4/
 130      4   ' Total system size is neq =',i7//)
 131       write(6,25) aa,ee,gg,bb,dprey,dpred,alph
 132  25   format(' Matrix parameters:  a =',d12.4,'   e =',d12.4,
 133      1   '   g =',d12.4/20x,' b =',d12.4//
 134      2   ' Diffusion coefficients: dprey =',d12.4,'   dpred =',d12.4/
 135      3   ' Rate parameter alpha =',d12.4//)
 136 c
 137 c Set remaining method parameters.
 138       jpre = 3
 139       jacflg = 1
 140       iwork(3) = jpre
 141       iwork(4) = jacflg
 142       iopt = 0
 143       mp = ns
 144       mq = mx*my
 145       mpsq = ns*ns
 146       uround = dumach()
 147       srur = sqrt(uround)
 148       meshx = mx
 149       meshy = my
 150       mxmp = meshx*mp
 151       ngx = 2
 152       ngy = 2
 153       ngrp = ngx*ngy
 154       call gset (meshx, ngx, jgx, jigx, jxr)
 155       call gset (meshy, ngy, jgy, jigy, jyr)
 156       iwork(1) = mpsq*ngrp
 157       iwork(2) = mp*ngrp
 158       itmax = 5
 159       itol = 1
 160       rtol(1) = 1.0d-5
 161       atol(1) = rtol(1)
 162       itask = 1
 163       write(6,30)ngrp,ngx,ngy,itmax,rtol(1),atol(1)
 164  30   format(' Preconditioning uses interaction-only block-diagonal',
 165      1   ' matrix'/' with block-grouping, and Gauss-Seidel iterations'//
 166      2   ' Number of diagonal block groups = ngrp =',i4,
 167      3   '   (ngx by ngy, ngx =',i2,'  ngy =',i2,' )'//
 168      4   ' G-S preconditioner uses itmax iterations, itmax =',i3//
 169      5   ' Tolerance parameters: rtol =',d10.2,'   atol =',d10.2)
 170 c
 171 c
 172 c Loop over mf values 10, 21, ..., 24, 29.
 173 c
 174       do 90 mf = 10,29
 175       if (mf .gt. 10 .and. mf .lt. 21) go to 90
 176       if (mf .gt. 24 .and. mf .lt. 29) go to 90
 177       write(6,40)mf
 178  40   format(//80('-')//' Solution with mf =',i3//
 179      1   '   t       nstep  nfe  nni  nli  npe  nq',
 180      2   4x,'h          avdim    ncf rate    lcf rate')
 181 c
 182       t = 0.0d0
 183       tout = 1.0d-8
 184       nout = 18
 185       if (mf .eq. 10) nout = 6
 186       if (mf .eq. 23 .or. mf .eq. 24) nout = 10
 187       call cinit (cc)
 188       if (mf .eq. 22) call outweb (t, cc, ns, mx, my, 8)
 189       istate = 1
 190       nli = 0
 191       nni = 0
 192       ncfn = 0
 193       ncfl = 0
 194       nst = 0
 195 c
 196 c Loop over output times and call DLSODPK.
 197 c
 198       do 70 iout = 1,nout
 199         call dlsodpk (fweb, neq, cc, t, tout, itol, rtol, atol, itask,
 200      1         istate, iopt, rwork, lrw, iwork, liw, jacbg, solsbg, mf)
 201         nsdif = iwork(11) - nst
 202         nst = iwork(11)
 203         nfe = iwork(12)
 204         npe = iwork(13)
 205         nnidif = iwork(19) - nni
 206         nni = iwork(19)
 207         nlidif = iwork(20) - nli
 208         nli = iwork(20)
 209         nfndif = iwork(22) - ncfn
 210         ncfn = iwork(22)
 211         nfldif = iwork(23) - ncfl
 212         ncfl = iwork(23)
 213         nqu = iwork(14)
 214         hu = rwork(11)
 215         avdim = 0.0d0
 216         rcfn = 0.0d0
 217         rcfl = 0.0d0
 218         if (nnidif .gt. 0) avdim = real(nlidif)/real(nnidif)
 219         if (nsdif .gt. 0) rcfn = real(nfndif)/real(nsdif)
 220         if (nnidif .gt. 0) rcfl = real(nfldif)/real(nnidif)
 221         write(6,50)t,nst,nfe,nni,nli,npe,nqu,hu,avdim,rcfn,rcfl
 222  50     format(d10.2,i5,i6,3i5,i4,2d11.2,d10.2,d12.2)
 223         imod3 = iout - 3*(iout/3)
 224         if (mf .eq. 22 .and. imod3 .eq. 0) call outweb (t,cc,ns,mx,my,8)
 225         if (istate .eq. 2) go to 65
 226         write(6,60)t
 227  60     format(//' final time reached =',d12.4//)
 228         go to 75
 229  65     continue
 230         if (tout .gt. 0.9d0) tout = tout + 1.0d0
 231         if (tout .lt. 0.9d0) tout = tout*10.0d0
 232  70     continue
 233 c
 234  75   continue
 235       nst = iwork(11)
 236       nfe = iwork(12)
 237       npe = iwork(13)
 238       lenrw = iwork(17)
 239       leniw = iwork(18)
 240       nni = iwork(19)
 241       nli = iwork(20)
 242       nps = iwork(21)
 243       if (nni .gt. 0) avdim = real(nli)/real(nni)
 244       ncfn = iwork(22)
 245       ncfl = iwork(23)
 246       write (6,80) lenrw,leniw,nst,nfe,npe,nps,nni,nli,avdim,
 247      1               ncfn,ncfl
 248  80   format(//' Final statistics for this run:'/
 249      1   ' rwork size =',i8,'   iwork size =',i6/
 250      2   ' number of time steps            =',i5/
 251      3   ' number of f evaluations         =',i5/
 252      4   ' number of preconditioner evals. =',i5/
 253      4   ' number of preconditioner solves =',i5/
 254      5   ' number of nonlinear iterations  =',i5/
 255      5   ' number of linear iterations     =',i5/
 256      6   ' average subspace dimension  =',f8.4/
 257      7   i5,' nonlinear conv. failures,',i5,' linear conv. failures')
 258 c
 259  90   continue
 260 c      stop
 261 c------  end of main program for DLSODPK demonstration program ----------
 262       end
 263
 264       subroutine setpar
 265 c-----------------------------------------------------------------------
 266 c This routine sets the problem parameters.
 267 c It set ns, mx, my, and problem coefficients acoef, bcoef, diff, alph,
 268 c using parameters np, aa, ee, gg, bb, dprey, dpred.
 269 c-----------------------------------------------------------------------
 270       integer ns, mx, my, mxns
 271       integer i, j, np
 272       double precision aa, ee, gg, bb, dprey, dpred,
 273      1     ax, ay, acoef, bcoef, dx, dy, alph, diff, cox, coy
 274       common /pcom0/ aa, ee, gg, bb, dprey, dpred
 275       common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph,
 276      1               diff(20), cox(20), coy(20), ns, mx, my, mxns
 277 c
 278       np = 3
 279       mx = 6
 280       my = 6
 281       aa = 1.0d0
 282       ee = 1.0d4
 283       gg = 0.5d-6
 284       bb = 1.0d0
 285       dprey = 1.0d0
 286       dpred = 0.5d0
 287       alph = 1.0d0
 288       ns = 2*np
 289       do 70 j = 1,np
 290         do 60 i = 1,np
 291           acoef(np+i,j) = ee
 292           acoef(i,np+j) = -gg
 293  60       continue
 294         acoef(j,j) = -aa
 295         acoef(np+j,np+j) = -aa
 296         bcoef(j) = bb
 297         bcoef(np+j) = -bb
 298         diff(j) = dprey
 299         diff(np+j) = dpred
 300  70     continue
 301 c
 302       return
 303 c------------  end of subroutine setpar  -------------------------------
 304       end
 305
 306       subroutine gset (m, ng, jg, jig, jr)
 307 c-----------------------------------------------------------------------
 308 c This routine sets arrays jg, jig, and jr describing
 309 c a uniform partition of (1,2,...,m) into ng groups.
 310 c-----------------------------------------------------------------------
 311       integer m, ng, jg, jig, jr
 312       dimension jg(*), jig(*), jr(*)
 313       integer ig, j, len1, mper, ngm1
 314 c
 315       mper = m/ng
 316       do 10 ig = 1,ng
 317  10     jg(ig) = 1 + (ig - 1)*mper
 318       jg(ng+1) = m + 1
 319 c
 320       ngm1 = ng - 1
 321       len1 = ngm1*mper
 322       do 20 j = 1,len1
 323  20     jig(j) = 1 + (j-1)/mper
 324       len1 = len1 + 1
 325       do 25 j = len1,m
 326  25     jig(j) = ng
 327 c
 328       do 30 ig = 1,ngm1
 329  30     jr(ig) = 0.5d0 + (ig - 0.5d0)*mper
 330       jr(ng) = 0.5d0*(1 + ngm1*mper + m)
 331 c
 332       return
 333 c------------  end of subroutine gset  ---------------------------------
 334       end
 335
 336       subroutine cinit (cc)
 337 c-----------------------------------------------------------------------
 338 c This routine computes and loads the vector of initial values.
 339 c-----------------------------------------------------------------------
 340       double precision cc
 341       dimension cc(*)
 342       integer ns, mx, my, mxns
 343       integer i, ici, ioff, iyoff, jx, jy
 344       double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy
 345       double precision argx, argy, x, y
 346       common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph,
 347      1               diff(20), cox(20), coy(20), ns, mx, my, mxns
 348 c
 349         do 20 jy = 1,my
 350           y = (jy-1)*dy
 351           argy = 16.0d0*y*y*(ay-y)*(ay-y)
 352           iyoff = mxns*(jy-1)
 353           do 10 jx = 1,mx
 354             x = (jx-1)*dx
 355             argx = 16.0d0*x*x*(ax-x)*(ax-x)
 356             ioff = iyoff + ns*(jx-1)
 357             do 5 i = 1,ns
 358               ici = ioff + i
 359               cc(ici) = 10.0d0 + i*argx*argy
 360   5           continue
 361  10         continue
 362  20       continue
 363       return
 364 c------------  end of subroutine cinit  --------------------------------
 365       end
 366
 367       subroutine outweb (t, c, ns, mx, my, lun)
 368 c-----------------------------------------------------------------------
 369 c This routine prints the values of the individual species densities
 370 c at the current time t.  The write statements use unit lun.
 371 c-----------------------------------------------------------------------
 372       integer ns, mx, my, lun
 373       double precision t, c
 374       dimension c(ns,mx,my)
 375       integer i, jx, jy
 376 c
 377       write(lun,10) t
 378  10   format(/80('-')/30x,'At time t = ',d16.8/80('-') )
 379 c
 380       do 40 i = 1,ns
 381         write(lun,20) i
 382  20     format(' the species c(',i2,') values are:')
 383         do 30 jy = my,1,-1
 384           write(lun,25) (c(i,jx,jy),jx=1,mx)
 385  25       format(6(1x,g12.6))
 386  30       continue
 387         write(lun,35)
 388  35     format(80('-'),/)
 389  40     continue
 390 c
 391       return
 392 c------------  end of subroutine outweb  -------------------------------
 393       end
 394
 395       subroutine fweb (neq, t, cc, cdot)
 396 c-----------------------------------------------------------------------
 397 c This routine computes the derivative of cc and returns it in cdot.
 398 c The interaction rates are computed by calls to webr, and these are
 399 c saved in cc(neq+1),...,cc(2*neq) for use in preconditioning.
 400 c-----------------------------------------------------------------------
 401       integer neq(1)
 402       double precision t, cc, cdot
 403       dimension cc(*), cdot(*)
 404       integer ns, mx, my, mxns
 405       integer i, ic, ici, idxl, idxu, idyl, idyu, iyoff, jx, jy
 406       double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy
 407       double precision dcxli, dcxui, dcyli, dcyui, x, y
 408       common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph,
 409      1               diff(20), cox(20), coy(20), ns, mx, my, mxns
 410 c
 411       do 100 jy = 1,my
 412         y = (jy-1)*dy
 413         iyoff = mxns*(jy-1)
 414         idyu = mxns
 415         if (jy .eq. my) idyu = -mxns
 416         idyl = mxns
 417         if (jy .eq. 1) idyl = -mxns
 418         do 90 jx = 1,mx
 419           x = (jx-1)*dx
 420           ic = iyoff + ns*(jx-1) + 1
 421 c Get interaction rates at one point (x,y).
 422           call webr (x, y, t, cc(ic), cc(neq(1)+ic))
 423           idxu = ns
 424           if (jx .eq. mx) idxu = -ns
 425           idxl = ns
 426           if (jx .eq. 1) idxl = -ns
 427           do 80 i = 1,ns
 428             ici = ic + i - 1
 429 c Do differencing in y.
 430             dcyli = cc(ici) - cc(ici-idyl)
 431             dcyui = cc(ici+idyu) - cc(ici)
 432 c Do differencing in x.
 433             dcxli = cc(ici) - cc(ici-idxl)
 434             dcxui = cc(ici+idxu) - cc(ici)
 435 c Collect terms and load cdot elements.
 436             cdot(ici) = coy(i)*(dcyui - dcyli) + cox(i)*(dcxui - dcxli)
 437      1                  + cc(neq(1)+ici)
 438  80         continue
 439  90       continue
 440  100    continue
 441       return
 442 c------------  end of subroutine fweb  ---------------------------------
 443       end
 444
 445       subroutine webr (x, y, t, c, rate)
 446 c-----------------------------------------------------------------------
 447 c This routine computes the interaction rates for the species
 448 c c(1),...,c(ns), at one spatial point and at time t.
 449 c-----------------------------------------------------------------------
 450       double precision x, y, t, c, rate
 451       dimension c(*), rate(*)
 452       integer ns, mx, my, mxns
 453       integer i
 454       double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy
 455       double precision fac
 456       common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph,
 457      1               diff(20), cox(20), coy(20), ns, mx, my, mxns
 458 c
 459       do 10 i = 1,ns
 460  10     rate(i) = 0.0d0
 461       do 15 i = 1,ns
 462         call daxpy (ns, c(i), acoef(1,i), 1, rate, 1)
 463  15     continue
 464       fac = 1.0d0 + alph*x*y
 465       do 20 i = 1,ns
 466  20     rate(i) = c(i)*(bcoef(i)*fac + rate(i))
 467       return
 468 c------------  end of subroutine webr  ---------------------------------
 469       end
 470
 471       subroutine jacbg (f, neq, t, cc, ccsv, rewt, f0, f1, hl0,
 472      1                  bd, ipbd, ier)
 473 c-----------------------------------------------------------------------
 474 c This routine generates part of the block-diagonal part of the
 475 c Jacobian, multiplies by -hl0, adds the identity matrix,
 476 c and calls DGEFA to do LU decomposition of each diagonal block.
 477 c The computation of the diagonal blocks uses the block and grouping
 478 c information in /pcom1/ and /pcom2/.  One block per group is computed.
 479 c The Jacobian elements are generated by difference quotients
 480 c using calls to the routine fbg.
 481 c-----------------------------------------------------------------------
 482 c The two Common blocks below are used for internal communication.
 483 c The variables used are:
 484 c   mp     = size of blocks in block-diagonal preconditioning matrix.
 485 c   mq     = number of blocks in each direction (neq = mp*mq).
 486 c   mpsq   = mp*mp.
 487 c   uround = unit roundoff, generated by a call uround = dumach().
 488 c   srur   = sqrt(uround).
 489 c   meshx  = x mesh size
 490 c   meshy  = y mesh size (mesh is meshx by meshy)
 491 c   ngx    = no. groups in x direction in block-grouping scheme.
 492 c   ngy    = no. groups in y direction in block-grouping scheme.
 493 c   ngrp   = total number of groups = ngx*ngy.
 494 c   mxmp   = meshx*mp.
 495 c   jgx    = length ngx+1 array of group boundaries in x direction.
 496 c            group igx has x indices jx = jgx(igx),...,jgx(igx+1)-1.
 497 c   jigx   = length meshx array of x group indices vs x node index.
 498 c            x node index jx is in x group jigx(jx).
 499 c   jxr    = length ngx array of x indices representing the x groups.
 500 c            the index for x group igx is jx = jxr(igx).
 501 c   jgy, jigy, jyr = analogous arrays for grouping in y direction.
 502 c-----------------------------------------------------------------------
 503       external f
 504       integer neq, ipbd, ier
 505       double precision t, cc, ccsv, rewt, f0, f1, hl0, bd
 506       dimension cc(*), ccsv(*), rewt(*), f0(*), f1(*),
 507      1          bd(*), ipbd(*)
 508       dimension neq(1)
 509       integer mp, mq, mpsq, itmax,
 510      2        meshx,meshy,ngx,ngy,ngrp,mxmp,jgx,jgy,jigx,jigy,jxr,jyr
 511       integer i, ibd, idiag, if0, if00, ig, igx, igy, iip,
 512      1   j, jj, jx, jy, n
 513       double precision uround, srur
 514       double precision fac, r, r0, dvnorm
 515 c
 516       common /pcom2/ uround, srur, mp, mq, mpsq, itmax
 517       common /pcom3/ meshx, meshy, ngx, ngy, ngrp, mxmp,
 518      1   jgx(21), jgy(21), jigx(50), jigy(50), jxr(20), jyr(20)
 519 c
 520       n = neq(1)
 521 c
 522 c-----------------------------------------------------------------------
 523 c Make mp calls to fbg to approximate each diagonal block of Jacobian.
 524 c Here cc(neq+1),...,cc(2*neq) contains the base fb value.
 525 c r0 is a minimum increment factor for the difference quotient.
 526 c-----------------------------------------------------------------------
 527  200  fac = dvnorm (n, f0, rewt)
 528       r0 = 1000.0d0*abs(hl0)*uround*n*fac
 529       if (r0 .eq. 0.0d0) r0 = 1.0d0
 530       ibd = 0
 531       do 240 igy = 1,ngy
 532         jy = jyr(igy)
 533         if00 = (jy - 1)*mxmp
 534         do 230 igx = 1,ngx
 535           jx = jxr(igx)
 536           if0 = if00 + (jx - 1)*mp
 537           do 220 j = 1,mp
 538             jj = if0 + j
 539             r = max(srur*abs(cc(jj)),r0/rewt(jj))
 540             cc(jj) = cc(jj) + r
 541             fac = -hl0/r
 542             call fbg (neq, t, cc, jx, jy, f1)
 543             do 210 i = 1,mp
 544  210          bd(ibd+i) = (f1(i) - cc(neq(1)+if0+i))*fac
 545             cc(jj) = ccsv(jj)
 546             ibd = ibd + mp
 547  220        continue
 548  230      continue
 549  240    continue
 550 c
 551 c Add identity matrix and do LU decompositions on blocks. --------------
 552  260  continue
 553       ibd = 1
 554       iip = 1
 555       do 280 ig = 1,ngrp
 556         idiag = ibd
 557         do 270 i = 1,mp
 558           bd(idiag) = bd(idiag) + 1.0d0
 559  270      idiag = idiag + (mp + 1)
 560         call dgefa (bd(ibd), mp, mp, ipbd(iip), ier)
 561         if (ier .ne. 0) go to 290
 562         ibd = ibd + mpsq
 563         iip = iip + mp
 564  280    continue
 565  290  return
 566 c------------  end of subroutine jacbg  --------------------------------
 567       end
 568
 569       subroutine fbg (neq, t, cc, jx, jy, cdot)
 570 c-----------------------------------------------------------------------
 571 c This routine computes one block of the interaction terms of the
 572 c system, namely block (jx,jy), for use in preconditioning.
 573 c-----------------------------------------------------------------------
 574       integer neq, jx, jy
 575       double precision t, cc, cdot
 576       dimension cc(neq), cdot(neq)
 577       integer ns, mx, my, mxns
 578       integer iblok, ic
 579       double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy
 580       double precision x, y
 581 c
 582       common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph,
 583      1               diff(20), cox(20), coy(20), ns, mx, my, mxns
 584 c
 585       iblok = jx + (jy-1)*mx
 586         y = (jy-1)*dy
 587           x = (jx-1)*dx
 588           ic = ns*(iblok-1) + 1
 589           call webr (x, y, t, cc(ic), cdot)
 590       return
 591 c------------  end of subroutine fbg  ----------------------------------
 592       end
 593
 594       subroutine solsbg (n, t, cc, f0, wk, hl0, bd, ipbd, v, lr, ier)
 595 c-----------------------------------------------------------------------
 596 c This routine applies one or two inverse preconditioner matrices
 597 c to the array v, using the interaction-only block-diagonal Jacobian
 598 c with block-grouping, and Gauss-Seidel applied to the diffusion terms.
 599 c When lr = 1 or 3, it calls gs for a Gauss-Seidel approximation
 600 c to ((I-hl0*Jd)-inverse)*v, and stores the result in v.
 601 c When lr = 2 or 3, it computes ((I-hl0*dg/dc)-inverse)*v, using LU
 602 c factors of the blocks in bd, and pivot information in ipbd.
 603 c In both cases, the array v is overwritten with the solution.
 604 c-----------------------------------------------------------------------
 605       integer n, ipbd, lr, ier
 606       double precision t, cc, f0, wk, hl0, bd, v
 607       dimension cc(*), f0(*), wk(*), bd(*), ipbd(*), v(*)
 608       dimension n(1)
 609       integer mp, mq, mpsq, itmax,
 610      2        meshx,meshy,ngx,ngy,ngrp,mxmp,jgx,jgy,jigx,jigy,jxr,jyr
 611       integer ibd, ig0, igm1, igx, igy, iip, iv, jx, jy
 612       double precision uround, srur
 613 c
 614       common /pcom2/ uround, srur, mp, mq, mpsq, itmax
 615       common /pcom3/ meshx, meshy, ngx, ngy, ngrp, mxmp,
 616      1   jgx(21), jgy(21), jigx(50), jigy(50), jxr(20), jyr(20)
 617 c
 618       ier = 0
 619 c
 620       if (lr.eq.0 .or. lr.eq.1 .or. lr.eq.3) call gs (n, hl0, v, wk)
 621       if (lr.eq.0 .or. lr.eq.2 .or. lr.eq.3) then
 622         iv = 1
 623         do 20 jy = 1,meshy
 624           igy = jigy(jy)
 625           ig0 = (igy - 1)*ngx
 626           do 10 jx = 1,meshx
 627             igx = jigx(jx)
 628             igm1 = igx - 1 + ig0
 629             ibd = 1 + igm1*mpsq
 630             iip = 1 + igm1*mp
 631             call dgesl (bd(ibd), mp, mp, ipbd(iip), v(iv), 0)
 632             iv = iv + mp
 633  10         continue
 634  20       continue
 635         endif
 636 c
 637       return
 638 c------------  end of subroutine solsbg  -------------------------------
 639       end
 640
 641       subroutine gs (n, hl0, z, x)
 642 c-----------------------------------------------------------------------
 643 c This routine performs itmax Gauss-Seidel iterations
 644 c to compute an approximation to P-inverse*z, where P = I - hl0*Jd, and
 645 c Jd represents the diffusion contributions to the Jacobian.
 646 c z contains the answer on return.
 647 c The dimensions below assume ns .le. 20.
 648 c-----------------------------------------------------------------------
 649       integer n
 650       double precision hl0, z, x
 651       dimension z(*), x(*)
 652       dimension n(1)
 653       integer ns, mx, my, mxns,
 654      1        mp, mq, mpsq, itmax
 655       integer i, ic, ici, iter, iyoff, jx, jy
 656       double precision ax,ay,acoef,bcoef,dx,dy,alph,diff,cox,coy,
 657      2     uround, srur
 658       double precision beta,beta2,cof1,elamda,gamma,gamma2
 659       dimension beta(20), gamma(20), beta2(20), gamma2(20), cof1(20)
 660       common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph,
 661      1               diff(20), cox(20), coy(20), ns, mx, my, mxns
 662       common /pcom2/ uround, srur, mp, mq, mpsq, itmax
 663 c
 664 c-----------------------------------------------------------------------
 665 c Write matrix as P = D - L - U.
 666 c Load local arrays beta, beta2, gamma, gamma2, and cof1.
 667 c-----------------------------------------------------------------------
 668       do 10 i = 1,ns
 669         elamda = 1.d0/(1.d0 + 2.d0*hl0*(cox(i) + coy(i)))
 670         beta(i) = hl0*cox(i)*elamda
 671         beta2(i) = 2.d0*beta(i)
 672         gamma(i) = hl0*coy(i)*elamda
 673         gamma2(i) = 2.d0*gamma(i)
 674         cof1(i) = elamda
 675  10     continue
 676 c-----------------------------------------------------------------------
 677 c Begin iteration loop.
 678 c Load array x with (D-inverse)*z for first iteration.
 679 c-----------------------------------------------------------------------
 680       iter = 1
 681 c
 682       do 50 jy = 1,my
 683         iyoff = mxns*(jy-1)
 684         do 40 jx = 1,mx
 685           ic = iyoff + ns*(jx-1)
 686           do 30 i = 1,ns
 687             ici = ic + i
 688             x(ici) = cof1(i)*z(ici)
 689             z(ici) = 0.d0
 690  30         continue
 691  40       continue
 692  50     continue
 693       go to 160
 694 c-----------------------------------------------------------------------
 695 c Calculate (D-inverse)*U*x.
 696 c-----------------------------------------------------------------------
 697  70   continue
 698       iter = iter + 1
 699       jy = 1
 700         jx = 1
 701         ic = ns*(jx-1)
 702         do 75 i = 1,ns
 703           ici = ic + i
 704  75       x(ici) = beta2(i)*x(ici+ns) + gamma2(i)*x(ici+mxns)
 705         do 85 jx = 2,mx-1
 706           ic = ns*(jx-1)
 707           do 80 i = 1,ns
 708             ici = ic + i
 709  80         x(ici) = beta(i)*x(ici+ns) + gamma2(i)*x(ici+mxns)
 710  85       continue
 711         jx = mx
 712         ic = ns*(jx-1)
 713         do 90 i = 1,ns
 714           ici = ic + i
 715  90       x(ici) = gamma2(i)*x(ici+mxns)
 716       do 115 jy = 2,my-1
 717         iyoff = mxns*(jy-1)
 718           jx = 1
 719           ic = iyoff
 720           do 95 i = 1,ns
 721             ici = ic + i
 722  95         x(ici) = beta2(i)*x(ici+ns) + gamma(i)*x(ici+mxns)
 723           do 105 jx = 2,mx-1
 724             ic = iyoff + ns*(jx-1)
 725             do 100 i = 1,ns
 726               ici = ic + i
 727  100          x(ici) = beta(i)*x(ici+ns) + gamma(i)*x(ici+mxns)
 728  105        continue
 729           jx = mx
 730           ic = iyoff + ns*(jx-1)
 731           do 110 i = 1,ns
 732             ici = ic + i
 733  110        x(ici) = gamma(i)*x(ici+mxns)
 734  115      continue
 735       jy = my
 736       iyoff = mxns*(jy-1)
 737         jx = 1
 738         ic = iyoff
 739         do 120 i = 1,ns
 740           ici = ic + i
 741  120      x(ici) = beta2(i)*x(ici+ns)
 742         do 130 jx = 2,mx-1
 743           ic = iyoff + ns*(jx-1)
 744           do 125 i = 1,ns
 745             ici = ic + i
 746  125      x(ici) = beta(i)*x(ici+ns)
 747  130      continue
 748         jx = mx
 749         ic = iyoff + ns*(jx-1)
 750         do 135 i = 1,ns
 751           ici = ic + i
 752  135      x(ici) = 0.0d0
 753 c-----------------------------------------------------------------------
 754 c Calculate (I - (D-inverse)*L)-inverse * x.
 755 c-----------------------------------------------------------------------
 756  160  continue
 757       jy = 1
 758         do 175 jx = 2,mx-1
 759           ic = ns*(jx-1)
 760           do 170 i = 1,ns
 761             ici = ic + i
 762  170        x(ici) = x(ici) + beta(i)*x(ici-ns)
 763  175      continue
 764         jx = mx
 765         ic = ns*(jx-1)
 766         do 180 i = 1,ns
 767           ici = ic + i
 768  180      x(ici) = x(ici) + beta2(i)*x(ici-ns)
 769       do 210 jy = 2,my-1
 770         iyoff = mxns*(jy-1)
 771           jx = 1
 772           ic = iyoff
 773           do 185 i = 1,ns
 774             ici = ic + i
 775  185        x(ici) = x(ici) + gamma(i)*x(ici-mxns)
 776           do 200 jx = 2,mx-1
 777             ic = iyoff + ns*(jx-1)
 778             do 195 i = 1,ns
 779               ici = ic + i
 780               x(ici) = (x(ici) + beta(i)*x(ici-ns))
 781      1             + gamma(i)*x(ici-mxns)
 782  195          continue
 783  200        continue
 784             jx = mx
 785             ic = iyoff + ns*(jx-1)
 786             do 205 i = 1,ns
 787               ici = ic + i
 788               x(ici) = (x(ici) + beta2(i)*x(ici-ns))
 789      1             + gamma(i)*x(ici-mxns)
 790  205          continue
 791  210        continue
 792       jy = my
 793       iyoff = mxns*(jy-1)
 794         jx = 1
 795         ic = iyoff
 796         do 215 i = 1,ns
 797           ici = ic + i
 798  215      x(ici) = x(ici) + gamma2(i)*x(ici-mxns)
 799         do 225 jx = 2,mx-1
 800           ic = iyoff + ns*(jx-1)
 801           do 220 i = 1,ns
 802             ici = ic + i
 803             x(ici) = (x(ici) + beta(i)*x(ici-ns))
 804      1           + gamma2(i)*x(ici-mxns)
 805  220        continue
 806  225      continue
 807         jx = mx
 808         ic = iyoff + ns*(jx-1)
 809         do 230 i = 1,ns
 810           ici = ic + i
 811           x(ici) = (x(ici) + beta2(i)*x(ici-ns))
 812      1         + gamma2(i)*x(ici-mxns)
 813  230      continue
 814 c-----------------------------------------------------------------------
 815 c Add increment x to z.
 816 c-----------------------------------------------------------------------
 817       do 300 i = 1,n(1)
 818  300    z(i) = z(i) + x(i)
 819 c
 820       if (iter .lt. itmax) go to 70
 821       return
 822 c------------  end of subroutine gs  -----------------------------------
 823       end