relibc_openlibm/slatec/dslvs.f

*DECK DSLVS
      SUBROUTINE DSLVS (WM, IWM, X, TEM)
C***BEGIN PROLOGUE  DSLVS
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DDEBDF
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (SLVS-S, DSLVS-D)
C***AUTHOR  Watts, H. A., (SNLA)
C***DESCRIPTION
C
C   DSLVS solves the linear system in the iteration scheme for the
C   integrator package DDEBDF.
C
C***SEE ALSO  DDEBDF
C***ROUTINES CALLED  DGBSL, DGESL
C***COMMON BLOCKS    DDEBD1
C***REVISION HISTORY  (YYMMDD)
C   820301  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910722  Updated AUTHOR section.  (ALS)
C   920422  Changed DIMENSION statement.  (WRB)
C***END PROLOGUE  DSLVS
C
      INTEGER I, IER, IOWND, IOWNS, IWM, JSTART, KFLAG, L, MAXORD,
     1      MEBAND, METH, MITER, ML, MU, N, NFE, NJE, NQ, NQU, NST
      DOUBLE PRECISION DI, EL0, H, HL0, HMIN, HMXI, HU, PHL0,
     1      R, ROWND, ROWNS, TEM, TN, UROUND, WM, X
      DIMENSION WM(*), IWM(*), X(*), TEM(*)
      COMMON /DDEBD1/ ROWND,ROWNS(210),EL0,H,HMIN,HMXI,HU,TN,UROUND,
     1                IOWND(14),IOWNS(6),IER,JSTART,KFLAG,L,METH,MITER,
     2                MAXORD,N,NQ,NST,NFE,NJE,NQU
C     ------------------------------------------------------------------
C      THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING
C      FROM A CHORD ITERATION.  IT IS CALLED BY DSTOD  IF MITER .NE. 0.
C      IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS.
C      IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
C      MATRIX, AND THEN COMPUTES THE SOLUTION.
C      IF MITER IS 4 OR 5, IT CALLS DGBSL.
C      COMMUNICATION WITH DSLVS USES THE FOLLOWING VARIABLES..
C      WM  = DOUBLE PRECISION WORK SPACE CONTAINING THE INVERSE DIAGONAL
C      MATRIX IF MITER
C            IS 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE.
C            STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
C            WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
C            WM(1) = SQRT(UROUND) (NOT USED HERE),
C            WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER =
C            3.
C      IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING
C            AT IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS
C            THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS
C            4 OR 5.
C      X   = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION
C            VECTOR ON OUTPUT, OF LENGTH N.
C      TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION.
C      IER = OUTPUT FLAG (IN COMMON).  IER = 0 IF NO TROUBLE OCCURRED.
C            IER = -1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
C      THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
C-----------------------------------------------------------------------
C     BEGIN BLOCK PERMITTING ...EXITS TO 80
C        BEGIN BLOCK PERMITTING ...EXITS TO 60
C***FIRST EXECUTABLE STATEMENT  DSLVS
            IER = 0
            GO TO (10,10,20,70,70), MITER
   10       CONTINUE
            CALL DGESL(WM(3),N,N,IWM(21),X,0)
C     ......EXIT
            GO TO 80
C
   20       CONTINUE
            PHL0 = WM(2)
            HL0 = H*EL0
            WM(2) = HL0
            IF (HL0 .EQ. PHL0) GO TO 40
               R = HL0/PHL0
               DO 30 I = 1, N
                  DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2))
C        .........EXIT
                  IF (ABS(DI) .EQ. 0.0D0) GO TO 60
                  WM(I+2) = 1.0D0/DI
   30          CONTINUE
   40       CONTINUE
            DO 50 I = 1, N
               X(I) = WM(I+2)*X(I)
   50       CONTINUE
C     ......EXIT
            GO TO 80
   60    CONTINUE
         IER = -1
C     ...EXIT
         GO TO 80
C
   70    CONTINUE
         ML = IWM(1)
         MU = IWM(2)
         MEBAND = 2*ML + MU + 1
         CALL DGBSL(WM(3),MEBAND,N,ML,MU,IWM(21),X,0)
   80 CONTINUE
      RETURN
C     ----------------------- END OF SUBROUTINE DSLVS
C     -----------------------
      END