share/odepack/fortran/dsolsy.f

   1 *DECK DSOLSY
   2       SUBROUTINE DSOLSY (WM, IWM, X, TEM)
   3 C***BEGIN PROLOGUE  DSOLSY
   4 C***SUBSIDIARY
   5 C***PURPOSE  ODEPACK linear system solver.
   6 C***TYPE      DOUBLE PRECISION (SSOLSY-S, DSOLSY-D)
   7 C***AUTHOR  Hindmarsh, Alan C., (LLNL)
   8 C***DESCRIPTION
   9 C
  10 C  This routine manages the solution of the linear system arising from
  11 C  a chord iteration.  It is called if MITER .ne. 0.
  12 C  If MITER is 1 or 2, it calls DGESL to accomplish this.
  13 C  If MITER = 3 it updates the coefficient h*EL0 in the diagonal
  14 C  matrix, and then computes the solution.
  15 C  If MITER is 4 or 5, it calls DGBSL.
  16 C  Communication with DSOLSY uses the following variables:
  17 C  WM    = real work space containing the inverse diagonal matrix if
  18 C          MITER = 3 and the LU decomposition of the matrix otherwise.
  19 C          Storage of matrix elements starts at WM(3).
  20 C          WM also contains the following matrix-related data:
  21 C          WM(1) = SQRT(UROUND) (not used here),
  22 C          WM(2) = HL0, the previous value of h*EL0, used if MITER = 3.
  23 C  IWM   = integer work space containing pivot information, starting at
  24 C          IWM(21), if MITER is 1, 2, 4, or 5.  IWM also contains band
  25 C          parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
  26 C  X     = the right-hand side vector on input, and the solution vector
  27 C          on output, of length N.
  28 C  TEM   = vector of work space of length N, not used in this version.
  29 C  IERSL = output flag (in COMMON).  IERSL = 0 if no trouble occurred.
  30 C          IERSL = 1 if a singular matrix arose with MITER = 3.
  31 C  This routine also uses the COMMON variables EL0, H, MITER, and N.
  32 C
  33 C***SEE ALSO  DLSODE
  34 C***ROUTINES CALLED  DGBSL, DGESL
  35 C***COMMON BLOCKS    DLS001
  36 C***REVISION HISTORY  (YYMMDD)
  37 C   791129  DATE WRITTEN
  38 C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
  39 C   890503  Minor cosmetic changes.  (FNF)
  40 C   930809  Renamed to allow single/double precision versions. (ACH)
  41 C   010418  Reduced size of Common block /DLS001/. (ACH)
  42 C   031105  Restored 'own' variables to Common block /DLS001/, to
  43 C           enable interrupt/restart feature. (ACH)
  44 C***END PROLOGUE  DSOLSY
  45 C**End
  46       INTEGER IWM
  47       DOUBLE PRECISION WM, X, TEM
  48       DIMENSION WM(*), IWM(*), X(*), TEM(*)
  49       INTEGER IOWND, IOWNS,
  50      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  51      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  52      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  53       DOUBLE PRECISION ROWNS,
  54      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  55       COMMON /DLS001/ ROWNS(209),
  56      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  57      2   IOWND(6), IOWNS(6),
  58      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  59      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  60      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  61       INTEGER I, MEBAND, ML, MU
  62       DOUBLE PRECISION DI, HL0, PHL0, R
  63 C
  64 C***FIRST EXECUTABLE STATEMENT  DSOLSY
  65       IERSL = 0
  66       GO TO (100, 100, 300, 400, 400), MITER
  67  100  CALL DGESL (WM(3), N, N, IWM(21), X, 0)
  68       RETURN
  69 C
  70  300  PHL0 = WM(2)
  71       HL0 = H*EL0
  72       WM(2) = HL0
  73       IF (HL0 .EQ. PHL0) GO TO 330
  74       R = HL0/PHL0
  75       DO 320 I = 1,N
  76         DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2))
  77         IF (ABS(DI) .EQ. 0.0D0) GO TO 390
  78  320    WM(I+2) = 1.0D0/DI
  79  330  DO 340 I = 1,N
  80  340    X(I) = WM(I+2)*X(I)
  81       RETURN
  82  390  IERSL = 1
  83       RETURN
  84 C
  85  400  ML = IWM(1)
  86       MU = IWM(2)
  87       MEBAND = 2*ML + MU + 1
  88       CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0)
  89       RETURN
  90 C----------------------- END OF SUBROUTINE DSOLSY ----------------------
  91       END