relibc_openlibm/slatec/lsod.f

*DECK LSOD
      SUBROUTINE LSOD (F, NEQ, T, Y, TOUT, RTOL, ATOL, IDID, YPOUT, YH,
     +   YH1, EWT, SAVF, ACOR, WM, IWM, JAC, INTOUT, TSTOP, TOLFAC,
     +   DELSGN, RPAR, IPAR)
C***BEGIN PROLOGUE  LSOD
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DEBDF
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (LSOD-S, DLSOD-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C   DEBDF  merely allocates storage for  LSOD  to relieve the user of
C   the inconvenience of a long call list.  Consequently  LSOD  is used
C   as described in the comments for  DEBDF .
C
C***SEE ALSO  DEBDF
C***ROUTINES CALLED  HSTART, INTYD, R1MACH, STOD, VNWRMS, XERMSG
C***COMMON BLOCKS    DEBDF1
C***REVISION HISTORY  (YYMMDD)
C   800901  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
C***END PROLOGUE  LSOD
C
C
      LOGICAL INTOUT
C
      DIMENSION Y(*),YPOUT(*),YH(NEQ,6),YH1(*),EWT(*),SAVF(*),
     1          ACOR(*),WM(*),IWM(*),RTOL(*),ATOL(*),RPAR(*),IPAR(*)
      CHARACTER*8 XERN1
      CHARACTER*16 XERN3, XERN4
C
      COMMON /DEBDF1/ TOLD, ROWNS(210),
     1   EL0, H, HMIN, HMXI, HU, X, U,
     2   IQUIT, INIT, LYH, LEWT, LACOR, LSAVF, LWM, KSTEPS,
     3   IBEGIN, ITOL, IINTEG, ITSTOP, IJAC, IBAND, IOWNS(6),
     4   IER, JSTART, KFLAG, LDUM, METH, MITER, MAXORD, N, NQ, NST,
     5   NFE, NJE, NQU
C
      EXTERNAL F, JAC
C
C.......................................................................
C
C  THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE
C  NUMBER OF  STEPS ATTEMPTED. WHEN THIS EXCEEDS  MAXNUM, THE COUNTER
C  IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE EXCESSIVE
C  WORK.
C
      SAVE MAXNUM
      DATA MAXNUM/500/
C
C.......................................................................
C
C***FIRST EXECUTABLE STATEMENT  LSOD
      IF (IBEGIN .EQ. 0) THEN
C
C        ON THE FIRST CALL , PERFORM INITIALIZATION --
C        DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY  U  BY CALLING THE
C        FUNCTION ROUTINE R1MACH. THE USER MUST MAKE SURE THAT THE
C        VALUES SET IN R1MACH ARE RELEVANT TO THE COMPUTER BEING USED.
C
         U = R1MACH(4)
C                          -- SET ASSOCIATED MACHINE DEPENDENT PARAMETER
         WM(1) = SQRT(U)
C                          -- SET TERMINATION FLAG
         IQUIT = 0
C                          -- SET INITIALIZATION INDICATOR
         INIT = 0
C                          -- SET COUNTER FOR ATTEMPTED STEPS
         KSTEPS = 0
C                          -- SET INDICATOR FOR INTERMEDIATE-OUTPUT
         INTOUT = .FALSE.
C                          -- SET START INDICATOR FOR STOD CODE
         JSTART = 0
C                          -- SET BDF METHOD INDICATOR
         METH = 2
C                          -- SET MAXIMUM ORDER FOR BDF METHOD
         MAXORD = 5
C                          -- SET ITERATION MATRIX INDICATOR
C
         IF (IJAC .EQ. 0 .AND. IBAND .EQ. 0) MITER = 2
         IF (IJAC .EQ. 1 .AND. IBAND .EQ. 0) MITER = 1
         IF (IJAC .EQ. 0 .AND. IBAND .EQ. 1) MITER = 5
         IF (IJAC .EQ. 1 .AND. IBAND .EQ. 1) MITER = 4
C
C                          -- SET OTHER NECESSARY ITEMS IN COMMON BLOCK
         N = NEQ
         NST = 0
         NJE = 0
         HMXI = 0.
         NQ = 1
         H = 1.
C                          -- RESET IBEGIN FOR SUBSEQUENT CALLS
         IBEGIN=1
      ENDIF
C
C.......................................................................
C
C      CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY
C
      IF (NEQ .LT. 1) THEN
         WRITE (XERN1, '(I8)') NEQ
         CALL XERMSG ('SLATEC', 'LSOD',
     *      'IN DEBDF, THE NUMBER OF EQUATIONS MUST BE A POSITIVE ' //
     *      'INTEGER.$$YOU HAVE CALLED THE CODE WITH NEQ = ' // XERN1,
     *      6, 1)
         IDID=-33
      ENDIF
C
      NRTOLP = 0
      NATOLP = 0
      DO 60 K = 1,NEQ
         IF (NRTOLP .LE. 0) THEN
            IF (RTOL(K) .LT. 0.) THEN
               WRITE (XERN1, '(I8)') K
               WRITE (XERN3, '(1PE15.6)') RTOL(K)
               CALL XERMSG ('SLATEC', 'LSOD',
     *            'IN DEBDF, THE RELATIVE ERROR TOLERANCES MUST ' //
     *            'BE NON-NEGATIVE.$$YOU HAVE CALLED THE CODE WITH ' //
     *            'RTOL(' // XERN1 // ') = ' // XERN3 // '$$IN THE ' //
     *            'CASE OF VECTOR ERROR TOLERANCES, NO FURTHER ' //
     *            'CHECKING OF RTOL COMPONENTS IS DONE.', 7, 1)
               IDID = -33
               IF (NATOLP .GT. 0) GO TO 70
               NRTOLP = 1
            ELSEIF (NATOLP .GT. 0) THEN
               GO TO 50
            ENDIF
         ENDIF
C
         IF (ATOL(K) .LT. 0.) THEN
            WRITE (XERN1, '(I8)') K
            WRITE (XERN3, '(1PE15.6)') ATOL(K)
            CALL XERMSG ('SLATEC', 'LSOD',
     *         'IN DEBDF, THE ABSOLUTE ERROR ' //
     *         'TOLERANCES MUST BE NON-NEGATIVE.$$YOU HAVE CALLED ' //
     *         'THE CODE WITH ATOL(' // XERN1 // ') = ' // XERN3 //
     *         '$$IN THE CASE OF VECTOR ERROR TOLERANCES, NO FURTHER '
     *         // 'CHECKING OF ATOL COMPONENTS IS DONE.', 8, 1)
            IDID=-33
            IF (NRTOLP .GT. 0) GO TO 70
            NATOLP=1
         ENDIF
   50    IF (ITOL .EQ. 0) GO TO 70
   60 CONTINUE
C
   70 IF (ITSTOP .EQ. 1) THEN
         IF (SIGN(1.,TOUT-T) .NE. SIGN(1.,TSTOP-T) .OR.
     1      ABS(TOUT-T) .GT. ABS(TSTOP-T)) THEN
            WRITE (XERN3, '(1PE15.6)') TOUT
            WRITE (XERN4, '(1PE15.6)') TSTOP
            CALL XERMSG ('SLATEC', 'LSOD',
     *         'IN DEBDF, YOU HAVE CALLED THE ' //
     *         'CODE WITH TOUT = ' // XERN3 // '$$BUT YOU HAVE ' //
     *         'ALSO TOLD THE CODE NOT TO INTEGRATE PAST THE POINT ' //
     *         'TSTOP = ' // XERN4 // ' BY SETTING INFO(4) = 1.  ' //
     *         'THESE INSTRUCTIONS CONFLICT.', 14, 1)
            IDID=-33
         ENDIF
      ENDIF
C
C        CHECK SOME CONTINUATION POSSIBILITIES
C
      IF (INIT .NE. 0) THEN
         IF (T .EQ. TOUT) THEN
            WRITE (XERN3, '(1PE15.6)') T
            CALL XERMSG ('SLATEC', 'LSOD',
     *         'IN DEBDF, YOU HAVE CALLED THE CODE WITH T = TOUT = ' //
     *         XERN3 // '  THIS IS NOT ALLOWED ON CONTINUATION CALLS.',
     *         9, 1)
            IDID=-33
         ENDIF
C
         IF (T .NE. TOLD) THEN
            WRITE (XERN3, '(1PE15.6)') TOLD
            WRITE (XERN4, '(1PE15.6)') T
            CALL XERMSG ('SLATEC', 'LSOD',
     *         'IN DEBDF, YOU HAVE CHANGED THE VALUE OF T FROM ' //
     *         XERN3 // ' TO ' // XERN4 //
     *         '  THIS IS NOT ALLOWED ON CONTINUATION CALLS.', 10, 1)
            IDID=-33
         ENDIF
C
         IF (INIT .NE. 1) THEN
            IF (DELSGN*(TOUT-T) .LT. 0.) THEN
               WRITE (XERN3, '(1PE15.6)') TOUT
               CALL XERMSG ('SLATEC', 'LSOD',
     *            'IN DEBDF, BY CALLING THE CODE WITH TOUT = ' //
     *            XERN3 // ' YOU ARE ATTEMPTING TO CHANGE THE ' //
     *            'DIRECTION OF INTEGRATION.$$' //
     *            'THIS IS NOT ALLOWED WITHOUT RESTARTING.', 11, 1)
               IDID=-33
            ENDIF
         ENDIF
      ENDIF
C
      IF (IDID .EQ. (-33)) THEN
         IF (IQUIT .NE. (-33)) THEN
C                       INVALID INPUT DETECTED
            IQUIT=-33
            IBEGIN=-1
         ELSE
            CALL XERMSG ('SLATEC', 'LSOD',
     *         'IN DEBDF, INVALID INPUT WAS ' //
     *         'DETECTED ON SUCCESSIVE ENTRIES.  IT IS IMPOSSIBLE ' //
     *         'TO PROCEED BECAUSE YOU HAVE NOT CORRECTED THE ' //
     *         'PROBLEM, SO EXECUTION IS BEING TERMINATED.', 12, 2)
         ENDIF
         RETURN
      ENDIF
C
C.......................................................................
C
C     RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND INTERPRETED AS
C     ASKING FOR THE MOST ACCURATE SOLUTION POSSIBLE. IN THIS CASE,
C     THE RELATIVE ERROR TOLERANCE RTOL IS RESET TO THE SMALLEST VALUE
C     100*U WHICH IS LIKELY TO BE REASONABLE FOR THIS METHOD AND MACHINE
C
      DO 170 K=1,NEQ
        IF (RTOL(K)+ATOL(K) .GT. 0.) GO TO 160
        RTOL(K)=100.*U
        IDID=-2
  160   IF (ITOL .EQ. 0) GO TO 180
  170   CONTINUE
C
  180 IF (IDID .NE. (-2)) GO TO 190
C                       RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A
C                                                SMALL POSITIVE VALUE
      IBEGIN=-1
      RETURN
C
C     BRANCH ON STATUS OF INITIALIZATION INDICATOR
C            INIT=0 MEANS INITIAL DERIVATIVES AND NOMINAL STEP SIZE
C                   AND DIRECTION NOT YET SET
C            INIT=1 MEANS NOMINAL STEP SIZE AND DIRECTION NOT YET SET
C            INIT=2 MEANS NO FURTHER INITIALIZATION REQUIRED
C
  190 IF (INIT .EQ. 0) GO TO 200
      IF (INIT .EQ. 1) GO TO 220
      GO TO 240
C
C.......................................................................
C
C     MORE INITIALIZATION --
C                         -- EVALUATE INITIAL DERIVATIVES
C
  200 INIT=1
      CALL F(T,Y,YH(1,2),RPAR,IPAR)
      NFE=1
      IF (T .NE. TOUT) GO TO 220
      IDID=2
      DO 210 L = 1,NEQ
  210    YPOUT(L) = YH(L,2)
      TOLD=T
      RETURN
C
C                         -- COMPUTE INITIAL STEP SIZE
C                         -- SAVE SIGN OF INTEGRATION DIRECTION
C                         -- SET INDEPENDENT AND DEPENDENT VARIABLES
C                                              X AND YH(*) FOR STOD
C
  220 LTOL = 1
      DO 225 L=1,NEQ
        IF (ITOL .EQ. 1) LTOL = L
        TOL = RTOL(LTOL)*ABS(Y(L)) + ATOL(LTOL)
        IF (TOL .EQ. 0.) GO TO 380
  225   EWT(L) = TOL
C
      BIG = SQRT(R1MACH(2))
      CALL HSTART (F,NEQ,T,TOUT,Y,YH(1,2),EWT,1,U,BIG,
     1             YH(1,3),YH(1,4),YH(1,5),YH(1,6),RPAR,IPAR,H)
C
      DELSGN = SIGN(1.0,TOUT-T)
      X = T
      DO 230 L = 1,NEQ
        YH(L,1) = Y(L)
  230   YH(L,2) = H*YH(L,2)
      INIT = 2
C
C.......................................................................
C
C   ON EACH CALL SET INFORMATION WHICH DETERMINES THE ALLOWED INTERVAL
C   OF INTEGRATION BEFORE RETURNING WITH AN ANSWER AT TOUT
C
  240 DEL = TOUT - T
      ABSDEL = ABS(DEL)
C
C.......................................................................
C
C   IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN
C
  250 IF (ABS(X-T) .LT. ABSDEL) GO TO 270
      CALL INTYD(TOUT,0,YH,NEQ,Y,INTFLG)
      CALL INTYD(TOUT,1,YH,NEQ,YPOUT,INTFLG)
      IDID = 3
      IF (X .NE. TOUT) GO TO 260
      IDID = 2
      INTOUT = .FALSE.
  260 T = TOUT
      TOLD = T
      RETURN
C
C   IF CANNOT GO PAST TSTOP AND SUFFICIENTLY CLOSE,
C   EXTRAPOLATE AND RETURN
C
  270 IF (ITSTOP .NE. 1) GO TO 290
      IF (ABS(TSTOP-X) .GE. 100.*U*ABS(X)) GO TO 290
      DT = TOUT - X
      DO 280 L = 1,NEQ
  280   Y(L) = YH(L,1) + (DT/H)*YH(L,2)
      CALL F(TOUT,Y,YPOUT,RPAR,IPAR)
      NFE = NFE + 1
      IDID = 3
      T = TOUT
      TOLD = T
      RETURN
C
  290 IF (IINTEG .EQ. 0  .OR.  .NOT.INTOUT) GO TO 300
C
C   INTERMEDIATE-OUTPUT MODE
C
      IDID = 1
      GO TO 500
C
C.......................................................................
C
C     MONITOR NUMBER OF STEPS ATTEMPTED
C
  300 IF (KSTEPS .LE. MAXNUM) GO TO 330
C
C                       A SIGNIFICANT AMOUNT OF WORK HAS BEEN EXPENDED
      IDID=-1
      KSTEPS=0
      IBEGIN = -1
      GO TO 500
C
C.......................................................................
C
C   LIMIT STEP SIZE AND SET WEIGHT VECTOR
C
  330 HMIN = 100.*U*ABS(X)
      HA = MAX(ABS(H),HMIN)
      IF (ITSTOP .NE. 1) GO TO 340
      HA = MIN(HA,ABS(TSTOP-X))
  340 H = SIGN(HA,H)
      LTOL = 1
      DO 350 L = 1,NEQ
        IF (ITOL .EQ. 1) LTOL = L
        EWT(L) = RTOL(LTOL)*ABS(YH(L,1)) + ATOL(LTOL)
        IF (EWT(L) .LE. 0.0) GO TO 380
  350   CONTINUE
      TOLFAC = U*VNWRMS(NEQ,YH,EWT)
      IF (TOLFAC .LE. 1.) GO TO 400
C
C                       TOLERANCES TOO SMALL
      IDID = -2
      TOLFAC = 2.*TOLFAC
      RTOL(1) = TOLFAC*RTOL(1)
      ATOL(1) = TOLFAC*ATOL(1)
      IF (ITOL .EQ. 0) GO TO 370
      DO 360 L = 2,NEQ
        RTOL(L) = TOLFAC*RTOL(L)
  360   ATOL(L) = TOLFAC*ATOL(L)
  370 IBEGIN = -1
      GO TO 500
C
C                       RELATIVE ERROR CRITERION INAPPROPRIATE
  380 IDID = -3
      IBEGIN = -1
      GO TO 500
C
C.......................................................................
C
C     TAKE A STEP
C
  400 CALL STOD(NEQ,Y,YH,NEQ,YH1,EWT,SAVF,ACOR,WM,IWM,F,JAC,RPAR,IPAR)
C
      JSTART = -2
      INTOUT = .TRUE.
      IF (KFLAG .EQ. 0) GO TO 250
C
C.......................................................................
C
      IF (KFLAG .EQ. -1) GO TO 450
C
C                       REPEATED CORRECTOR CONVERGENCE FAILURES
      IDID = -6
      IBEGIN = -1
      GO TO 500
C
C                       REPEATED ERROR TEST FAILURES
  450 IDID = -7
      IBEGIN = -1
C
C.......................................................................
C
C                       STORE VALUES BEFORE RETURNING TO DEBDF
  500 DO 555 L = 1,NEQ
        Y(L) = YH(L,1)
  555   YPOUT(L) = YH(L,2)/H
      T = X
      TOLD = T
      INTOUT = .FALSE.
      RETURN
      END