C
C Main program to execute SDARDM demonstration program.
C
C***BEGIN PROLOGUE  SDARDM
C***DATE WRITTEN   900214     (YYMMDD)
C***REVISION DATE  900214     (YYMMDD)
C***CATEGORY NO.   D2A2
C***KEYWORDS  LIBRARY=SLATEC,QUICK CHECK,SDARDM,
C             TYPE=REAL
C***AUTHOR  PETZOLD,LINDA R.,COMPUTING AND MATHEMATICS RESEARCH DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             L - 316, P.O. BOX 808,
C             LIVERMORE, CA.    94550
C***PURPOSE  Quick check for SLATEC routine SDASRT.
c***DESCRIPTION
C       Demonstration program for SDASRT.
C       This version is in double precision.
C
C       SDASRT is used to solve two simple problems,
C       one nonstiff and one intermittently stiff.
C       If the errors are too large, or other difficulty occurs,
C       a warning message is printed.  All output is on unit LUN.
C
C       To run the demonstration problems with full printing,
C       set KPRINT = 3.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  SDASRT,RES1,GR1,RES2,JAC2,GR2
C***END PROLOGUE SDARDM
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL (A-H,O-Z)
      EXTERNAL RES1, GR1, RES2, JAC2, GR2
      INTEGER IDID, IOUT, IWORK, JROOT, INFO, JTYPE,
     *   KROOT, LENIW, LENRW, LIW, LRW, LUN, NEQ, NERR, NG,
     *   NRE, NREA, NGE, NJE, NST
      REAL ATOL, ER, ERO, ERRT, RTOL, RWORK,
     *   T, TOUT, TZERO, Y, YT
      DIMENSION Y(2), YPRIME(2), RTOL(2), ATOL(2),
     *   RWORK(100), IWORK(100), JROOT(2), INFO(15)
      COMMON /LOCAL/NEQ
C
      LUN = 6
      KPRINT = 3
      IPASS = 1
      NERR = 0
C-----------------------------------------------------------------------
C First problem.
C The initial value problem is..
C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
C The two root functions are..
C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
C-----------------------------------------------------------------------
C SET ALL INPUT PARAMETERS AND PRINT HEADING.
      DO 10 I = 1,15
 10      INFO(I) = 0
      NEQ = 1
      T = 1.0E0
      Y(1) = 1.0E0
      TOUT = 2.0E0
      RTOL(1) = 0.0E0
      ATOL(1) = 1.0E-3
      LRW = 100
      LIW = 100
      IDID = 0
C
C SET INFO(11) TO ONE IF SDASRT IS TO COMPUTE THE INITIAL
C YPRIME, AND GENERATE AN INITIAL GUESS FOR YPRIME.
C OTHERWISE, SET INFO(11) = 0, AND SUPPLY THE CORRECT INITIAL
C VALUE FOR YPRIME.
C
      YPRIME(1) = 3.0E0
      INFO(11) = 0
