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

Go to the documentation of this file.
00001       SUBROUTINE SPREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
00002      1   F, JAC)
00003 C***BEGIN PROLOGUE  SPREPJ
00004 C***SUBSIDIARY
00005 C***PURPOSE  Compute and process Newton iteration matrix.
00006 C***TYPE      SINGLE PRECISION (SPREPJ-S, DPREPJ-D)
00007 C***AUTHOR  Hindmarsh, Alan C., (LLNL)
00008 C***DESCRIPTION
00009 C
00010 C  SPREPJ is called by SSTODE to compute and process the matrix
00011 C  P = I - h*el(1)*J , where J is an approximation to the Jacobian.
00012 C  Here J is computed by the user-supplied routine JAC if
00013 C  MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.
00014 C  If MITER = 3, a diagonal approximation to J is used.
00015 C  J is stored in WM and replaced by P.  If MITER .ne. 3, P is then
00016 C  subjected to LU decomposition in preparation for later solution
00017 C  of linear systems with P as coefficient matrix.  This is done
00018 C  by SGETRF if MITER = 1 or 2, and by SGBTRF if MITER = 4 or 5.
00019 C
00020 C  In addition to variables described in SSTODE and SLSODE prologues,
00021 C  communication with SPREPJ uses the following:
00022 C  Y     = array containing predicted values on entry.
00023 C  FTEM  = work array of length N (ACOR in SSTODE).
00024 C  SAVF  = array containing f evaluated at predicted y.
00025 C  WM    = real work space for matrices.  On output it contains the
00026 C          inverse diagonal matrix if MITER = 3 and the LU decomposition
00027 C          of P if MITER is 1, 2 , 4, or 5.
00028 C          Storage of matrix elements starts at WM(3).
00029 C          WM also contains the following matrix-related data:
00030 C          WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
00031 C          WM(2) = H*EL0, saved for later use if MITER = 3.
00032 C  IWM   = integer work space containing pivot information, starting at
00033 C          IWM(21), if MITER is 1, 2, 4, or 5.  IWM also contains band
00034 C          parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
00035 C  EL0   = EL(1) (input).
00036 C  IERPJ = output error flag,  = 0 if no trouble, .gt. 0 if
00037 C          P matrix found to be singular.
00038 C  JCUR  = output flag = 1 to indicate that the Jacobian matrix
00039 C          (or approximation) is now current.
00040 C  This routine also uses the COMMON variables EL0, H, TN, UROUND,
00041 C  MITER, N, NFE, and NJE.
00042 C
00043 C***SEE ALSO  SLSODE
00044 C***ROUTINES CALLED  SGBTRF, SGETRF, SVNORM
00045 C***COMMON BLOCKS    SLS001
00046 C***REVISION HISTORY  (YYMMDD)
00047 C   791129  DATE WRITTEN
00048 C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
00049 C   890504  Minor cosmetic changes.  (FNF)
00050 C   930809  Renamed to allow single/double precision versions. (ACH)
00051 C   010412  Reduced size of Common block /SLS001/. (ACH)
00052 C   031105  Restored 'own' variables to Common block /SLS001/, to
00053 C           enable interrupt/restart feature. (ACH)
00054 C***END PROLOGUE  SPREPJ
00055 C**End
00056       EXTERNAL F, JAC
00057       INTEGER NEQ, NYH, IWM
00058       REAL Y, YH, EWT, FTEM, SAVF, WM
00059       DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
00060      1   WM(*), IWM(*)
00061       INTEGER IOWND, IOWNS,
00062      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
00063      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
00064      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
00065       REAL ROWNS,
00066      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
00067       COMMON /SLS001/ ROWNS(209),
00068      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
00069      2   IOWND(6), IOWNS(6),
00070      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
00071      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
00072      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
00073       INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
00074      1   MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
00075       REAL CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
00076      1   SVNORM
00077 C
00078 C***FIRST EXECUTABLE STATEMENT  SPREPJ
00079       NJE = NJE + 1
00080       IERPJ = 0
00081       JCUR = 1
00082       HL0 = H*EL0
00083       GO TO (100, 200, 300, 400, 500), MITER
00084 C If MITER = 1, call JAC and multiply by scalar. -----------------------
00085  100  LENP = N*N
00086       DO 110 I = 1,LENP
00087  110    WM(I+2) = 0.0E0
00088       CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N)
00089       CON = -HL0
00090       DO 120 I = 1,LENP
00091  120    WM(I+2) = WM(I+2)*CON
00092       GO TO 240
00093 C If MITER = 2, make N calls to F to approximate J. --------------------
00094  200  FAC = SVNORM (N, SAVF, EWT)
00095       R0 = 1000.0E0*ABS(H)*UROUND*N*FAC
00096       IF (R0 .EQ. 0.0E0) R0 = 1.0E0
00097       SRUR = WM(1)
00098       J1 = 2
00099       DO 230 J = 1,N
00100         YJ = Y(J)
00101         R = MAX(SRUR*ABS(YJ),R0/EWT(J))
00102         Y(J) = Y(J) + R
00103         FAC = -HL0/R
00104         CALL F (NEQ, TN, Y, FTEM)
00105         DO 220 I = 1,N
00106  220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
00107         Y(J) = YJ
00108         J1 = J1 + N
00109  230    CONTINUE
00110       NFE = NFE + N
00111 C Add identity matrix. -------------------------------------------------
00112  240  J = 3
00113       NP1 = N + 1
00114       DO 250 I = 1,N
00115         WM(J) = WM(J) + 1.0E0
00116  250    J = J + NP1
00117 C Do LU decomposition on P. --------------------------------------------
00118       CALL SGETRF (N, N, WM(3), N, IWM(21), IER)
00119       IF (IER .NE. 0) IERPJ = 1
00120       RETURN
00121 C If MITER = 3, construct a diagonal approximation to J and P. ---------
00122  300  WM(2) = HL0
00123       R = EL0*0.1E0
00124       DO 310 I = 1,N
00125  310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
00126       CALL F (NEQ, TN, Y, WM(3))
00127       NFE = NFE + 1
00128       DO 320 I = 1,N
00129         R0 = H*SAVF(I) - YH(I,2)
00130         DI = 0.1E0*R0 - H*(WM(I+2) - SAVF(I))
00131         WM(I+2) = 1.0E0
00132         IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320
00133         IF (ABS(DI) .EQ. 0.0E0) GO TO 330
00134         WM(I+2) = 0.1E0*R0/DI
00135  320    CONTINUE
00136       RETURN
00137  330  IERPJ = 1
00138       RETURN
00139 C If MITER = 4, call JAC and multiply by scalar. -----------------------
00140  400  ML = IWM(1)
00141       MU = IWM(2)
00142       ML3 = ML + 3
00143       MBAND = ML + MU + 1
00144       MEBAND = MBAND + ML
00145       LENP = MEBAND*N
00146       DO 410 I = 1,LENP
00147  410    WM(I+2) = 0.0E0
00148       CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
00149       CON = -HL0
00150       DO 420 I = 1,LENP
00151  420    WM(I+2) = WM(I+2)*CON
00152       GO TO 570
00153 C If MITER = 5, make MBAND calls to F to approximate J. ----------------
00154  500  ML = IWM(1)
00155       MU = IWM(2)
00156       MBAND = ML + MU + 1
00157       MBA = MIN(MBAND,N)
00158       MEBAND = MBAND + ML
00159       MEB1 = MEBAND - 1
00160       SRUR = WM(1)
00161       FAC = SVNORM (N, SAVF, EWT)
00162       R0 = 1000.0E0*ABS(H)*UROUND*N*FAC
00163       IF (R0 .EQ. 0.0E0) R0 = 1.0E0
00164       DO 560 J = 1,MBA
00165         DO 530 I = J,N,MBAND
00166           YI = Y(I)
00167           R = MAX(SRUR*ABS(YI),R0/EWT(I))
00168  530      Y(I) = Y(I) + R
00169         CALL F (NEQ, TN, Y, FTEM)
00170         DO 550 JJ = J,N,MBAND
00171           Y(JJ) = YH(JJ,1)
00172           YJJ = Y(JJ)
00173           R = MAX(SRUR*ABS(YJJ),R0/EWT(JJ))
00174           FAC = -HL0/R
00175           I1 = MAX(JJ-MU,1)
00176           I2 = MIN(JJ+ML,N)
00177           II = JJ*MEB1 - ML + 2
00178           DO 540 I = I1,I2
00179  540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
00180  550      CONTINUE
00181  560    CONTINUE
00182       NFE = NFE + MBA
00183 C Add identity matrix. -------------------------------------------------
00184  570  II = MBAND + 2
00185       DO 580 I = 1,N
00186         WM(II) = WM(II) + 1.0E0
00187  580    II = II + MEBAND
00188 C Do LU decomposition of P. --------------------------------------------
00189       CALL SGBTRF ( N, N, ML, MU, WM(3), MEBAND, IWM(21), IER)
00190       IF (IER .NE. 0) IERPJ = 1
00191       RETURN
00192 C----------------------- END OF SUBROUTINE SPREPJ ----------------------
00193       END