relibc_openlibm/slatec/dpjac.f

*DECK DPJAC
      SUBROUTINE DPJAC (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, DF,
     +   DJAC, RPAR, IPAR)
C***BEGIN PROLOGUE  DPJAC
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DDEBDF
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (PJAC-S, DPJAC-D)
C***AUTHOR  Watts, H. A., (SNLA)
C***DESCRIPTION
C
C   DPJAC sets up the iteration matrix (involving the Jacobian) for the
C   integration package DDEBDF.
C
C***SEE ALSO  DDEBDF
C***ROUTINES CALLED  DGBFA, DGEFA, DVNRMS
C***COMMON BLOCKS    DDEBD1
C***REVISION HISTORY  (YYMMDD)
C   820301  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910722  Updated AUTHOR section.  (ALS)
C   920422  Changed DIMENSION statement.  (WRB)
C***END PROLOGUE  DPJAC
C
      INTEGER I, I1, I2, IER, II, IOWND, IOWNS, IPAR, IWM, J, J1,
     1      JJ, JSTART, KFLAG, L, LENP, MAXORD, MBA, MBAND,
     2      MEB1, MEBAND, METH, MITER, ML, ML3, MU, N, NEQ,
     3      NFE, NJE, NQ, NQU, NST, NYH
      DOUBLE PRECISION CON, DI, DVNRMS, EL0, EWT,
     1      FAC, FTEM, H, HL0, HMIN, HMXI, HU, R, R0, ROWND, ROWNS,
     2      RPAR, SAVF, SRUR, TN, UROUND, WM, Y, YH, YI, YJ, YJJ
      EXTERNAL DF, DJAC
      DIMENSION Y(*),YH(NYH,*),EWT(*),FTEM(*),SAVF(*),WM(*),IWM(*),
     1          RPAR(*),IPAR(*)
      COMMON /DDEBD1/ ROWND,ROWNS(210),EL0,H,HMIN,HMXI,HU,TN,UROUND,
     1                IOWND(14),IOWNS(6),IER,JSTART,KFLAG,L,METH,MITER,
     2                MAXORD,N,NQ,NST,NFE,NJE,NQU
C     ------------------------------------------------------------------
C      DPJAC IS CALLED BY DSTOD  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 DJAC 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 DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5.
C
C      IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C      WITH DPJAC USES THE FOLLOWING..
C      Y    = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
C      FTEM = WORK ARRAY OF LENGTH N (ACOR IN DSTOD ).
C      SAVF = ARRAY CONTAINING DF EVALUATED AT PREDICTED Y.
C      WM   = DOUBLE PRECISION WORK SPACE FOR MATRICES.  ON OUTPUT IT
C      CONTAINS THE
C             INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU
C             DECOMPOSITION 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
C             INCREMENTS.  WM(2) = H*EL0, SAVED FOR LATER USE IF MITER =
C             3.
C      IWM  = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING
C             AT IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS
C             THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER
C             IS 4 OR 5.
C      EL0  = EL(1) (INPUT).
C      IER  = OUTPUT ERROR FLAG,  = 0 IF NO TROUBLE, .NE. 0 IF
C             P MATRIX FOUND TO BE SINGULAR.
C      THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
C      MITER, N, NFE, AND NJE.
C-----------------------------------------------------------------------
C     BEGIN BLOCK PERMITTING ...EXITS TO 240
C        BEGIN BLOCK PERMITTING ...EXITS TO 220
C           BEGIN BLOCK PERMITTING ...EXITS TO 130
C              BEGIN BLOCK PERMITTING ...EXITS TO 70
C***FIRST EXECUTABLE STATEMENT  DPJAC
                  NJE = NJE + 1
                  HL0 = H*EL0
                  GO TO (10,40,90,140,170), MITER
C                 IF MITER = 1, CALL DJAC AND MULTIPLY BY SCALAR.
C                 -----------------------
   10             CONTINUE
                  LENP = N*N
                  DO 20 I = 1, LENP
                     WM(I+2) = 0.0D0
   20             CONTINUE
                  CALL DJAC(TN,Y,WM(3),N,RPAR,IPAR)
                  CON = -HL0
                  DO 30 I = 1, LENP
                     WM(I+2) = WM(I+2)*CON
   30             CONTINUE
C              ...EXIT
                  GO TO 70
