share/odepack/fortran/dainvg.f

   1 *DECK DAINVG
   2       SUBROUTINE DAINVG (RES, ADDA, NEQ, T, Y, YDOT, MITER,
   3      1                   ML, MU, PW, IPVT, IER )
   4       EXTERNAL RES, ADDA
   5       INTEGER NEQ, MITER, ML, MU, IPVT, IER
   6       INTEGER I, LENPW, MLP1, NROWPW
   7       DOUBLE PRECISION T, Y, YDOT, PW
   8       DIMENSION Y(*), YDOT(*), PW(*), IPVT(*)
   9       dimension neq(1)
  10 C-----------------------------------------------------------------------
  11 C This subroutine computes the initial value
  12 C of the vector YDOT satisfying
  13 C     A * YDOT = g(t,y)
  14 C when A is nonsingular.  It is called by DLSODI for
  15 C initialization only, when ISTATE = 0 .
  16 C DAINVG returns an error flag IER:
  17 C   IER  =  0  means DAINVG was successful.
  18 C   IER .ge. 2 means RES returned an error flag IRES = IER.
  19 C   IER .lt. 0 means the a-matrix was found to be singular.
  20 C-----------------------------------------------------------------------
  21 C
  22       IF (MITER .GE. 4)  GO TO 100
  23 C
  24 C Full matrix case -----------------------------------------------------
  25 C
  26       LENPW = NEQ(1)*NEQ(1)
  27       DO 10  I = 1, LENPW
  28    10    PW(I) = 0.0D0
  29 C
  30       IER = 1
  31       CALL RES ( NEQ, T, Y, PW, YDOT, IER )
  32       IF (IER .GT. 1) RETURN
  33 C
  34       itemp = neq(1)
  35       CALL ADDA ( NEQ, T, Y, 0, 0, PW, itemp )
  36       CALL DGEFA ( PW, itemp, itemp, IPVT, IER )
  37       IF (IER .EQ. 0) GO TO 20
  38          IER = -IER
  39          RETURN
  40    20 CALL DGESL ( PW, itemp, itemp, IPVT, YDOT, 0 )
  41       RETURN
  42 C
  43 C Band matrix case -----------------------------------------------------
  44 C
  45   100 CONTINUE
  46       NROWPW = 2*ML + MU + 1
  47       LENPW = NEQ(1) * NROWPW
  48       DO 110  I = 1, LENPW
  49   110    PW(I) = 0.0D0
  50 C
  51       IER = 1
  52       CALL RES ( NEQ, T, Y, PW, YDOT, IER )
  53       IF (IER .GT. 1) RETURN
  54 C
  55       MLP1 = ML + 1
  56       CALL ADDA ( NEQ, T, Y, ML, MU, PW(MLP1), NROWPW )
  57       CALL DGBFA ( PW, NROWPW, NEQ, ML, MU, IPVT, IER )
  58       IF (IER .EQ. 0) GO TO 120
  59          IER = -IER
  60          RETURN
  61   120 CALL DGBSL ( PW, NROWPW, NEQ, ML, MU, IPVT, YDOT, 0 )
  62       RETURN
  63 C----------------------- End of Subroutine DAINVG ----------------------
  64       END