share/odepack/fortran/dsolss.f

   1 *DECK DSOLSS
   2       SUBROUTINE DSOLSS (WK, IWK, X, TEM)
   3       INTEGER IWK
   4       DOUBLE PRECISION WK, X, TEM
   5       DIMENSION WK(*), IWK(*), X(*), TEM(*)
   6       INTEGER IOWND, IOWNS,
   7      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
   8      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
   9      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  10       INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
  11      1   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
  12      2   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
  13      3   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
  14       DOUBLE PRECISION ROWNS,
  15      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  16       DOUBLE PRECISION RLSS
  17       COMMON /DLS001/ ROWNS(209),
  18      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  19      2   IOWND(6), IOWNS(6),
  20      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  21      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  22      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  23       COMMON /DLSS01/ RLSS(6),
  24      1   IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
  25      2   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
  26      3   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
  27      4   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
  28       INTEGER I
  29       DOUBLE PRECISION DI, HL0, PHL0, R
  30 C-----------------------------------------------------------------------
  31 C This routine manages the solution of the linear system arising from
  32 C a chord iteration.  It is called if MITER .ne. 0.
  33 C If MITER is 1 or 2, it calls CDRV to accomplish this.
  34 C If MITER = 3 it updates the coefficient H*EL0 in the diagonal
  35 C matrix, and then computes the solution.
  36 C communication with DSOLSS uses the following variables:
  37 C WK    = real work space containing the inverse diagonal matrix if
  38 C         MITER = 3 and the LU decomposition of the matrix otherwise.
  39 C         Storage of matrix elements starts at WK(3).
  40 C         WK also contains the following matrix-related data:
  41 C         WK(1) = SQRT(UROUND) (not used here),
  42 C         WK(2) = HL0, the previous value of H*EL0, used if MITER = 3.
  43 C IWK   = integer work space for matrix-related data, assumed to
  44 C         be equivalenced to WK.  In addition, WK(IPRSP) and IWK(IPISP)
  45 C         are assumed to have identical locations.
  46 C X     = the right-hand side vector on input, and the solution vector
  47 C         on output, of length N.
  48 C TEM   = vector of work space of length N, not used in this version.
  49 C IERSL = output flag (in Common).
  50 C         IERSL = 0  if no trouble occurred.
  51 C         IERSL = -1 if CDRV returned an error flag (MITER = 1 or 2).
  52 C                    This should never occur and is considered fatal.
  53 C         IERSL = 1  if a singular matrix arose with MITER = 3.
  54 C This routine also uses other variables in Common.
  55 C-----------------------------------------------------------------------
  56       IERSL = 0
  57       GO TO (100, 100, 300), MITER
  58  100  CALL CDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN),
  59      1   WK(IPA),X,X,NSP,IWK(IPISP),WK(IPRSP),IESP,4,IERSL)
  60       IF (IERSL .NE. 0) IERSL = -1
  61       RETURN
  62 C
  63  300  PHL0 = WK(2)
  64       HL0 = H*EL0
  65       WK(2) = HL0
  66       IF (HL0 .EQ. PHL0) GO TO 330
  67       R = HL0/PHL0
  68       DO 320 I = 1,N
  69         DI = 1.0D0 - R*(1.0D0 - 1.0D0/WK(I+2))
  70         IF (ABS(DI) .EQ. 0.0D0) GO TO 390
  71  320    WK(I+2) = 1.0D0/DI
  72  330  DO 340 I = 1,N
  73  340    X(I) = WK(I+2)*X(I)
  74       RETURN
  75  390  IERSL = 1
  76       RETURN
  77 C
  78 C----------------------- End of Subroutine DSOLSS ----------------------
  79       END