prepj.f

Go to the documentation of this file.

      SUBROUTINE prepj (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
     1   f, jac, ierr)
CLLL. OPTIMIZE
      EXTERNAL f, jac
      INTEGER neq, nyh, iwm
      INTEGER iownd, iowns,
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     2   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
      DOUBLE PRECISION y, yh, ewt, ftem, savf, wm 
      DOUBLE PRECISION rowns, 
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      DOUBLE PRECISION con, di, fac, hl0, r, r0, srur, yi, yj, yjj,
     1   vnorm
      dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
     1   wm(*), iwm(*)
      COMMON /ls0001/ rowns(209),
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     3   iownd(14), iowns(6), 
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
C-----------------------------------------------------------------------
C PREPJ IS CALLED BY STODE 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 DGETRF IF MITER = 1 OR 2, AND BY DGBTRF IF MITER = 4 OR 5.
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C WITH PREPJ USES THE FOLLOWING..
C Y     = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
C FTEM  = WORK ARRAY OF LENGTH N (ACOR IN STODE). 
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-----------------------------------------------------------------------
      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.0d0
      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 = vnorm(n, savf, ewt)
      r0 = 1000.0d0*dabs(h)*uround*dble(n)*fac  
      IF (r0 .EQ. 0.0d0) r0 = 1.0d0
      srur = wm(1)
      j1 = 2
      DO 230 j = 1,n
        yj = y(j)
        r = dmax1(srur*dabs(yj),r0/ewt(j))
        y(j) = y(j) + r
        fac = -hl0/r
        ierr = 0
        CALL f(neq, tn, y, ftem, ierr)
        IF (ierr .LT. 0) RETURN
        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.0d0 
 250    j = j + np1 
C DO LU DECOMPOSITION ON P. --------------------------------------------
      CALL dgetrf( 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.1d0 
      DO 310 i = 1,n
 310    y(i) = y(i) + r*(h*savf(i) - yh(i,2))
      ierr = 0
      CALL f(neq, tn, y, wm(3), ierr)
      IF (ierr .LT. 0) RETURN
      nfe = nfe + 1 
      DO 320 i = 1,n
        r0 = h*savf(i) - yh(i,2)
        di = 0.1d0*r0 - h*(wm(i+2) - savf(i))
        wm(i+2) = 1.0d0
        IF (dabs(r0) .LT. uround/ewt(i)) go to 320
        IF (dabs(di) .EQ. 0.0d0) go to 330
        wm(i+2) = 0.1d0*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.0d0
      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 = min0(mband,n)
      meband = mband + ml
      meb1 = meband - 1
      srur = wm(1)
      fac = vnorm(n, savf, ewt)
      r0 = 1000.0d0*dabs(h)*uround*dble(n)*fac  
      IF (r0 .EQ. 0.0d0) r0 = 1.0d0
      DO 560 j = 1,mba
        DO 530 i = j,n,mband
          yi = y(i) 
          r = dmax1(srur*dabs(yi),r0/ewt(i))
 530      y(i) = y(i) + r
        ierr = 0
        CALL f(neq, tn, y, ftem, ierr)
        IF (ierr .LT. 0) RETURN
        DO 550 jj = j,n,mband 
          y(jj) = yh(jj,1)
          yjj = y(jj)
          r = dmax1(srur*dabs(yjj),r0/ewt(jj))
          fac = -hl0/r
          i1 = max0(jj-mu,1)
          i2 = min0(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.0d0
 580    ii = ii + meband
C DO LU DECOMPOSITION OF P. --------------------------------------------
      CALL dgbtrf( n, n, ml, mu, wm(3), meband, iwm(21), ier)
      IF (ier .NE. 0) ierpj = 1
      RETURN
C----------------------- END OF SUBROUTINE PREPJ -----------------------
      END