share/odepack/fortran/dainvgs.f

   1 *DECK DAINVGS
   2       SUBROUTINE DAINVGS (NEQ, T, Y, WK, IWK, TEM, YDOT, IER, RES, ADDA)
   3       EXTERNAL RES, ADDA
   4       INTEGER NEQ, IWK, IER
   5       INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
   6      1   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
   7      2   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
   8      3   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
   9       INTEGER I, IMUL, J, K, KMIN, KMAX
  10       DOUBLE PRECISION T, Y, WK, TEM, YDOT
  11       DOUBLE PRECISION RLSS
  12       DIMENSION Y(*), WK(*), IWK(*), TEM(*), YDOT(*)
  13       COMMON /DLSS01/ RLSS(6),
  14      1   IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
  15      2   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
  16      3   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
  17      4   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
  18 C-----------------------------------------------------------------------
  19 C This subroutine computes the initial value of the vector YDOT
  20 C satisfying
  21 C     A * YDOT = g(t,y)
  22 C when A is nonsingular.  It is called by DLSODIS for initialization
  23 C only, when ISTATE = 0.  The matrix A is subjected to LU
  24 C decomposition in CDRV.  Then the system A*YDOT = g(t,y) is solved
  25 C in CDRV.
  26 C In addition to variables described previously, communication
  27 C with DAINVGS uses the following:
  28 C Y     = array of initial values.
  29 C WK    = real work space for matrices.  On output it contains A and
  30 C         its LU decomposition.  The LU decomposition is not entirely
  31 C         sparse unless the structure of the matrix A is identical to
  32 C         the structure of the Jacobian matrix dr/dy.
  33 C         Storage of matrix elements starts at WK(3).
  34 C         WK(1) = SQRT(UROUND), not used here.
  35 C IWK   = integer work space for matrix-related data, assumed to
  36 C         be equivalenced to WK.  In addition, WK(IPRSP) and WK(IPISP)
  37 C         are assumed to have identical locations.
  38 C TEM   = vector of work space of length N (ACOR in DSTODI).
  39 C YDOT  = output vector containing the initial dy/dt. YDOT(i) contains
  40 C         dy(i)/dt when the matrix A is non-singular.
  41 C IER   = output error flag with the following values and meanings:
  42 C       = 0  if DAINVGS was successful.
  43 C       = 1  if the A-matrix was found to be singular.
  44 C       = 2  if RES returned an error flag IRES = IER = 2.
  45 C       = 3  if RES returned an error flag IRES = IER = 3.
  46 C       = 4  if insufficient storage for CDRV (should not occur here).
  47 C       = 5  if other error found in CDRV (should not occur here).
  48 C-----------------------------------------------------------------------
  49 C
  50       DO 10 I = 1,NNZ
  51  10     WK(IBA+I) = 0.0D0
  52 C
  53       IER = 1
  54       CALL RES (NEQ, T, Y, WK(IPA), YDOT, IER)
  55       IF (IER .GT. 1) RETURN
  56 C
  57       KMIN = IWK(IPIAN)
  58       DO 30 J = 1,NEQ
  59         KMAX = IWK(IPIAN+J) - 1
  60         DO 15 K = KMIN,KMAX
  61           I = IWK(IBJAN+K)
  62  15       TEM(I) = 0.0D0
  63         CALL ADDA (NEQ, T, Y, J, IWK(IPIAN), IWK(IPJAN), TEM)
  64         DO 20 K = KMIN,KMAX
  65           I = IWK(IBJAN+K)
  66  20       WK(IBA+K) = TEM(I)
  67         KMIN = KMAX + 1
  68  30   CONTINUE
  69       NLU = NLU + 1
  70       IER = 0
  71       DO 40 I = 1,NEQ
  72  40     TEM(I) = 0.0D0
  73 C
  74 C Numerical factorization of matrix A. ---------------------------------
  75       CALL CDRV (NEQ,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN),
  76      1  WK(IPA),TEM,TEM,NSP,IWK(IPISP),WK(IPRSP),IESP,2,IYS)
  77       IF (IYS .EQ. 0) GO TO 50
  78       IMUL = (IYS - 1)/NEQ
  79       IER = 5
  80       IF (IMUL .EQ. 8) IER = 1
  81       IF (IMUL .EQ. 10) IER = 4
  82       RETURN
  83 C
  84 C Solution of the linear system. ---------------------------------------
  85  50   CALL CDRV (NEQ,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN),
  86      1  WK(IPA),YDOT,YDOT,NSP,IWK(IPISP),WK(IPRSP),IESP,4,IYS)
  87       IF (IYS .NE. 0) IER = 5
  88       RETURN
  89 C----------------------- End of Subroutine DAINVGS ---------------------
  90       END