      PROGRAM DASSLT
C...
C...  AN INTERFACE BETWEEN DSS/2 APPLICATIONS AND INTEGRATOR DASSL
C...
C...  THE FOLLOWING PROGRAM CALLS A DSS/2 APPLICATION, DEFINED AS THE
C...  USER-SUPPLIED SUBROUTINES INITAL, DERV AND PRINT, PLUS DATA.
C...  THE SYSTEM OF ODES PROGRAMMED IN SUBROUTINE DERV IS INTEGRATED
C...  BY DASSL (A SOLVER FOR ALGEBRAIC/ODE SYSTEMS DEVELOPED BY LINDA
C...  R. PETZOLD, LAWRENCE LIVERMORE NATIONAL LABORATORY).
C...
C...  THE MODEL INITIAL CONDITIONS ARE SET IN SUBROUTINE INITAL, AND
C...  THE MODEL DERIVATIVES ARE PROGRAMMED IN SUBROUTINE DERV.  THE
C...  NUMERICAL SOLUTION IS PRINTED AND PLOTTED IN SUBROUTINE PRINT.
C...
C...  TYPE ALL REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  SUBROUTINES INITAL, DERV AND PRINT.  THE FOLLOWING CODING IS FOR
C...  250 ORDINARY DIFFERENTIAL EQUATIONS (ODES).  IF MORE ODES ARE TO
C...  BE INTEGRATED, ALL OF THE 250*S SHOULD BE CHANGED TO THE REQUIRED
C...  NUMBER
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    Y(250)
     2      /F/    F(250)
C...
C...  THE NUMBER OF DIFFERENTIAL EQUATIONS IS AVAILABLE THROUGH COMMON
C...  /N/ FOR USE IN SUBROUTINE RES
      COMMON/N/      NEQN
C...
C...  COMMON AREA TO PROVIDE THE INPUT/OUTPUT UNIT NUMBERS TO OTHER
C...  SUBROUTINES
      COMMON/IO/       NI,        NO
C...
C...  ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY DASSL.  ARRAY
C...  WORK IS SIZED BY THE FOLLOWING FORMULA FOR BANDED PROCESSING
C...  OF THE ODE JACOBIAN MATRIX VIA FINITE DIFFERENCES: 40+(MAXORD
C...  +4)*NEQN+(2*ML+MU+1)*NEQN+2*(NEQN/(ML+MU+1)+1).  FOR MAXORD = 5,
C...  NEQN = 250, ML = MU = 3, 40+(5+4)*250+(2*3+3+1)*250+2*(250/(3+3
C...  +1)+1) = 5000 (APPROXIMATELY).  IF FEWER THAN 250 EQUATIONS ARE
C...  TO BE SOLVED, THEN ML AND MU CAN BE LARGER FOR THE GIVEN SIZE OF
C...  WORK.  IF NEQN GT 250 AND/OR ML AND MU GT 3, WORK MUST BE DIMEN-
C...  SIONED TO A LARGER SIZE
      DIMENSION YV(250), YPRIME(250), DELTA(250),
     1      RWORK(5000),  IWORK(270),   INFO(15)
C...
C...  EXTERNAL THE DERIVATIVE AND ODE JACOBIAN MATRIX ROUTINES CALLED BY
C...  DASSL
      EXTERNAL RES, JAC
C...
C...  ARRAY FOR THE TITLE (FIRST LINE OF DATA), CHARACTERS  END OF RUNS
      CHARACTER TITLE(20)*4, ENDRUN(3)*4
C...
C...  VARIABLE FOR THE TYPE OF ERROR CRITERION
      CHARACTER*3 ABSREL
C...
C...  DEFINE THE CHARACTERS  END OF RUNS
      DATA ENDRUN/'END ','OF RU','NS  '/
C...
C...  DEFINE THE INPUT/OUTPUT UNIT NUMBERS
      NI=5
      NO=6
C...
C...  OPEN INPUT AND OUTPUT FILES
      OPEN(NI,FILE='DATA',  STATUS='OLD')
      OPEN(NO,FILE='OUTPUT',STATUS='NEW')
C...
C...  INITIALIZE THE RUN COUNTER
      NORUN=0
C...
C...  BEGIN A RUN
1     NORUN=NORUN+1
C...
C...  INITIALIZE THE RUN TERMINATION VARIABLE
      NSTOP=0
C...
C...  READ THE FIRST LINE OF DATA
      READ(NI,1000,END=999)(TITLE(I),I=1,20)
