GNU Octave: libcruft/odepack/prepj.f Source File

Go to the documentation of this file.
00001       SUBROUTINE PREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
00002      1   F, JAC, IERR)
00003 CLLL. OPTIMIZE
00004       EXTERNAL F, JAC
00005       INTEGER NEQ, NYH, IWM
00006       INTEGER IOWND, IOWNS,
00007      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
00008      2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
00009       INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
00010      1   MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
00011       DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM 
00012       DOUBLE PRECISION ROWNS, 
00013      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
00014       DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
00015      1   VNORM
00016       DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
00017      1   WM(*), IWM(*)
00018       COMMON /LS0001/ ROWNS(209),
00019      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
00020      3   IOWND(14), IOWNS(6), 
00021      4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
00022      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
00023 C-----------------------------------------------------------------------
00024 C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX
00025 C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
00026 C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
00027 C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
00028 C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
00029 C J IS STORED IN WM AND REPLACED BY P.  IF MITER .NE. 3, P IS THEN
00030 C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
00031 C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
00032 C BY DGETRF IF MITER = 1 OR 2, AND BY DGBTRF IF MITER = 4 OR 5.
00033 C
00034 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
00035 C WITH PREPJ USES THE FOLLOWING..
00036 C Y     = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
00037 C FTEM  = WORK ARRAY OF LENGTH N (ACOR IN STODE). 
00038 C SAVF  = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
00039 C WM    = REAL WORK SPACE FOR MATRICES.  ON OUTPUT IT CONTAINS THE
00040 C         INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
00041 C         OF P IF MITER IS 1, 2 , 4, OR 5.
00042 C         STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
00043 C         WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
00044 C         WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
00045 C         WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
00046 C IWM   = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
00047 C         IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS BAND 
00048 C         PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
00049 C EL0   = EL(1) (INPUT).
00050 C IERPJ = OUTPUT ERROR FLAG,  = 0 IF NO TROUBLE, .GT. 0 IF
00051 C         P MATRIX FOUND TO BE SINGULAR.
00052 C JCUR  = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
00053 C         (OR APPROXIMATION) IS NOW CURRENT.
00054 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
00055 C MITER, N, NFE, AND NJE.
00056 C-----------------------------------------------------------------------
00057       NJE = NJE + 1 
00058       IERPJ = 0
00059       JCUR = 1
00060       HL0 = H*EL0
00061       GO TO (100, 200, 300, 400, 500), MITER
00062 C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
00063  100  LENP = N*N
00064       DO 110 I = 1,LENP
00065  110    WM(I+2) = 0.0D0
00066       CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N)
00067       CON = -HL0
00068       DO 120 I = 1,LENP
00069  120    WM(I+2) = WM(I+2)*CON 
00070       GO TO 240
00071 C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
00072  200  FAC = VNORM (N, SAVF, EWT)
00073       R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC  
00074       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
00075       SRUR = WM(1)
00076       J1 = 2
00077       DO 230 J = 1,N
00078         YJ = Y(J)
00079         R = DMAX1(SRUR*DABS(YJ),R0/EWT(J))
00080         Y(J) = Y(J) + R
00081         FAC = -HL0/R
00082         IERR = 0
00083         CALL F (NEQ, TN, Y, FTEM, IERR)
00084         IF (IERR .LT. 0) RETURN
00085         DO 220 I = 1,N
00086  220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
00087         Y(J) = YJ
00088         J1 = J1 + N 
00089  230    CONTINUE
00090       NFE = NFE + N 
00091 C ADD IDENTITY MATRIX. -------------------------------------------------
00092  240  J = 3
00093       NP1 = N + 1
00094       DO 250 I = 1,N
00095         WM(J) = WM(J) + 1.0D0 
00096  250    J = J + NP1 
00097 C DO LU DECOMPOSITION ON P. --------------------------------------------
00098       CALL DGETRF ( N, N, WM(3), N, IWM(21), IER)
00099       IF (IER .NE. 0) IERPJ = 1
00100       RETURN
00101 C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
00102  300  WM(2) = HL0
00103       R = EL0*0.1D0 
00104       DO 310 I = 1,N
00105  310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
00106       IERR = 0
00107       CALL F (NEQ, TN, Y, WM(3), IERR)
00108       IF (IERR .LT. 0) RETURN
00109       NFE = NFE + 1 
00110       DO 320 I = 1,N
00111         R0 = H*SAVF(I) - YH(I,2)
00112         DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
00113         WM(I+2) = 1.0D0
00114         IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320
00115         IF (DABS(DI) .EQ. 0.0D0) GO TO 330
00116         WM(I+2) = 0.1D0*R0/DI 
00117  320    CONTINUE
00118       RETURN
00119  330  IERPJ = 1
00120       RETURN
00121 C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
00122  400  ML = IWM(1)
00123       MU = IWM(2)
00124       ML3 = ML + 3
00125       MBAND = ML + MU + 1
00126       MEBAND = MBAND + ML
00127       LENP = MEBAND*N
00128       DO 410 I = 1,LENP
00129  410    WM(I+2) = 0.0D0
00130       CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
00131       CON = -HL0
00132       DO 420 I = 1,LENP
00133  420    WM(I+2) = WM(I+2)*CON 
00134       GO TO 570
00135 C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
00136  500  ML = IWM(1)
00137       MU = IWM(2)
00138       MBAND = ML + MU + 1
00139       MBA = MIN0(MBAND,N)
00140       MEBAND = MBAND + ML
00141       MEB1 = MEBAND - 1
00142       SRUR = WM(1)
00143       FAC = VNORM (N, SAVF, EWT)
00144       R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC  
00145       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
00146       DO 560 J = 1,MBA
00147         DO 530 I = J,N,MBAND
00148           YI = Y(I) 
00149           R = DMAX1(SRUR*DABS(YI),R0/EWT(I))
00150  530      Y(I) = Y(I) + R
00151         IERR = 0
00152         CALL F (NEQ, TN, Y, FTEM, IERR)
00153         IF (IERR .LT. 0) RETURN
00154         DO 550 JJ = J,N,MBAND 
00155           Y(JJ) = YH(JJ,1)
00156           YJJ = Y(JJ)
00157           R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ))
00158           FAC = -HL0/R
00159           I1 = MAX0(JJ-MU,1)
00160           I2 = MIN0(JJ+ML,N)
00161           II = JJ*MEB1 - ML + 2
00162           DO 540 I = I1,I2
00163  540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
00164  550      CONTINUE
00165  560    CONTINUE
00166       NFE = NFE + MBA
00167 C ADD IDENTITY MATRIX. -------------------------------------------------
00168  570  II = MBAND + 2
00169       DO 580 I = 1,N
00170         WM(II) = WM(II) + 1.0D0
00171  580    II = II + MEBAND
00172 C DO LU DECOMPOSITION OF P. --------------------------------------------
00173       CALL DGBTRF ( N, N, ML, MU, WM(3), MEBAND, IWM(21), IER)
00174       IF (IER .NE. 0) IERPJ = 1
00175       RETURN
00176 C----------------------- END OF SUBROUTINE PREPJ -----------------------
00177       END