C
C Main program to execute DDEMO demonstration program.
C
C***BEGIN PROLOGUE  DDEMO
C***DATE WRITTEN   851210     (YYMMDD)
C***REVISION DATE  900904     (YYMMDD)
C***CATEGORY NO.   D2A2
C***KEYWORDS  LIBRARY=SLATEC,QUICK CHECK,DDEMO,
C             TYPE=DOUBLE
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 DDASSL.
c***DESCRIPTION
C       Demonstration program for DDASSL.
C       This version is in double precision.
C
C       DDASSL is used to solve two simple problems,
C       one with a full jacobian, the other with a banded jacobian,
C       with the 2 appropriate values of info(5) in each case.
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  DDASSL,DRES1,DJAC1,DRES2,DJAC2
C***END PROLOGUE DDEMO
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      EXTERNAL DRES1, DJAC1, DRES2, DJAC2
      DIMENSION Y(25), RWORK(550), IWORK(45), YPRIME(25), DELTA(25)
      DIMENSION INFO(15)
      DATA TOUT1/1.0D0/, DTOUT/1.0D0/
C
      LUN = 6
      KPRINT = 3
      IPASS = 1
      NERR = 0
      RTOL = 0.0D0
      ATOL = 1.0D-6
      LRW = 550
      LIW = 45
C
C FIRST PROBLEM
C
      NEQ = 2
      NOUT = 10
      IF (KPRINT .GE. 2) THEN
         WRITE (LUN,110) NEQ,RTOL,ATOL
      ENDIF
C
 110  FORMAT('1'/1X,' DEMONSTRATION PROGRAM FOR DDASSL',///
     1  1X,' PROBLEM 1..   LINEAR DIFFERENTIAL/ALGEBRAIC SYSTEM..',/
     2  1X,' X1DOT + 10.0*X1 = 0,  X1 + X2 = 1',/
     3  1X,' X1(0) = 1.0, X2(0) = 0.0',/
     4  1X,' NEQ =',I2/
     5  1X,' RTOL =',E10.1,'   ATOL =',E10.1)
C
      DO 190 J190 = 1,2
      DO 115 I = 1,15
 115    INFO(I) = 0
      IF(J190 .EQ. 2) INFO(5) = 1
      IF (KPRINT .GT. 2) THEN
         WRITE (LUN,120) INFO(5)
      ENDIF
 120  FORMAT(////1X,' INFO(5) =',I3//
     1  6X,'T',15X,'X1',14X,'X2',12X,'NQ',6X,'H',12X/)
      T = 0.0D0
      Y(1) = 1.0D0
      Y(2) = 0.0D0
      YPRIME(1) = -10.0D0
      YPRIME(2) =  10.0D0
      TOUT = TOUT1
      ERO = 0.0D0
      DO 170 IOUT = 1,NOUT
        CALL DDASSL(DRES1,NEQ,T,Y,YPRIME,TOUT,INFO,RTOL,ATOL,IDID,
     1     RWORK,LRW,IWORK,LIW,RPAR,IPAR,DJAC1)
        HU = RWORK(7)
        NQU = IWORK(8)
        IF (KPRINT .GT. 2) THEN
          WRITE (LUN,140) T,Y(1),Y(2),NQU,HU
        ENDIF
C
 140    FORMAT(1X,E15.5,E16.5,E16.5,I6,E14.3)
        IF (IDID .LT. 0) GO TO 175
        YT1 = DEXP(-10.0D0*T)
        YT2 = 1.0D0 - YT1
        ER1 = DABS(YT1 - Y(1))
        ER2 = DABS(YT2 - Y(2))
        ER = DMAX1(ER1,ER2)/ATOL
        ERO = DMAX1(ERO,ER)
        IF (ER .LT. 1000.0D0) GO TO 170
        IF (KPRINT .GE. 2) THEN
          IPASS = 0
          WRITE (LUN,150)
        ENDIF
