share/odepack/fortran/dsolpk.f

   1 *DECK DSOLPK
   2       SUBROUTINE DSOLPK (NEQ, Y, SAVF, X, EWT, WM, IWM, F, PSOL)
   3       EXTERNAL F, PSOL
   4       INTEGER NEQ, IWM
   5       DOUBLE PRECISION Y, SAVF, X, EWT, WM
   6       DIMENSION NEQ(*), Y(*), SAVF(*), X(*), EWT(*), WM(*), IWM(*)
   7       INTEGER IOWND, IOWNS,
   8      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
   9      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  10      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  11       INTEGER JPRE, JACFLG, LOCWP, LOCIWP, LSAVX, KMP, MAXL, MNEWT,
  12      1   NNI, NLI, NPS, NCFN, NCFL
  13       DOUBLE PRECISION ROWNS,
  14      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  15       DOUBLE PRECISION DELT, EPCON, SQRTN, RSQRTN
  16       COMMON /DLS001/ ROWNS(209),
  17      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  18      2   IOWND(6), IOWNS(6),
  19      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  20      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  21      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  22       COMMON /DLPK01/ DELT, EPCON, SQRTN, RSQRTN,
  23      1   JPRE, JACFLG, LOCWP, LOCIWP, LSAVX, KMP, MAXL, MNEWT,
  24      2   NNI, NLI, NPS, NCFN, NCFL
  25 C-----------------------------------------------------------------------
  26 C This routine interfaces to one of DSPIOM, DSPIGMR, DPCG, DPCGS, or
  27 C DUSOL, for the solution of the linear system arising from a Newton
  28 C iteration.  It is called if MITER .ne. 0.
  29 C In addition to variables described elsewhere,
  30 C communication with DSOLPK uses the following variables:
  31 C WM    = real work space containing data for the algorithm
  32 C         (Krylov basis vectors, Hessenberg matrix, etc.)
  33 C IWM   = integer work space containing data for the algorithm
  34 C X     = the right-hand side vector on input, and the solution vector
  35 C         on output, of length N.
  36 C IERSL = output flag (in Common):
  37 C         IERSL =  0 means no trouble occurred.
  38 C         IERSL =  1 means the iterative method failed to converge.
  39 C                    If the preconditioner is out of date, the step
  40 C                    is repeated with a new preconditioner.
  41 C                    Otherwise, the stepsize is reduced (forcing a
  42 C                    new evaluation of the preconditioner) and the
  43 C                    step is repeated.
  44 C         IERSL = -1 means there was a nonrecoverable error in the
  45 C                    iterative solver, and an error exit occurs.
  46 C This routine also uses the Common variables TN, EL0, H, N, MITER,
  47 C DELT, EPCON, SQRTN, RSQRTN, MAXL, KMP, MNEWT, NNI, NLI, NPS, NCFL,
  48 C LOCWP, LOCIWP.
  49 C-----------------------------------------------------------------------
  50       INTEGER IFLAG, LB, LDL, LHES, LIOM, LGMR, LPCG, LP, LQ, LR,
  51      1   LV, LW, LWK, LZ, MAXLP1, NPSL
  52       DOUBLE PRECISION DELTA, HL0
  53 C
  54       IERSL = 0
  55       HL0 = H*EL0
  56       DELTA = DELT*EPCON
  57       GO TO (100, 200, 300, 400, 900, 900, 900, 900, 900), MITER
  58 C-----------------------------------------------------------------------
  59 C Use the SPIOM algorithm to solve the linear system P*x = -f.
  60 C-----------------------------------------------------------------------
  61  100  CONTINUE
  62       LV = 1
  63       LB = LV + N*MAXL
  64       LHES = LB + N
  65       LWK = LHES + MAXL*MAXL
  66       CALL DCOPY (N, X, 1, WM(LB), 1)
  67       CALL DSCAL (N, RSQRTN, EWT, 1)
  68       CALL DSPIOM (NEQ, TN, Y, SAVF, WM(LB), EWT, N, MAXL, KMP, DELTA,
  69      1   HL0, JPRE, MNEWT, F, PSOL, NPSL, X, WM(LV), WM(LHES), IWM,
  70      2   LIOM, WM(LOCWP), IWM(LOCIWP), WM(LWK), IFLAG)
  71       NNI = NNI + 1
  72       NLI = NLI + LIOM
  73       NPS = NPS + NPSL
  74       CALL DSCAL (N, SQRTN, EWT, 1)
  75       IF (IFLAG .NE. 0) NCFL = NCFL + 1
  76       IF (IFLAG .GE. 2) IERSL = 1
  77       IF (IFLAG .LT. 0) IERSL = -1
  78       RETURN
  79 C-----------------------------------------------------------------------
  80 C Use the SPIGMR algorithm to solve the linear system P*x = -f.
  81 C-----------------------------------------------------------------------
  82  200  CONTINUE
  83       MAXLP1 = MAXL + 1
  84       LV = 1
  85       LB = LV + N*MAXL
  86       LHES = LB + N + 1
  87       LQ = LHES + MAXL*MAXLP1
  88       LWK = LQ + 2*MAXL
  89       LDL = LWK + MIN(1,MAXL-KMP)*N
  90       CALL DCOPY (N, X, 1, WM(LB), 1)
  91       CALL DSCAL (N, RSQRTN, EWT, 1)
  92       CALL DSPIGMR (NEQ, TN, Y, SAVF, WM(LB), EWT, N, MAXL, MAXLP1, KMP,
  93      1   DELTA, HL0, JPRE, MNEWT, F, PSOL, NPSL, X, WM(LV), WM(LHES),
  94      2   WM(LQ), LGMR, WM(LOCWP), IWM(LOCIWP), WM(LWK), WM(LDL), IFLAG)
  95       NNI = NNI + 1
  96       NLI = NLI + LGMR
  97       NPS = NPS + NPSL
  98       CALL DSCAL (N, SQRTN, EWT, 1)
  99       IF (IFLAG .NE. 0) NCFL = NCFL + 1
 100       IF (IFLAG .GE. 2) IERSL = 1
 101       IF (IFLAG .LT. 0) IERSL = -1
 102       RETURN
 103 C-----------------------------------------------------------------------
 104 C Use DPCG to solve the linear system P*x = -f
 105 C-----------------------------------------------------------------------
 106  300  CONTINUE
 107       LR = 1
 108       LP = LR + N
 109       LW = LP + N
 110       LZ = LW + N
 111       LWK = LZ + N
 112       CALL DCOPY (N, X, 1, WM(LR), 1)
 113       CALL DPCG (NEQ, TN, Y, SAVF, WM(LR), EWT, N, MAXL, DELTA, HL0,
 114      1          JPRE, MNEWT, F, PSOL, NPSL, X, WM(LP), WM(LW), WM(LZ),
 115      2          LPCG, WM(LOCWP), IWM(LOCIWP), WM(LWK), IFLAG)
 116       NNI = NNI + 1
 117       NLI = NLI + LPCG
 118       NPS = NPS + NPSL
 119       IF (IFLAG .NE. 0) NCFL = NCFL + 1
 120       IF (IFLAG .GE. 2) IERSL = 1
 121       IF (IFLAG .LT. 0) IERSL = -1
 122       RETURN
 123 C-----------------------------------------------------------------------
 124 C Use DPCGS to solve the linear system P*x = -f
 125 C-----------------------------------------------------------------------
 126  400  CONTINUE
 127       LR = 1
 128       LP = LR + N
 129       LW = LP + N
 130       LZ = LW + N
 131       LWK = LZ + N
 132       CALL DCOPY (N, X, 1, WM(LR), 1)
 133       CALL DPCGS (NEQ, TN, Y, SAVF, WM(LR), EWT, N, MAXL, DELTA, HL0,
 134      1           JPRE, MNEWT, F, PSOL, NPSL, X, WM(LP), WM(LW), WM(LZ),
 135      2           LPCG, WM(LOCWP), IWM(LOCIWP), WM(LWK), IFLAG)
 136       NNI = NNI + 1
 137       NLI = NLI + LPCG
 138       NPS = NPS + NPSL
 139       IF (IFLAG .NE. 0) NCFL = NCFL + 1
 140       IF (IFLAG .GE. 2) IERSL = 1
 141       IF (IFLAG .LT. 0) IERSL = -1
 142       RETURN
 143 C-----------------------------------------------------------------------
 144 C Use DUSOL, which interfaces to PSOL, to solve the linear system
 145 C (no Krylov iteration).
 146 C-----------------------------------------------------------------------
 147  900  CONTINUE
 148       LB = 1
 149       LWK = LB + N
 150       CALL DCOPY (N, X, 1, WM(LB), 1)
 151       CALL DUSOL (NEQ, TN, Y, SAVF, WM(LB), EWT, N, DELTA, HL0, MNEWT,
 152      1   PSOL, NPSL, X, WM(LOCWP), IWM(LOCIWP), WM(LWK), IFLAG)
 153       NNI = NNI + 1
 154       NPS = NPS + NPSL
 155       IF (IFLAG .NE. 0) NCFL = NCFL + 1
 156       IF (IFLAG .EQ. 3) IERSL = 1
 157       IF (IFLAG .LT. 0) IERSL = -1
 158       RETURN
 159 C----------------------- End of Subroutine DSOLPK ----------------------
 160       END