relibc_openlibm/slatec/dlsod.f

*DECK DLSOD
      SUBROUTINE DLSOD (DF, NEQ, T, Y, TOUT, RTOL, ATOL, IDID, YPOUT,
     +   YH, YH1, EWT, SAVF, ACOR, WM, IWM, DJAC, INTOUT, TSTOP, TOLFAC,
     +   DELSGN, RPAR, IPAR)
C***BEGIN PROLOGUE  DLSOD
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DDEBDF
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (LSOD-S, DLSOD-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C   DDEBDF  merely allocates storage for  DLSOD  to relieve the user of
C   the inconvenience of a long call list.  Consequently  DLSOD  is used
C   as described in the comments for  DDEBDF .
C
C***SEE ALSO  DDEBDF
C***ROUTINES CALLED  D1MACH, DHSTRT, DINTYD, DSTOD, DVNRMS, XERMSG
C***COMMON BLOCKS    DDEBD1
C***REVISION HISTORY  (YYMMDD)
C   820301  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  DLSOD
C
      INTEGER IBAND, IBEGIN, IDID, IER, IINTEG, IJAC, INIT, INTFLG,
     1      IOWNS, IPAR, IQUIT, ITOL, ITSTOP, IWM, JSTART, K, KFLAG,
     2      KSTEPS, L, LACOR, LDUM, LEWT, LSAVF, LTOL, LWM, LYH, MAXNUM,
     3      MAXORD, METH, MITER, N, NATOLP, NEQ, NFE, NJE, NQ, NQU,
     4      NRTOLP, NST
      DOUBLE PRECISION ABSDEL, ACOR, ATOL, BIG, D1MACH, DEL,
     1      DELSGN, DT, DVNRMS, EL0, EWT,
     2      H, HA, HMIN, HMXI, HU, ROWNS, RPAR, RTOL, SAVF, T, TOL,
     3      TOLD, TOLFAC, TOUT, TSTOP, U, WM, X, Y, YH, YH1, YPOUT
      LOGICAL INTOUT
      CHARACTER*8 XERN1
      CHARACTER*16 XERN3, XERN4
C
      DIMENSION Y(*),YPOUT(*),YH(NEQ,6),YH1(*),EWT(*),SAVF(*),
     1          ACOR(*),WM(*),IWM(*),RTOL(*),ATOL(*),RPAR(*),IPAR(*)
C
C
      COMMON /DDEBD1/ TOLD,ROWNS(210),EL0,H,HMIN,HMXI,HU,X,U,IQUIT,INIT,
     1                LYH,LEWT,LACOR,LSAVF,LWM,KSTEPS,IBEGIN,ITOL,
     2                IINTEG,ITSTOP,IJAC,IBAND,IOWNS(6),IER,JSTART,
     3                KFLAG,LDUM,METH,MITER,MAXORD,N,NQ,NST,NFE,NJE,NQU
C
      EXTERNAL DF, DJAC
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
C       COUNTER IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE
C       EXCESSIVE WORK.
      SAVE MAXNUM
C
      DATA MAXNUM /500/
C
C     ..................................................................
C
C***FIRST EXECUTABLE STATEMENT  DLSOD
      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 D1MACH. THE USER MUST MAKE SURE THAT THE
C        VALUES SET IN D1MACH ARE RELEVANT TO THE COMPUTER BEING USED.
C
         U = D1MACH(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 DSTOD 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.0D0
         NQ = 1
         H = 1.0D0
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', 'DLSOD',
     *      'IN DDEBDF, 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', 'DLSOD',
     *            'IN DDEBDF, 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', 'DLSOD',
     *         'IN DDEBDF, 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.0D0,TOUT-T) .NE. SIGN(1.0D0,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', 'DLSOD',
     *         'IN DDEBDF, 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', 'DLSOD',
     *         'IN DDEBDF, 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', 'DLSOD',
     *         'IN DDEBDF, 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.0D0) THEN
               WRITE (XERN3, '(1PE15.6)') TOUT
               CALL XERMSG ('SLATEC', 'DLSOD',
     *            'IN DDEBDF, 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', 'DLSOD',
     *         'IN DDEBDF, 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
C             AS ASKING FOR THE MOST ACCURATE SOLUTION POSSIBLE. IN THIS
C             CASE, THE RELATIVE ERROR TOLERANCE RTOL IS RESET TO THE
C             SMALLEST VALUE 100*U WHICH IS LIKELY TO BE REASONABLE FOR
C             THIS METHOD AND MACHINE
C
      DO 180 K = 1, NEQ
         IF (RTOL(K) + ATOL(K) .GT. 0.0D0) GO TO 170
            RTOL(K) = 100.0D0*U
            IDID = -2
  170    CONTINUE
C     ...EXIT
         IF (ITOL .EQ. 0) GO TO 190
  180 CONTINUE
  190 CONTINUE
C
      IF (IDID .NE. (-2)) GO TO 200
C        RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A
C                                 SMALL POSITIVE VALUE
         IBEGIN = -1
      GO TO 460
  200 CONTINUE
C        BEGIN BLOCK PERMITTING ...EXITS TO 450
C           BEGIN BLOCK PERMITTING ...EXITS TO 430
C              BEGIN BLOCK PERMITTING ...EXITS TO 260
C                 BEGIN BLOCK PERMITTING ...EXITS TO 230
C
C                    BRANCH ON STATUS OF INITIALIZATION INDICATOR
C                           INIT=0 MEANS INITIAL DERIVATIVES AND
C                           NOMINAL STEP SIZE
C                                  AND DIRECTION NOT YET SET
C                           INIT=1 MEANS NOMINAL STEP SIZE AND
C                           DIRECTION NOT YET SET INIT=2 MEANS NO
C                           FURTHER INITIALIZATION REQUIRED
C
                     IF (INIT .EQ. 0) GO TO 210
C                 ......EXIT
                        IF (INIT .EQ. 1) GO TO 230
C              .........EXIT
                        GO TO 260
  210                CONTINUE
C
C                    ................................................
C
C                         MORE INITIALIZATION --
C                                             -- EVALUATE INITIAL
C                                             DERIVATIVES
C
                     INIT = 1
                     CALL DF(T,Y,YH(1,2),RPAR,IPAR)
                     NFE = 1
C                 ...EXIT
                     IF (T .NE. TOUT) GO TO 230
                     IDID = 2
                     DO 220 L = 1, NEQ
                        YPOUT(L) = YH(L,2)
  220                CONTINUE
                     TOLD = T
C        ............EXIT
                     GO TO 450
  230             CONTINUE
C
C                 -- COMPUTE INITIAL STEP SIZE
C                 -- SAVE SIGN OF INTEGRATION DIRECTION
C                 -- SET INDEPENDENT AND DEPENDENT VARIABLES
C                                      X AND YH(*) FOR DSTOD
C
                  LTOL = 1
                  DO 240 L = 1, NEQ
                     IF (ITOL .EQ. 1) LTOL = L
                     TOL = RTOL(LTOL)*ABS(Y(L)) + ATOL(LTOL)
                     IF (TOL .EQ. 0.0D0) GO TO 390
                     EWT(L) = TOL
  240             CONTINUE
C
                  BIG = SQRT(D1MACH(2))
                  CALL DHSTRT(DF,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,
     2                        IPAR,H)
C
                  DELSGN = SIGN(1.0D0,TOUT-T)
                  X = T
                  DO 250 L = 1, NEQ
                     YH(L,1) = Y(L)
                     YH(L,2) = H*YH(L,2)
  250             CONTINUE
                  INIT = 2
  260          CONTINUE
C
C              ......................................................
C
C                 ON EACH CALL SET INFORMATION WHICH DETERMINES THE
C                 ALLOWED INTERVAL OF INTEGRATION BEFORE RETURNING
C                 WITH AN ANSWER AT TOUT
C
               DEL = TOUT - T
               ABSDEL = ABS(DEL)
C
C              ......................................................
C
C                 IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND
C                 RETURN
C
  270          CONTINUE
C                 BEGIN BLOCK PERMITTING ...EXITS TO 400
C                    BEGIN BLOCK PERMITTING ...EXITS TO 380
                        IF (ABS(X-T) .LT. ABSDEL) GO TO 290
                           CALL DINTYD(TOUT,0,YH,NEQ,Y,INTFLG)
                           CALL DINTYD(TOUT,1,YH,NEQ,YPOUT,INTFLG)
                           IDID = 3
                           IF (X .NE. TOUT) GO TO 280
                              IDID = 2
                              INTOUT = .FALSE.
  280                      CONTINUE
                           T = TOUT
                           TOLD = T
C        ..................EXIT
                           GO TO 450
  290                   CONTINUE
C
C                       IF CANNOT GO PAST TSTOP AND SUFFICIENTLY
C                       CLOSE, EXTRAPOLATE AND RETURN
C
                        IF (ITSTOP .NE. 1) GO TO 310
                        IF (ABS(TSTOP-X) .GE. 100.0D0*U*ABS(X))
     1                     GO TO 310
                           DT = TOUT - X
                           DO 300 L = 1, NEQ
                              Y(L) = YH(L,1) + (DT/H)*YH(L,2)
  300                      CONTINUE
                           CALL DF(TOUT,Y,YPOUT,RPAR,IPAR)
                           NFE = NFE + 1
                           IDID = 3
                           T = TOUT
                           TOLD = T
C        ..................EXIT
                           GO TO 450
  310                   CONTINUE
C
                        IF (IINTEG .EQ. 0 .OR. .NOT.INTOUT) GO TO 320
C
C                          INTERMEDIATE-OUTPUT MODE
C
                           IDID = 1
                        GO TO 370
  320                   CONTINUE
C
C                       .............................................
C
C                            MONITOR NUMBER OF STEPS ATTEMPTED
C
                        IF (KSTEPS .LE. MAXNUM) GO TO 330
C
C                          A SIGNIFICANT AMOUNT OF WORK HAS BEEN
C                          EXPENDED
                           IDID = -1
                           KSTEPS = 0
                           IBEGIN = -1
                        GO TO 370
  330                   CONTINUE
C
C                          ..........................................
C
C                             LIMIT STEP SIZE AND SET WEIGHT VECTOR
C
                           HMIN = 100.0D0*U*ABS(X)
                           HA = MAX(ABS(H),HMIN)
                           IF (ITSTOP .EQ. 1)
     1                        HA = MIN(HA,ABS(TSTOP-X))
                           H = SIGN(HA,H)
                           LTOL = 1
                           DO 340 L = 1, NEQ
                              IF (ITOL .EQ. 1) LTOL = L
                              EWT(L) = RTOL(LTOL)*ABS(YH(L,1))
     1                                 + ATOL(LTOL)
C                    .........EXIT
                              IF (EWT(L) .LE. 0.0D0) GO TO 380
  340                      CONTINUE
                           TOLFAC = U*DVNRMS(NEQ,YH,EWT)
C                 .........EXIT
                           IF (TOLFAC .LE. 1.0D0) GO TO 400
C
C                          TOLERANCES TOO SMALL
                           IDID = -2
                           TOLFAC = 2.0D0*TOLFAC
                           RTOL(1) = TOLFAC*RTOL(1)
                           ATOL(1) = TOLFAC*ATOL(1)
                           IF (ITOL .EQ. 0) GO TO 360
                              DO 350 L = 2, NEQ
                                 RTOL(L) = TOLFAC*RTOL(L)
                                 ATOL(L) = TOLFAC*ATOL(L)
  350                         CONTINUE
  360                      CONTINUE
                           IBEGIN = -1
  370                   CONTINUE
C           ............EXIT
                        GO TO 430
  380                CONTINUE
C
C                    RELATIVE ERROR CRITERION INAPPROPRIATE
  390                CONTINUE
                     IDID = -3
                     IBEGIN = -1
C           .........EXIT
                     GO TO 430
  400             CONTINUE
C
C                 ...................................................
C
C                      TAKE A STEP
C
                  CALL DSTOD(NEQ,Y,YH,NEQ,YH1,EWT,SAVF,ACOR,WM,IWM,
     1                       DF,DJAC,RPAR,IPAR)
C
                  JSTART = -2
                  INTOUT = .TRUE.
               IF (KFLAG .EQ. 0) GO TO 270
C
C              ......................................................
C
               IF (KFLAG .EQ. -1) GO TO 410
C
C                 REPEATED CORRECTOR CONVERGENCE FAILURES
                  IDID = -6
                  IBEGIN = -1
               GO TO 420
  410          CONTINUE
C
C                 REPEATED ERROR TEST FAILURES
                  IDID = -7
                  IBEGIN = -1
  420          CONTINUE
  430       CONTINUE
C
C           .........................................................
C
C                                  STORE VALUES BEFORE RETURNING TO
C                                  DDEBDF
            DO 440 L = 1, NEQ
               Y(L) = YH(L,1)
               YPOUT(L) = YH(L,2)/H
  440       CONTINUE
            T = X
            TOLD = T
            INTOUT = .FALSE.
  450    CONTINUE
  460 CONTINUE
      RETURN
      END