sprepj.f

Go to the documentation of this file.

      SUBROUTINE sprepj (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
     1   f, jac)
C***BEGIN PROLOGUE  SPREPJ
C***SUBSIDIARY
C***PURPOSE  Compute and process Newton iteration matrix.
C***TYPE      SINGLE PRECISION (SPREPJ-S, DPREPJ-D)
C***AUTHOR  Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C  SPREPJ is called by SSTODE to compute and process the matrix
C  P = I - h*el(1)*J , where J is an approximation to the Jacobian.
C  Here J is computed by the user-supplied routine JAC if
C  MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.
C  If MITER = 3, a diagonal approximation to J is used.
C  J is stored in WM and replaced by P.  If MITER .ne. 3, P is then
C  subjected to LU decomposition in preparation for later solution
C  of linear systems with P as coefficient matrix.  This is done
C  by SGETRF if MITER = 1 or 2, and by SGBTRF if MITER = 4 or 5.
C
C  In addition to variables described in SSTODE and SLSODE prologues,
C  communication with SPREPJ uses the following:
C  Y     = array containing predicted values on entry.
C  FTEM  = work array of length N (ACOR in SSTODE).
C  SAVF  = array containing f evaluated at predicted y.
C  WM    = real work space for matrices.  On output it contains the
C          inverse diagonal matrix if MITER = 3 and the LU decomposition
C          of P if MITER is 1, 2 , 4, or 5.
C          Storage of matrix elements starts at WM(3).
C          WM also contains the following matrix-related data:
C          WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
C          WM(2) = H*EL0, saved for later use if MITER = 3.
C  IWM   = integer work space containing pivot information, starting at
C          IWM(21), if MITER is 1, 2, 4, or 5.  IWM also contains band
C          parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
C  EL0   = EL(1) (input).
C  IERPJ = output error flag,  = 0 if no trouble, .gt. 0 if
C          P matrix found to be singular.
C  JCUR  = output flag = 1 to indicate that the Jacobian matrix
C          (or approximation) is now current.
C  This routine also uses the COMMON variables EL0, H, TN, UROUND,
C  MITER, N, NFE, and NJE.
C
C***SEE ALSO  SLSODE
C***ROUTINES CALLED  SGBTRF, SGETRF, SVNORM
C***COMMON BLOCKS    SLS001
C***REVISION HISTORY  (YYMMDD)
C   791129  DATE WRITTEN
C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
C   890504  Minor cosmetic changes.  (FNF)
C   930809  Renamed to allow single/double precision versions. (ACH)
C   010412  Reduced size of Common block /SLS001/. (ACH)
C   031105  Restored 'own' variables to Common block /SLS001/, to
C           enable interrupt/restart feature. (ACH)
C***END PROLOGUE  SPREPJ
C**End
      EXTERNAL f, jac
      INTEGER NEQ, NYH, IWM
      REAL Y, YH, EWT, FTEM, SAVF, WM
      dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
     1   wm(*), iwm(*)
      INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH,
     1   ialth, ipup, lmax, meo, nqnyh, nslp,
     1   icf, ierpj, iersl, jcur, jstart, kflag, l,
     2   lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
     3   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      COMMON /sls001/ conit, crate, el(13), elco(13,12),
     1   hold, rmax, tesco(3,12),
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     2   init, mxstep, mxhnil, nhnil, nslast, cnyh,
     3   ialth, ipup, lmax, meo, nqnyh, nslp,
     3   icf, ierpj, iersl, jcur, jstart, kflag, l,
     4   lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
     1   mba, mband, meb1, meband, ml, ml3, mu, np1
      REAL CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
     1   svnorm
C
C***FIRST EXECUTABLE STATEMENT  SPREPJ
      nje = nje + 1
      ierpj = 0
      jcur = 1
      hl0 = h*el0
      go to(100, 200, 300, 400, 500), miter
C If MITER = 1, call JAC and multiply by scalar. -----------------------
 100  lenp = n*n
      DO 110 i = 1,lenp
 110    wm(i+2) = 0.0e0
      CALL jac(neq, tn, y, 0, 0, wm(3), n)
      con = -hl0
      DO 120 i = 1,lenp
 120    wm(i+2) = wm(i+2)*con
      go to 240