C...
C...  TEST FOR  END OF RUNS  IN THE DATA
      DO 2 I=1,3
      IF(TITLE(I).NE.ENDRUN(I))GO TO 3
2     CONTINUE
C...
C...  AN END OF RUNS HAS BEEN READ, SO TERMINATE EXECUTION
999   STOP
C...
C...  READ THE SECOND LINE OF DATA
3     READ(NI,1001,END=999)T0,TF,TP
C...
C...  READ THE THIRD LINE OF DATA
      READ(NI,1002,END=999)NEQN,NMAX,INT,IERROR,ABSREL,ERROR
C...
C...  PRINT A DATA SUMMARY
      WRITE(NO,1003)NORUN,(TITLE(I),I=1,20),
     1                     T0,TF,TP,
     2                     NEQN,NMAX,INT,IERROR,ABSREL,ERROR
C...
C...  INITIALIZE TIME
      T=T0
C...
C...  SET THE INITIAL CONDITIONS
      CALL INITAL
C...
C...  PRINT THE INITIAL CONDITIONS
      CALL PRINT(NI,NO)
C...
C...  SET THE INITIAL CONDITIONS FOR SUBROUTINE DASSL
      CALL DERV
      DO 5 I=1,NEQN
      YV(I)=Y(I)
      YPRIME(I)=F(I)
5     CONTINUE
C...
C...  SET THE PARAMETERS FOR SUBROUTINE DASSL
      TV=T0
      RTOL=ERROR
      ATOL=1.0D-10
      LRW=5000
      LIW=270
      DO 7 I=1,15
      INFO(I)=0
7     CONTINUE
C...
C...  FULL JACOBIAN MATRIX
      IF(NORUN.EQ.1)THEN
C...
C...  BANDED JACOBIAN MATRIX
      ELSE IF(NORUN.EQ.2)THEN
         INFO(6)=1
         IWORK(1)=1
         IWORK(2)=1
      END IF
C...
C...  INITIATE THE INTEGRATION
      TOUT=TV+TP
C...
C...  CALL SUBROUTINE DASSL TO COVER ONE PRINT INTERVAL
4     CALL DASSL(RES,NEQN,TV,YV,YPRIME,TOUT,INFO,RTOL,ATOL,IDID,
     1           RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC)
C...
C...  PRINT THE SOLUTION
      T=TV
      DO 6 I=1,NEQN
      Y(I)=YV(I)
6     CONTINUE
      CALL PRINT(NI,NO)
C...
C...  TEST FOR AN ERROR CONDITION
      IF(IDID.LT.0)THEN
C...
C...     PRINT A MESSAGE INDICATING AN ERROR CONDITION
         WRITE(NO,1004)IDID
C...
C...     GO ON TO THE NEXT RUN
         GO TO 1
      END IF
C...
C...  CHECK FOR A RUN TERMINATION
      IF(NSTOP.NE.0)GO TO 1
C...
C...  THE INTEGRATION HAS PROCEEDED SATISFACTORILY, SO PREPARE FOR THE
C...  NEXT INTERVAL IN T.  NOTE THAT THE FOLLOWING CALCULATION OF THE
C...  OUTPUT TIME, TOUT, PRODUCES THREE OUTPUT POINTS FOR EACH DECADE
C...  IN T
      TOUT=TV*(1.0D+01)**(1.0D+00/3.0D+00)
C...
C...  CHECK FOR THE END OF THE RUN
      IF(TV.LT.(TF-0.5D+00*TP))GO TO 4
C...
C...  THE CURRENT RUN IS COMPLETE, SO PRINT THE COMPUTATIONAL STAT-
C...  ISTICS FOR DASSL AND GO ON TO THE NEXT RUN
      WRITE(NO,8)RWORK(7),IWORK(8),IWORK(11),IWORK(12),IWORK(13)
8     FORMAT(1H ,//,' COMPUTATIONAL STATISTICS'                    ,//,
     1 ' LAST STEP SIZE                       ',              D10.3,//,
     2 ' LAST ORDER OF THE METHOD             ',                I10,//,
     3 ' TOTAL NUMBER OF STEPS TAKEN          ',                I10,//,
     4 ' NUMBER OF FUNCTION EVALUATIONS       ',                I10,//,
     5 ' NUMBER OF JACOBIAN EVALUATIONS       ',                I10,/)
      GO TO 1
