C C Main program to execute SDEMO demonstration program. C C C***BEGIN PROLOGUE SDEMO C***DATE WRITTEN 851113 (YYMMDD) C***REVISION DATE 900904 (YYMMDD) C***CATEGORY NO. D2A2 C***KEYWORDS LIBRARY=SLATEC,QUICK CHECK,SDEMO, C TYPE=SINGLE 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 SDASSL. c***DESCRIPTION C Demonstration program for SDASSL. C This version is in single precision. C C SDASSL 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 SDASSL,RES1,JAC1,RES2,JAC2 C***END PROLOGUE SDEMO C EXTERNAL RES1, JAC1, RES2, JAC2 DIMENSION Y(25), RWORK(550), IWORK(45), YPRIME(25), DELTA(25) DIMENSION INFO(15) DATA TOUT1/1.0E0/, DTOUT/1.0E0/ C LUN = 6 KPRINT = 3 IPASS = 1 NERR = 0 RTOL = 0.0E0 ATOL = 1.0E-3 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 SDASSL',/// 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 C 120 FORMAT(////1X,' INFO(5) =',I3// 1 6X,'T',15X,'X1',14X,'X2',12X,'NQ',6X,'H',12X/) T = 0.0E0 Y(1) = 1.0E0 Y(2) = 0.0E0 YPRIME(1) = -10.0E0 YPRIME(2) = 10.0E0 TOUT = TOUT1 ERO = 0.0E0 DO 170 IOUT = 1,NOUT CALL SDASSL(RES1,NEQ,T,Y,YPRIME,TOUT,INFO,RTOL,ATOL,IDID, 1 RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC1) 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 = EXP(-10.0E0*T) YT2 = 1.0E0 - YT1 ER1 = ABS(YT1 - Y(1)) ER2 = ABS(YT2 - Y(2)) ER = AMAX1(ER1,ER2)/ATOL ERO = AMAX1(ERO,ER) IF (ER .LT. 1000.0E0) GO TO 170 IF (KPRINT .GE. 2) THEN WRITE (LUN,150) IPASS = 0 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 SDASSL',/// 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.0E0 DO 230 I = 2,NEQ 230 Y(I) = 0.0E0 Y(1) = 1.0E0 DO 235 I = 1,NEQ 235 DELTA(I) = 0.0E0 CALL RES2(T,Y,DELTA,YPRIME,IRES,RPAR,IPAR) TOUT = 0.01E0 ERO = 0.0E0 DO 270 IOUT = 1,NOUT CALL SDASSL(RES2,NEQ,T,Y,YPRIME,TOUT,INFO,RTOL,ATOL,IDID, 1 RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC2) CALL EDIT2(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 = AMAX1(ERO,ER) IF (ER .LE. 1000.0E0) GO TO 270 IF (KPRINT .GE. 2) THEN WRITE (LUN,150) IPASS = 0 ENDIF C NERR = NERR + 1 270 TOUT = TOUT*10.0E0 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 C 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 (/,' **********SDASSL PASSED ALL TESTS**********') 800 FORMAT (/,' **********SDASSL FAILED SOME TESTS*********') STOP END SUBROUTINE RES1(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR) DIMENSION Y(1),YPRIME(1),DELTA(1) DELTA(1) = YPRIME(1) + 10.0E0*Y(1) DELTA(2) = Y(2) + Y(1) - 1.0E0 RETURN END SUBROUTINE JAC1(T,Y,YPRIME,PD,CJ,RPAR,IPAR) DIMENSION Y(1),YPRIME(1),PD(2,1) PD(1,1) = CJ + 10.0E0 PD(1,2) = 0.0E0 PD(2,1) = 1.0E0 PD(2,2) = 1.0E0 RETURN END SUBROUTINE RES2(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR) DIMENSION Y(1),YPRIME(1),DELTA(1) DATA ALPH1/1.0E0/, ALPH2/1.0E0/, NG/5/ DO 10 J = 1,NG DO 10 I = 1,NG K = I + (J - 1)*NG D = -2.0E0*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 JAC2(T,Y,YPRIME,PD,CJ,RPAR,IPAR) DIMENSION Y(1), PD(11,1), YPRIME(1) DATA ALPH1/1.0E0/, ALPH2/1.0E0/, 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.0E0 - CJ PD(MBANDP1,J) = ALPH1 PD(MBANDP2,J) = 0.0E0 PD(MBANDP3,J) = 0.0E0 PD(MBANDP4,J) = 0.0E0 10 PD(MBANDP5,J) = ALPH2 DO 20 J = 1,NEQ,NG 20 PD(MBANDP1,J) = 0.0E0 RETURN END SUBROUTINE EDIT2 (Y, T, ERM) DIMENSION Y(25) DATA ALPH1/1.0E0/, ALPH2/1.0E0/, NG/5/ ERM = 0.0E0 IF (T .EQ. 0.0E0) RETURN EX = 0.0E0 IF (T .LE. 30.0E0) EX = EXP(-2.0E0*T) A2 = 1.0E0 DO 60 J = 1,NG A1 = 1.0E0 DO 50 I = 1,NG K = I + (J - 1)*NG YT = T**(I+J-2)*EX*A1*A2 ER = ABS(Y(K)-YT) ERM = AMAX1(ERM,ER) A1 = A1*ALPH1/FLOAT(I) 50 CONTINUE A2 = A2*ALPH2/FLOAT(J) 60 CONTINUE RETURN END