C                 IF MITER = 2, MAKE N CALLS TO DF TO APPROXIMATE J.
C                 --------------------
   40             CONTINUE
                  FAC = DVNRMS(N,SAVF,EWT)
                  R0 = 1000.0D0*ABS(H)*UROUND*N*FAC
                  IF (R0 .EQ. 0.0D0) R0 = 1.0D0
                  SRUR = WM(1)
                  J1 = 2
                  DO 60 J = 1, N
                     YJ = Y(J)
                     R = MAX(SRUR*ABS(YJ),R0*EWT(J))
                     Y(J) = Y(J) + R
                     FAC = -HL0/R
                     CALL DF(TN,Y,FTEM,RPAR,IPAR)
                     DO 50 I = 1, N
                        WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
   50                CONTINUE
                     Y(J) = YJ
                     J1 = J1 + N
   60             CONTINUE
                  NFE = NFE + N
   70          CONTINUE
C              ADD IDENTITY MATRIX.
C              -------------------------------------------------
               J = 3
               DO 80 I = 1, N
                  WM(J) = WM(J) + 1.0D0
                  J = J + (N + 1)
   80          CONTINUE
C              DO LU DECOMPOSITION ON P.
C              --------------------------------------------
               CALL DGEFA(WM(3),N,N,IWM(21),IER)
C     .........EXIT
               GO TO 240
C              IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND
C              P. ---------
   90          CONTINUE
               WM(2) = HL0
               IER = 0
               R = EL0*0.1D0
               DO 100 I = 1, N
                  Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
  100          CONTINUE
               CALL DF(TN,Y,WM(3),RPAR,IPAR)
               NFE = NFE + 1
               DO 120 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 (ABS(R0) .LT. UROUND*EWT(I)) GO TO 110
C           .........EXIT
                     IF (ABS(DI) .EQ. 0.0D0) GO TO 130
                     WM(I+2) = 0.1D0*R0/DI
  110             CONTINUE
  120          CONTINUE
C     .........EXIT
               GO TO 240
  130       CONTINUE
            IER = -1
C     ......EXIT
            GO TO 240
C           IF MITER = 4, CALL DJAC AND MULTIPLY BY SCALAR.
C           -----------------------
  140       CONTINUE
            ML = IWM(1)
            MU = IWM(2)
            ML3 = 3
            MBAND = ML + MU + 1
            MEBAND = MBAND + ML
            LENP = MEBAND*N
            DO 150 I = 1, LENP
               WM(I+2) = 0.0D0
  150       CONTINUE
            CALL DJAC(TN,Y,WM(ML3),MEBAND,RPAR,IPAR)
            CON = -HL0
            DO 160 I = 1, LENP
               WM(I+2) = WM(I+2)*CON
  160       CONTINUE
C        ...EXIT
            GO TO 220
C           IF MITER = 5, MAKE MBAND CALLS TO DF TO APPROXIMATE J.
C           ----------------
  170       CONTINUE
            ML = IWM(1)
            MU = IWM(2)
            MBAND = ML + MU + 1
            MBA = MIN(MBAND,N)
            MEBAND = MBAND + ML
            MEB1 = MEBAND - 1
            SRUR = WM(1)
            FAC = DVNRMS(N,SAVF,EWT)
            R0 = 1000.0D0*ABS(H)*UROUND*N*FAC
            IF (R0 .EQ. 0.0D0) R0 = 1.0D0
            DO 210 J = 1, MBA
               DO 180 I = J, N, MBAND
                  YI = Y(I)
                  R = MAX(SRUR*ABS(YI),R0*EWT(I))
                  Y(I) = Y(I) + R
  180          CONTINUE
               CALL DF(TN,Y,FTEM,RPAR,IPAR)
               DO 200 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 190 I = I1, I2
                     WM(II+I) = (FTEM(I) - SAVF(I))*FAC
  190             CONTINUE
  200          CONTINUE
  210       CONTINUE
            NFE = NFE + MBA
  220    CONTINUE
C        ADD IDENTITY MATRIX.
C        -------------------------------------------------
         II = MBAND + 2
         DO 230 I = 1, N
            WM(II) = WM(II) + 1.0D0
            II = II + MEBAND
  230    CONTINUE
C        DO LU DECOMPOSITION OF P.
C        --------------------------------------------
         CALL DGBFA(WM(3),MEBAND,N,ML,MU,IWM(21),IER)
  240 CONTINUE
      RETURN
C     ----------------------- END OF SUBROUTINE DPJAC
C     -----------------------
      END