C...
C...  *****************************************************************
C...
C...  FORMATS
C...
1000  FORMAT(20A4)
1001  FORMAT(3D10.0)
1002  FORMAT(4I5,2X,A3,D10.0)
1003  FORMAT(1H1,
     1 ' RUN NO. - ',I3,2X,20A4,//,
     2 ' INITIAL T - ',D10.3,//,
     3 '   FINAL T - ',D10.3,//,
     4 '   PRINT T - ',D10.3,//,
     5 ' NUMBER OF DIFFERENTIAL EQUATIONS - ',I3,//,
     6 ' PRINT INTERVAL/MINIMUM INTEGRATION INTERVAL - ',I5,//,
     7 ' INTEGRATION ALGORITHM - ',I2,' - DASSL ',//,
     8 ' INTEGRATION ERROR MESSAGES - ',I1,//,
     9 ' ERROR CRITERION - ',A3,//,
     A ' MAXIMUM INTEGRATION ERROR - ',D10.3,//,
     B 1H1)
1004  FORMAT(1H ,//,' IDID = ',I3,//,
     1 ' INDICATING AN INTEGRATION ERROR, SO THE CURRENT RUN'     ,/,
     2 ' IS TERMINATED.  PLEASE REFER TO THE DOCUMENTATION FOR'   ,/,
     3 ' SUBROUTINE',//,25X,'DASSL',//,
     4 ' FOR AN EXPLANATION OF THESE ERROR INDICATORS'            )
      END
      SUBROUTINE RES(TV,YV,YPRIME,DELTA,IRES,RPAR,IPAR)
C...
C...  SUBROUTINE RES IS AN INTERFACE ROUTINE BETWEEN SUBROUTINES DASSL
C...  AND DERV
C...
C...  DSS/2 COMMON AREA
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(1)
     2      /F/      F(1)
C...
C...  THE NUMBER OF DIFFERENTIAL EQUATIONS IS AVAILABLE THROUGH COMMON
C...  /N/
      COMMON/N/      NEQN
C...
C...  ABSOLUTE DIMENSION THE DEPENDENT, DERIVATIVE AND RESIDUAL ARRAYS
      DIMENSION YV(250), YPRIME(250), DELTA(250)
C...
C...  TRANSFER THE INDEPENDENT VARIABLE, DEPENDENT VARIABLE VECTOR
C...  FOR USE IN SUBROUTINE DERV
      T=TV
      DO 1 I=1,NEQN
      Y(I)=YV(I)
1     CONTINUE
C...
C...  EVALUATE THE DERIVATIVE VECTOR
      CALL DERV
C...
C...  TRANSFER THE DERIVATIVE VECTOR FOR USE BY SUBROUTINE DASSL
      DO 2 I=1,NEQN
      DELTA(I)=YPRIME(I)-F(I)
2     CONTINUE
      RETURN
      END
      SUBROUTINE JAC
C...
C...  SUBROUTINE JAC IS A DUMMY ROUTINE TO SATISFY THE LOADER (SINCE JAC
C...  IS DECLARED AS AN EXTERNAL IN THE MAIN PROGRAM DASSLT).  JAC IS NOT
C...  ACTUALLY CALLED UNLESS AN OPTION OF DASSL IS SELECTED WHICH REQUIRES
C...  THE USER TO PROVIDE THE ODE ANALYTICAL JACOBIAN MATRIX.  THIS IS
C...  USUALLY NOT THE CASE SINCE FOR MOST PROBLEMS, THE JACOBIAN MATRIX
C...  IS CALCULATED INTERNALLY BY DASSL USING FINITE DIFFERENCES
      RETURN
      END
      SUBROUTINE INITAL