C
 150    FORMAT(//1X,' WARNING.. ERROR EXCEEDS 1000 * TOLERANCE',//)
        NERR = NERR + 1
 170    TOUT = TOUT + DTOUT
 175  CONTINUE
      IF (IDID .LT. 0) NERR = NERR + 1
      NST = IWORK(11)
      NFE = IWORK(12)
      NJE = IWORK(13)
      IF (KPRINT .GT. 2) THEN
        WRITE (LUN,180) NST,NFE,NJE,ERO
      ENDIF
C
 180  FORMAT(//1X,' FINAL STATISTICS FOR THIS RUN..',/
     1  1X,' NUMBER OF STEPS =',I5/
     2  1X,' NUMBER OF F-S   =',I5/
     4  1X,' NUMBER OF J-S   =',I5/
     5  1X,' ERROR OVERRUN =',E10.2)
 190  CONTINUE
C
C SECOND PROBLEM
C
      NEQ = 25
      ML = 5
      MU = 0
      IWORK(1) = ML
      IWORK(2) = MU
      NOUT = 5
      IF (KPRINT .GE. 2) THEN
        WRITE (LUN,210) NEQ,ML,MU,RTOL,ATOL
      ENDIF
C
 210  FORMAT('1'/1X,' DEMONSTRATION PROGRAM FOR DDASSL',///
     1  1X,' PROBLEM 2.. YDOT = A * Y , WHERE ',
     2  ' A IS A BANDED LOWER TRIANGULAR MATRIX',/
     2  1X,'  DERIVED FROM 2-D ADVECTION PDE',/
     3  1X,' NEQ =',I3,'   ML =',I2,'   MU =',I2/
     4  1X,' RTOL =',E10.1,'   ATOL =',E10.1)
      DO 290 J290 = 1,2
      DO 215 I = 1,15
 215    INFO(I) = 0
      INFO(6) = 1
      IF(J290 .EQ. 2) INFO(5) = 1
      IF (KPRINT .GT. 2) THEN
        WRITE (LUN,220) INFO(5)
      ENDIF
C
 220  FORMAT(////1X,' INFO(5) =',I3//
     1  6X,'T',14X,'MAX.ERR.',5X,'NQ',6X,'H'/)
      T = 0.0D0
      DO 230 I = 2,NEQ
 230    Y(I) = 0.0D0
      Y(1) = 1.0D0
      DO 235 I = 1,NEQ
 235    DELTA(I) = 0.0D0
      CALL DRES2(T,Y,DELTA,YPRIME,IRES,RPAR,IPAR)
      TOUT = 0.01D0
      ERO = 0.0D0
      DO 270 IOUT = 1,NOUT
        CALL DDASSL(DRES2,NEQ,T,Y,YPRIME,TOUT,INFO,RTOL,ATOL,IDID,
     1     RWORK,LRW,IWORK,LIW,RPAR,IPAR,DJAC2)
        CALL DEDIT2(Y,T,ERM)
        HU = RWORK(7)
        NQU = IWORK(8)
        IF (KPRINT .GT. 2) THEN
          WRITE (LUN,240) T,ERM,NQU,HU
        ENDIF
C
 240    FORMAT(1X,E15.5,E14.3,I6,E14.3)
        IF (IDID .LT. 0) GO TO 275
        ER = ERM/ATOL
        ERO = DMAX1(ERO,ER)
        IF (ER .LE. 1000.0D0) GO TO 270
        IF (KPRINT .GE. 2) THEN
          IPASS = 0
          WRITE (LUN,150)
        ENDIF
        NERR = NERR + 1
 270    TOUT = TOUT*10.0D0
 275  CONTINUE
      IF (IDID .LT. 0) NERR = NERR + 1
      NST = IWORK(11)
      NFE = IWORK(12)
      NJE = IWORK(13)
      IF (KPRINT .GT. 2) THEN
        WRITE (LUN,180) NST,NFE,NJE,ERO
      ENDIF
 290  CONTINUE
      IF (KPRINT .GE. 2) THEN
        WRITE (LUN,300) NERR
      ENDIF
C
 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 (/,' **********DDASSL PASSED ALL TESTS**********')
 800  FORMAT (/,' **********DDASSL FAILED SOME TESTS*********')
      STOP
      END
      SUBROUTINE DRES1(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Y(1),YPRIME(1),DELTA(1)
      DELTA(1) = YPRIME(1) + 10.0D0*Y(1)
      DELTA(2) = Y(2) + Y(1) - 1.0D0
      RETURN
      END
      SUBROUTINE DJAC1(T,Y,YPRIME,PD,CJ,RPAR,IPAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Y(1),YPRIME(1),PD(2,1)
      PD(1,1) = CJ + 10.0D0
      PD(1,2) = 0.0D0
      PD(2,1) = 1.0D0
      PD(2,2) = 1.0D0
      RETURN
      END
      SUBROUTINE DRES2(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Y(1),YPRIME(1),DELTA(1)
      DATA ALPH1/1.0D0/, ALPH2/1.0D0/, NG/5/
      DO 10 J = 1,NG
      DO 10 I = 1,NG
        K = I + (J - 1)*NG
        D = -2.0D0*Y(K)
        IF (I .NE. 1) D = D + Y(K-1)*ALPH1
        IF (J .NE. 1) D = D + Y(K-NG)*ALPH2
 10     DELTA(K) = D - YPRIME(K)
      RETURN
      END
      SUBROUTINE DJAC2(T,Y,YPRIME,PD,CJ,RPAR,IPAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Y(1), PD(11,1), YPRIME(1)
      DATA ALPH1/1.0D0/, ALPH2/1.0D0/, NG/5/
      DATA ML/5/, MU/0/, NEQ/25/
      MBAND = ML + MU + 1
      MBANDP1 = MBAND + 1
      MBANDP2 = MBAND + 2
      MBANDP3 = MBAND + 3
      MBANDP4 = MBAND + 4
      MBANDP5 = MBAND + 5
      DO 10 J = 1,NEQ
        PD(MBAND,J) = -2.0D0 - CJ
        PD(MBANDP1,J) = ALPH1
        PD(MBANDP2,J) = 0.0D0
        PD(MBANDP3,J) = 0.0D0
        PD(MBANDP4,J) = 0.0D0
 10     PD(MBANDP5,J) = ALPH2
      DO 20 J = 1,NEQ,NG
 20     PD(MBANDP1,J) = 0.0D0
      RETURN
      END
      SUBROUTINE DEDIT2 (Y, T, ERM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Y(25)
      DATA ALPH1/1.0D0/, ALPH2/1.0D0/, NG/5/
      ERM = 0.0D0
      IF (T .EQ. 0.0D0) RETURN
      EX = 0.0D0
      IF (T .LE. 30.0D0) EX = DEXP(-2.0D0*T)
      A2 = 1.0D0
      DO 60 J = 1,NG
        A1 = 1.0D0
        DO 50 I = 1,NG
          K = I + (J - 1)*NG
          YT = T**(I+J-2)*EX*A1*A2
          ER = DABS(Y(K)-YT)
          ERM = DMAX1(ERM,ER)
          A1 = A1*ALPH1/DFLOAT(I)
 50       CONTINUE
        A2 = A2*ALPH2/DFLOAT(J)
 60     CONTINUE
      RETURN
      END