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