C If MITER = 2, make N calls to F to approximate J. --------------------
 200  fac = svnorm(n, savf, ewt)
      r0 = 1000.0e0*abs(h)*uround*n*fac
      IF (r0 .EQ. 0.0e0) r0 = 1.0e0
      srur = wm(1)
      j1 = 2
      DO 230 j = 1,n
        yj = y(j)
        r = max(srur*abs(yj),r0/ewt(j))
        y(j) = y(j) + r
        fac = -hl0/r
        CALL f(neq, tn, y, ftem)
        DO 220 i = 1,n
 220      wm(i+j1) = (ftem(i) - savf(i))*fac
        y(j) = yj
        j1 = j1 + n
 230    CONTINUE
      nfe = nfe + n
C Add identity matrix. -------------------------------------------------
 240  j = 3
      np1 = n + 1
      DO 250 i = 1,n
        wm(j) = wm(j) + 1.0e0
 250    j = j + np1
C Do LU decomposition on P. --------------------------------------------
      CALL sgetrf(n, n, wm(3), n, iwm(21), ier)
      IF (ier .NE. 0) ierpj = 1
      RETURN
C If MITER = 3, construct a diagonal approximation to J and P. ---------
 300  wm(2) = hl0
      r = el0*0.1e0
      DO 310 i = 1,n
 310    y(i) = y(i) + r*(h*savf(i) - yh(i,2))
      CALL f(neq, tn, y, wm(3))
      nfe = nfe + 1
      DO 320 i = 1,n
        r0 = h*savf(i) - yh(i,2)
        di = 0.1e0*r0 - h*(wm(i+2) - savf(i))
        wm(i+2) = 1.0e0
        IF (abs(r0) .LT. uround/ewt(i)) go to 320
        IF (abs(di) .EQ. 0.0e0) go to 330
        wm(i+2) = 0.1e0*r0/di
 320    CONTINUE
      RETURN
 330  ierpj = 1
      RETURN
C If MITER = 4, call JAC and multiply by scalar. -----------------------
 400  ml = iwm(1)
      mu = iwm(2)
      ml3 = ml + 3
      mband = ml + mu + 1
      meband = mband + ml
      lenp = meband*n
      DO 410 i = 1,lenp
 410    wm(i+2) = 0.0e0
      CALL jac(neq, tn, y, ml, mu, wm(ml3), meband)
      con = -hl0
      DO 420 i = 1,lenp
 420    wm(i+2) = wm(i+2)*con
      go to 570
C If MITER = 5, make MBAND calls to F to approximate J. ----------------
 500  ml = iwm(1)
      mu = iwm(2)
      mband = ml + mu + 1
      mba = min(mband,n)
      meband = mband + ml
      meb1 = meband - 1
      srur = wm(1)
      fac = svnorm(n, savf, ewt)
      r0 = 1000.0e0*abs(h)*uround*n*fac
      IF (r0 .EQ. 0.0e0) r0 = 1.0e0
      DO 560 j = 1,mba
        DO 530 i = j,n,mband
          yi = y(i)
          r = max(srur*abs(yi),r0/ewt(i))
 530      y(i) = y(i) + r
        CALL f(neq, tn, y, ftem)
        DO 550 jj = j,n,mband
          y(jj) = yh(jj,1)
          yjj = y(jj)
          r = max(srur*abs(yjj),r0/ewt(jj))
          fac = -hl0/r
          i1 = max(jj-mu,1)
          i2 = min(jj+ml,n)
          ii = jj*meb1 - ml + 2
          DO 540 i = i1,i2
 540        wm(ii+i) = (ftem(i) - savf(i))*fac
 550      CONTINUE
 560    CONTINUE
      nfe = nfe + mba
C Add identity matrix. -------------------------------------------------
 570  ii = mband + 2
      DO 580 i = 1,n
        wm(ii) = wm(ii) + 1.0e0
 580    ii = ii + meband
C Do LU decomposition of P. --------------------------------------------
      CALL sgbtrf( n, n, ml, mu, wm(3), meband, iwm(21), ier)
      IF (ier .NE. 0) ierpj = 1
      RETURN
C----------------------- END OF SUBROUTINE SPREPJ ----------------------
      END