C...
C...  TWO LINEAR CONSTANT COEFFICIENT ORDINARY DIFFERENTIAL EQUATIONS
C...  WITH A STIFFNESS RATIO OF ONE MILLION - PRINTING OF THE NUMERI-
C...  CAL AND ANALYTICAL SOLUTIONS
C...
C...  THE TWO SIMULTANEOUS, LINEAR, CONSTANT COEFFICIENT ORDINARY
C...  DIFFERENTIAL EQUATIONS
C...
C...     DY1/DT = A11*Y1 + A12*Y2                                    (1)
C...
C...     DY2/DT = A21*Y1 + A22*Y2                                    (2)
C...
C...  ARE INTEGRATED SUBJECT TO THE INITIAL CONDITIONS
C...
C...     Y1(0) = 0, Y2(0) = 2                                     (3)(4)
C...
C...  THE ANALYTICAL SOLUTION TO EQUATIONS (1) TO (4) CAN EASILY BE
C...  DERIVED.  FOR EXAMPLE, IF Y1(T) AND Y2(T) ARE ASSUMED TO BE
C...
C...     Y1 = C1*EXP(L*T), Y2 = C2*EXP(L*T)
C...
C...  SUBSTITUTION IN EQUATIONS (1) AND (2) GIVES
C...
C...     L*C1*EXP(L*T) = A11*C1*EXP(L*T) + A12*C2*EXP(L*T)           (5)
C...
C...     L*C2*EXP(L*T) = A21*C1*EXP(L*T) + A22*C2*EX(L*T)            (6)
C...
C...  THE EXPONENTIAL EXP(L*T) CANCELS FROM EQUATIONS (5) AND (6), AND
C...  THE FOLLOWING HOMOGENEOUS LINEAR ALGEBRAIC EQUATIONS FOR C1 AND C2
C...  RESULT
C...
C...     (A11 - L)*C1 +       A12*C2 = 0                             (7)
C...
C...           A21*C1 + (A22 - L)*C2 = 0                             (8)
C...
C...  EQUATIONS (7) AND (8) CONSTITUTE THE CLASSICAL, SECOND-ORDER,
C...  LINEAR ALGERAIC EIGENVALUE PROBLEM WHICH CAN BE SOLVED FOR THE
C...  VALUES OF L THAT DO NOT GIVE THE TRIVIAL SOLUTION C1 = 0, C2 = 0.
C...  THIS WILL BE DONE BY SETTING THE DETERMINANT OF THE COEFFICIENT
C...  MATRIX EQUAL TO ZERO (OTHER ALGORITHMS FOR CALCULATING THE EIGEN-
C...  VALUES, L1 AND L2, ARE AVAILABLE)
C...
C...     (A11 - L)*(A22 - L) - A12*A21 = 0                           (9)
C...
C...  OR
C...
C...     L**2 - (A11 + A22)*L + (A11*A22 - A12*A21) = 0             (10)
C...
C...  SOLUTION BY THE QUADRATIC EQUATION GIVES THE TWO EIGENVALUES, L1
C...  AND L2.  SINCE EQUATIONS (1) AND (2) ARE LINEAR, THE GENERAL
C...  SOLUTION CONSISTS OF A LINEAR COMBINATION OF THE TWO EXPONENTIALS
C...  EXP(L1*T) AND EXP(L2*T), I.E.,
C...
C...     Y1 = C11*EXP(L1*T) + C12*EXP(L2*T)                         (11)
C...
C...     Y2 = C21*EXP(L1*T) + C22*EXP(L2*T)                         (12)
C...
C...  WHERE C11, C12, C21 AND C22 ARE CONSTANTS TO BE DETERMINED (NOTE
C...  THAT (C11,C21) IS THE EIGENVECTOR FOR EIGENVALUE L1, AND (C12,C22)
C...  IS THE EIGENVECTOR FOR L2).  APPLYING INITIAL CONDITIONS (3) AND
C...  (4), PLUS EITHER EQUATION (7) OR (8) FOR L1 AND L2, GIVES THE
C...  RESULT
C...
C...     Y1 = -EXP(L1*T) + EXP(L2*T)                                (13)
C...
C...     Y2 =  EXP(L1*T) + EXP(L2*T)                                (14)
C...
C...  FOR THE SPECIAL CASE A11 = A22 = -A, A12 = A21 = B, L1 = -(A + B),
C...  L2 = -(A - B).  THE FOLLOWING PARTICULAR VALUES OF A AND B ILLUS-
C...  TRATE VARYING DEGREES OF STIFFNESS
C...
C...     CASE 1   A = 50.5        L1 = -100        L1/L2 = 100
C...              B = 49.5        L2 =   -1        (NONSTIFF)
C...
C...     CASE 2   A = 500.5       L1 = -1000       L1/L2 = 1000
C...              B = 499.5       L2 =    -1    (MODERATELY STIFF)
C...
C...     CASE 3   A = 500000.5    L1 = -1000000    L1/L2 = 1000000
C...              B = 499999.5    L2 =       -1     (VERY STIFF)
C...
C...  THE THIRD CASE PROVIDES A STRINGENT TEST OF A NUMERICAL INTEGRATOR
C...  FOR STIFF DIFFERENTIAL EQUATIONS.  FOR EXAMPLE, IF AN EXPLICIT
C...  INTEGRATOR IS USED IN CASE THREE, THE MAXIMUM STEP SIZE DETERMINED
C...  BY THE STABILITY OF THE CALCULATION WILL BE GIVEN APPROXIMATELY BY
C...
C...      ABS(HMAX*L1) LT 2
C...
C...  AND WITH L1 = -1000000, HMAX LT 2/1000000 = 0.000002.  SINCE THE
C...  TIME FOR THE PROBLEM IS DETERMINED BY L2, THROUGH THE EXPONENTIAL
C...  EXP(L2*T) IN EQUATIONS (13) AND (14), THE NUMERICAL SOLUTION
C...  MUST BE COMPUTED APPROXIMATELY OVER THE INTERVAL O LE T LE 10.
C...  THUS, TO COMPUTE THE SOLUTION TO T = 10, 10/0.000002 = 5,000,000
C...  STEPS IN T WILL BE REQUIRED BY THE EXPLICIT INTEGRATOR.  AN
C...  IMPLICIT INTEGRATOR, HOWEVER, WHICH IS NOT LIMITED BY STABILITY
C...  SHOULD PRODUCE A SOLUTION OVER THE INTERVAL 0 LE T LE 10 WITH AT
C...  MOST A FEW HUNDRED STEPS (THE NUMBER OF STEPS  WILL BE DETERMINED
C...  BY THE REQUIRED ACCURACY AND NOT BY THE STABILITY OF THE CALCU-
C...  LATION).  THIS VASTLY IMPROVED EFFICIENCY OF IMPLICIT INTEGRATORS
C...  FOR THIS STIFF PROBLEM IS DEMONSTRATED BY THE FOLLOWING CODING
C...
C...  BECAUSE OF THE STIFFNESS OF THIS PROBLEM, DOUBLE PRECISION CAL-
C...  CULATIONS ARE ADVISED
      DOUBLE PRECISION
     1              T,
     2             Y1,    Y2,
     3          DY1DT, DY2DT,
     4              A,     B
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/    Y1,    Y2
     2      /F/ DY1DT, DY2DT
     3      /C/     A,     B