C
C NOTE: JTYPE INDICATES THE JACOBIAN TYPE:
C       JTYPE = 1 ==> JACOBIAN IS DENSE AND USER-SUPPLIED
C       JTYPE = 2 ==> JACOBIAN IS DENSE AND COMPUTED INTERNALLY
C
      JTYPE = 2
      INFO(5) = 2 - JTYPE
      NG = 2
      IF (KPRINT .GE. 2) THEN
         WRITE (LUN,110) RTOL(1),ATOL(1),JTYPE
      ENDIF
 110  FORMAT('1'/1X,1P,' DEMONSTRATION PROGRAM FOR SDASRT'/////
     1  1X,' PROBLEM 1..'///
     2  1X,' problem is  DY/DT = ((2*LOG(Y)+8)/T - 5)*Y,  Y(1) = 1'//
     3  1X,' solution is  Y(T) = EXP(-T**2 + 5*T - 4)'//
     4  1X,' root functions are..'/
     5  10X,' G1 = DY/DT  (root at T = 2.5).'/
     6  10X,' G2 = LOG(Y) - 6.2491  (roots at T = 2.47 and T = 2.53).'//
     7  1X,' RTOL =',E10.1,'   ATOL =',E10.1/
     8  1X,' JTYPE =',I3/////)
C
C CALL SDASRT IN LOOP OVER TOUT VALUES 2,3,4,5,6
      ERO = 0.0E0
      DO 180 IOUT = 1,5
 120    CONTINUE
C
        CALL SDASRT (RES1,NEQ,T,Y,YPRIME,TOUT,
     *  INFO,RTOL,ATOL,IDID,
     *  RWORK,LRW,IWORK,LIW,RPAR,IPAR,JDUM,
     *  GR1,NG,JROOT)
C
C PRINT Y AND ERROR IN Y, AND PRINT WARNING IF ERROR TOO LARGE.
        YT = EXP(-T*T + 5.0E0*T - 4.0E0)
        ER = Y(1) - YT
        IF (KPRINT .GT. 2) THEN
           WRITE (LUN,130) T,Y(1),ER
        ENDIF
 130    FORMAT(1X,1P,' at t =',E15.7,5X,'y =',E15.7,5X,'error =',E12.4)
        IF (IDID .LT. 0) GO TO 185
        ER = ABS(ER)/ATOL(1)
        ERO = AMAX1(ERO,ER)
        IF (ER .LT. 1000.0E0) GO TO 140
        IF (KPRINT .GE. 2) THEN
          IPASS = 0
          WRITE (LUN,135)
        ENDIF
 135    FORMAT(//1X,' WARNING.. ERROR EXCEEDS 1000 * TOLERANCE'//)
        NERR = NERR + 1
 140    CONTINUE
        IF (IDID .NE. 4) GO TO 175
C
C IF A ROOT WAS FOUND, WRITE RESULTS AND CHECK ROOT LOCATION.
C THEN RETURN TO SDASRT TO CONTINUE THE INTEGRATION.
        IF (KPRINT .GT. 2) THEN
           WRITE (LUN,150) T,JROOT(1),JROOT(2)
        ENDIF
 150    FORMAT(/1X,1P,' root found at t =',E15.7,5X,'jroot =',2I5)
        IF (JROOT(1) .EQ. 1) ERRT = T - 2.5E0
        IF (JROOT(2) .EQ. 1 .AND. T .LE. 2.5E0) ERRT = T - 2.47E0
        IF (JROOT(2) .EQ. 1 .AND. T .GT. 2.5E0) ERRT = T - 2.53E0
        IF (KPRINT .GT. 2) THEN
           WRITE (LUN,160) ERRT
        ENDIF
 160    FORMAT(1X,1P,' error in t location of root is',E12.4//)
        IF (ABS(ERRT) .LT. 1.0E-3) GO TO 170
        IF (KPRINT .GE. 2) THEN
           IPASS = 0
           WRITE (LUN,165)
        ENDIF
 165    FORMAT(//1X,' WARNING.. ROOT ERROR EXCEEDS 1.0E-3'//)
        NERR = NERR + 1
 170    CONTINUE
        GO TO 120
C
C IF NO ROOT FOUND, INCREMENT TOUT AND LOOP BACK.
 175    TOUT = TOUT + 1.0E0
 180    CONTINUE
C
C PROBLEM COMPLETE.  PRINT FINAL STATISTICS.
 185  CONTINUE
      IF (IDID .LT. 0) NERR = NERR + 1
      NST = IWORK(11)
      NJE = IWORK(13)
      NRE = IWORK(12)
      NGE = IWORK(16)
      LENRW = 0
      LENIW = 0
      NREA = NRE
      IF (JTYPE .EQ. 2) NRE = NRE + NEQ*NJE
C
      IF (KPRINT .GT. 2) THEN
         WRITE (LUN,190) NST,NRE,NREA,NJE,NGE,ERO
      ENDIF
 190  FORMAT(//1X,1P,' FINAL STATISTICS FOR THIS RUN..'/
     *  1X,' number of steps =',I5/
     *  1X,' number of r-s   =',I5/
     *  1X,' (excluding j-s) =',I5/
     *  1X,' number of j-s   =',I5/
     *  1X,' number of g-s   =',I5/
     *  1X,' error overrun   =',E10.2)
C
C-----------------------------------------------------------------------
C Second problem (Van Der Pol oscillator).
C The initial value problem is..
C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
C The root function is  G = Y1.
C An analytic solution is not known, but the zeros of Y1 are known
C to 15 figures for purposes of checking the accuracy.
C-----------------------------------------------------------------------
C
C RESET INFO ARRAY
C
      DO 195 I = 1,15
 195     INFO(I) = 0
C
C SET TOLERANCE PARAMETERS AND PRINT HEADING.
C NOTE THAT INFO(2) IS SET TO ONE, INDICATING THAT RTOL AND ATOL
C ARE ARRAYS.  EACH ENTRY OF RTOL AND ATOL MUST THEN BE DEFINED.
C
      INFO(2) = 1
      RTOL(1) = 1.0E-4
      RTOL(2) = 1.0E-4
      ATOL(1) = 1.0E-4
      ATOL(2) = 1.0E-3
      IF (KPRINT .GE. 2) THEN
         WRITE (LUN,200) RTOL(1),ATOL(1),ATOL(2)
      ENDIF
 200  FORMAT('1'/1X,1p,' DEMONSTRATION PROGRAM FOR SDASRT'///
     *  1X,' PROBLEM 2.. (Van Der Pol oscillator)'///
     *  1X,' problem is DY1/DT = Y2,  DY2/DT = 100*(1-Y1**2)*Y2 - Y1'/
     *  1X,'            Y1(0) = 2,  Y2(0) = 0'//
     *  1X,' root function is  G = Y1'//
     *  1X,'   RTOL =',E10.1,'   ATOL =',2E10.1///)
C
C NOTE: JTYPE INDICATES THE JACOBIAN TYPE:
C       JTYPE = 1 ==> JACOBIAN IS DENSE AND USER-SUPPLIED
C       JTYPE = 2 ==> JACOBIAN IS DENSE AND COMPUTED INTERNALLY
C
C LOOP OVER JTYPE = 1, 2.  SET REMAINING PARAMETERS AND PRINT JTYPE.
      DO 290 JTYPE = 1,2
C
C     SET INFO(1) = 0 TO INDICATE START OF A NEW PROBLEM
C     SET INFO(5) = 2-JTYPE TO TELL DASSL THE JACOBIAN TYPE.
C
      INFO(1) = 0
      INFO(5) = 2-JTYPE
      NEQ = 2
      T = 0.0E0
      Y(1) = 2.0E0
      Y(2) = 0.0E0
      YPRIME(1) = 0.E0
      YPRIME(2) = -2.0E0
      TOUT = 20.0E0
      LRW = 100
      LIW = 100
      NG = 1
      IF (KPRINT .GT. 2) THEN
        WRITE (LUN,210) JTYPE
      ENDIF
 210  FORMAT(////////1X,' SOLUTION WITH JTYPE =',I2////)
C
C CALL SDASRT IN LOOP OVER TOUT VALUES 20,40,...,200.
      DO 270 IOUT = 1,10
 220    CONTINUE
C
        CALL SDASRT (RES2,NEQ,T,Y,YPRIME,TOUT,
     *  INFO,RTOL,ATOL,IDID,
     *  RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC2,
     *  GR2,NG,JROOT)
C
C PRINT Y1 AND Y2.
        IF (KPRINT .GT. 2) THEN
           WRITE (LUN,230) T,Y(1),Y(2)
        ENDIF
 230    FORMAT(1X,1P,' at t =',E15.7,5X,'y1 =',E15.7,5X,'y2 =',E15.7)
        IF (IDID .LT. 0) GO TO 275
        IF (IDID .NE. 4) GO TO 265
C
C IF A ROOT WAS FOUND, WRITE RESULTS AND CHECK ROOT LOCATION.
C THEN RETURN TO SDASRT TO CONTINUE THE INTEGRATION.
        IF (KPRINT .GT. 2) THEN
           WRITE (LUN,240) T
        ENDIF
 240    FORMAT(/1X,1P,' root found at t =',E15.7)
        KROOT = INT(T/81.2E0 + 0.5E0)
        TZERO = 81.17237787055E0 + FLOAT(KROOT-1)*81.41853556212E0
        ERRT = T - TZERO
        IF (KPRINT .GT. 2) THEN
           WRITE (LUN,250) ERRT
        ENDIF
 250    FORMAT(1X,1P,' error in t location of root is',E12.4//)
        IF (ERRT .LT. 1.0E0) GO TO 260
        IF (KPRINT .GE. 2) THEN
          IPASS = 0
          WRITE (LUN,255)
        ENDIF
 255    FORMAT(//1X,' WARNING.. ROOT ERROR EXCEEDS 1.0'//)
        NERR = NERR + 1
 260    CONTINUE
        GO TO 220
C
C IF NO ROOT FOUND, INCREMENT TOUT AND LOOP BACK.
 265    TOUT = TOUT + 20.0E0
 270    CONTINUE
C
C PROBLEM COMPLETE.  PRINT FINAL STATISTICS.
 275  CONTINUE
      IF (IDID .LT. 0) NERR = NERR + 1
      NST = IWORK(11)
      NJE = IWORK(13)
      NRE = IWORK(12)
      NGE = IWORK(16)
      LENRW = 0
      LENIW = 0
      NREA = NRE
      IF (JTYPE .EQ. 2) NRE = NRE + NEQ*NJE
      IF (KPRINT .GE. 2) THEN
         WRITE (LUN,280) NST,NRE,NREA,NJE,NGE
      ENDIF
 280  FORMAT(//1X,1P,' FINAL STATISTICS FOR THIS RUN..'/
     *  1X,' number of steps =',I5/
     *  1X,' number of r-s   =',I5/
     *  1X,' (excluding j-s) =',I5/
     *  1X,' number of j-s   =',I5/
     *  1X,' number of g-s   =',I5)
 290  CONTINUE
C
C
      IF (KPRINT .GE. 2) THEN
         WRITE (LUN,300) NERR
      ENDIF
 300  FORMAT(////1X,' NUMBER OF ERRORS ENCOUNTERED =',I3)
C
      IF (NERR .GT. 0) THEN
        IPASS = 0
      ENDIF
      IF ((IPASS .EQ. 1) .AND. (KPRINT .GT. 1)) WRITE (LUN,700)
      IF ((IPASS .EQ. 0) .AND. (KPRINT .NE. 0)) WRITE (LUN,800)
 700  FORMAT (/,' **********SDASRT PASSED ALL TESTS**********')
 800  FORMAT (/,' **********SDASRT FAILED SOME TESTS*********')
      STOP
      END
C
      SUBROUTINE RES1(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
      IMPLICIT REAL (A-H,O-Z)
      COMMON /LOCAL/NEQ
      DIMENSION Y(NEQ), YPRIME(NEQ), DELTA(NEQ)
C
C     CHECK Y TO MAKE SURE THAT IT IS VALID INPUT.
C     IF Y IS LESS THAN OR EQUAL TO ZERO, THIS IS INVALID INPUT.
C
      IF (Y(1) .LE. 0.0E0) THEN
         IRES = -1
      ELSE
C
C        CALL F TO OBTAIN F(T,Y)
C
         CALL F1(NEQ,T,Y,DELTA)
C
C        FORM F = F'-F(T,Y)
C
         DO 10 I = 1,NEQ
            DELTA(I) = YPRIME(I) - DELTA(I)
 10      CONTINUE
      ENDIF
C
      RETURN
      END
C
      SUBROUTINE F1 (NEQ, T, Y, YDOT)
      IMPLICIT REAL (A-H,O-Z)
      INTEGER NEQ
      REAL T, Y, YDOT
      DIMENSION Y(1), YDOT(1)
      YDOT(1) = ((2.0E0*LOG(Y(1)) + 8.0E0)/T - 5.0E0)*Y(1)
      RETURN
      END
C
      SUBROUTINE GR1 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)
      IMPLICIT REAL (A-H,O-Z)
      INTEGER NEQ, NG
      REAL T, Y, GROOT
      DIMENSION Y(1), GROOT(2)
      GROOT(1) = ((2.0E0*LOG(Y(1)) + 8.0E0)/T - 5.0E0)*Y(1)
      GROOT(2) = LOG(Y(1)) - 2.2491E0
      RETURN
      END
C
      SUBROUTINE RES2(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
      IMPLICIT REAL (A-H,O-Z)
      COMMON /LOCAL/NEQ
      DIMENSION Y(NEQ), YPRIME(NEQ), DELTA(NEQ)
C
C     CALL F TO OBTAIN F(T,Y)
C
      CALL F2(NEQ,T,Y,DELTA)
C
C     FORM F = F'-F(T,Y)
C
      DO 10 I = 1,NEQ
         DELTA(I) = YPRIME(I) - DELTA(I)
 10   CONTINUE
C
      RETURN
      END
C
      SUBROUTINE F2 (NEQ, T, Y, YDOT)
      IMPLICIT REAL (A-H,O-Z)
      INTEGER NEQ
      REAL T, Y, YDOT
      DIMENSION Y(2), YDOT(2)
      YDOT(1) = Y(2)
      YDOT(2) = 100.0E0*(1.0E0 - Y(1)*Y(1))*Y(2) - Y(1)
      RETURN
      END
C
      SUBROUTINE JAC2 (T, Y, YPRIME, PD, CJ, RPAR, IPAR)

      IMPLICIT REAL (A-H,O-Z)
      INTEGER NEQ,  NROWPD
      REAL T, Y, PD
      PARAMETER (NROWPD=2)
      DIMENSION Y(2), PD(NROWPD,2)
      COMMON /LOCAL/NEQ
C
C FIRST DEFINE THE JACOBIAN MATRIX FOR THE RIGHT HAND SIDE
C OF THE ODE: F' = F(T,Y) , I.E. DF/DY)
C
      PD(1,1) = 0.0E0
      PD(1,2) = 1.0E0
      PD(2,1) = -200.0E0*Y(1)*Y(2) - 1.0E0
      PD(2,2) = 100.0E0*(1.0E0 - Y(1)*Y(1))
C
C NEXT UPDATE THE JACOBIAN WITH THE RIGHT HAND SIDE TO FORM THE
C DAE JACOBIAN: D(F'-F)/DY = DF'/DY - DF/DY = I - DF/DY
C
      PD(1,1) = CJ - PD(1,1)
      PD(1,2) =    - PD(1,2)
      PD(2,1) =    - PD(2,1)
      PD(2,2) = CJ - PD(2,2)
C
      RETURN
      END
C
      SUBROUTINE GR2 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)
      IMPLICIT REAL (A-H,O-Z)
      INTEGER NEQ, NG
      REAL T, Y, GROOT
      DIMENSION Y(2), GROOT(1)
      GROOT(1) = Y(1)
      RETURN
      END