C...
C...  INITIAL CONDITIONS, EQUATIONS (3) AND (4)
      Y1=0.D+00
      Y2=2.D+00
C...
C...  THE FOLLOWING CONSTANTS IN THE ORDINARY DIFFERENTIAL EQUATIONS
C...  (ODES) CORRESPOND TO A STIFFNESS RATIO OF ONE MILLION (CASE 3
C...  ABOVE)
      A=500000.5D+00
      B=499999.5D+00
      RETURN
      END
      SUBROUTINE DERV
      DOUBLE PRECISION
     1              T,
     2             Y1,    Y2,
     3          DY1DT, DY2DT,
     4              A,     B
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/    Y1,    Y2
     2      /F/ DY1DT, DY2DT
     3      /C/     A,     B
C...
C...  ODES, EQUATIONS (1) AND (2)
      DY1DT = -A*Y1 + B*Y2
      DY2DT =  B*Y1 - A*Y2
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      DOUBLE PRECISION
     1              T,
     2             Y1,    Y2,
     3          DY1DT, DY2DT,
     4              A,     B,   AMB,  EXP1,  EXP2,   Y1E,   Y2E,
     5            Y1D,   Y2D
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/    Y1,    Y2
     2      /F/ DY1DT, DY2DT
     3      /C/     A,     B
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(T.EQ.0.D+00)WRITE(NO,2)
2     FORMAT(1H ,//,9X,1HT,
     1 13X,2HY1,12X,3HY1E,12X,3HY1D,/,
     2 23X,2HY2,12X,3HY2E,12X,3HY2D)
C...
C...  EVALUATE THE EXACT (ANALYTICAL) SOLUTION, EQUATIONS (11) AND (12)
C...  WITH L1 = -(A + B), L2 = -(A - B), FOR COMPARISON WITH THE NUMERI-
C...  CAL SOLUTION
      AMB=A-B
      IF(AMB.EQ.0.D+00)AMB=1.D+00
      EXP1=DEXP(-AMB*T)
      EXP2=0.D+00
      IF(((A+B)*T).LT.50.D+00)EXP2=DEXP(-(A+B)*T)
      Y1E=EXP1-EXP2
      Y2E=EXP1+EXP2
C...
C...  COMPUTE THE DIFFERENCES BETWEEN THE NUMERICAL AND EXACT SOLUTIONS
      Y1D=Y1-Y1E
      Y2D=Y2-Y2E
C...
C...  PRINT THE INDEPENDENT VARIABLE, NUMERICAL AND EXACT SOLUTIONS, AND
C...  THE DIFFERENCES BETWEEN THE TWO SOLUTIONS
      WRITE(NO,1)T,Y1,Y1E,Y1D,Y2,Y2E,Y2D
1     FORMAT(1H ,/,D10.3,3D15.4,/,10X,3D15.4)
      RETURN
      END
