https://github.com/cran/fOptions
Revision e4884f7b9a4e5e51f1d92f988d3aa6a5e2531ae1 authored by Diethelm Wuertz on 08 August 1977, 00:00:00 UTC, committed by Gabor Csardi on 08 August 1977, 00:00:00 UTC
1 parent cdde1d3
Tip revision: e4884f7b9a4e5e51f1d92f988d3aa6a5e2531ae1 authored by Diethelm Wuertz on 08 August 1977, 00:00:00 UTC
version 201.10059
version 201.10059
Tip revision: e4884f7
E5-EBMAsianOptions.f
C ALGORITHM 540R (REMARK ON ALG.540), COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 343-344.
C
C USED FOR:
C PROGRAM RUNPDE
C FOR TESTING AND DEBUGGING UNDER FORTRAN
C CALL PDETEST()
C END
C MODEL 1: VECER's PDE
C MODEL 2: ZHANG's PDE
C ------------------------------------------------------------------------------
SUBROUTINE PDETEST()
CC NOT USED BY R
CC FOR TESTING AND DEBUGGING UNDER FORTRAN
IMPLICIT REAL*8 (A-H, O-Z)
PARAMETER(MNP=10)
DIMENSION PRICE(MNP+1), XBYS(MNP+1)
C
PARAMETER (MMF=12, MMX=1000)
PARAMETER (MNPDE=1, MKORD=4, MNINT=MMX, MNCC=2, MMAXDER=5)
C WORKING ARRAYS:
DIMENSION WORK
* (MKORD+MNPDE*(4+9*MNPDE)+(MKORD+(MNINT-1)*(MKORD-MNCC))*
* (3*MKORD+2+MNPDE*(3*(MKORD-1)*MNPDE+MMAXDER+4)))
DIMENSION IWORK((MNPDE+1)*(MKORD+(MNINT-1)*(MKORD-MNCC)))
DIMENSION XBKPT(MNINT+1)
C PDE PARAMETERS:
NP = MNP
MF = MMF
MX = MMX
NPDE = MNPDE
KORD = MKORD
NINT = MNINT
NCC = MNCC
MAXDER = MMAXDER
C OPTION SETTINGS:
SIGMA = 0.30D0
TIME = 1.00D0
RR = 0.09D0
XS = 100.00D0
XSMIN = 90.00D0
XSMAX = 110.00D0
SS = 100.00D0
DELTA = (XSMAX-XSMIN)/NP
DO I = 1, NP+1
XBYS(I) = (XSMIN +(I-1)*DELTA)/XS
ENDDO
C SET TIME POINTS:
C T0 = INITIAL VALUE OF T, THE INDEPENDENT VARIABLE
C TOUT = VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT
C DT = INITIAL STEP SIZE IN T
C EPS = RELATIVE TIME ERROR BOUND
T0 = 0.0D0
TOUT = 1.0D0
EPS = 1.0D-08
DT = 1.0D-10
C FURTHER PARAMETERS:
C NINT=1000 - NUMBER OF SUBINTERVALS (XLEFT,XRIGHT) IS TO BE DIVIDED
C KORD=4 - ORDER OF THE PIECEWISE POLYNOMIAL SPACE TO BE USED
C NCC=2 - NUMBER OF CONTINUITY CONDITIONS TO BE IMPOSED
C MF=12 - METHOD FLAG
C ADAMS METHODS - GENERALIZATIONS OF CRANK-NICOLSON AND
C CHORD METHOD WITH FINITE DIFFERENCES JACOBIAN
C INDEX - INTEGER USED ON INPUT TO INDICATE TYPE OF CALL
C WORK - WORKING ARRAY
C IWORK - SIZE OF WORKING ARRAY
C ASIAN CALL (1) AND PUT(2) VALUE:
Z = -1
DO IP = 1, 2
Z = -Z
C PDE PARAMETERS:
MODSEL = 1
SIGMAT = SIGMA * DSQRT(TIME)
RRT = RR*TIME
XM = 5.0D0 * SIGMAT
WRITE (*,*)
WRITE (*,*) " PDE - ASIAN OPTION SETTINGS"
WRITE (*,*) " MF KORD NCC : ", MF, KORD, NCC
WRITE (*,*) " SIGMA*TIME : ", SIGMAT
WRITE (*,*) " R*TIME : ", RRT
WRITE (*,*) " XM : ", XM
WRITE (*,*) " (XMIN,XMAX)/S : ", XSMIN/SS, XSMAX/SS
CALL ASIANVAL(
& Z, SS, XS, XSMIN, XSMAX, TIME, RR, SIGMA,
& T0, TOUT, EPS, DT, PRICE, NP, MODSEL,
& MF, NPDE, KORD, MX, NCC, MAXDER,
& XBYS, XBKPT, WORK, IWORK)
C OUTPUT U - NUMERICAL SOLUTION:
WRITE (*,*) " XLEFT XRIGHT : ", XBKPT(1), XBKPT(NINT+1)
WRITE (*,*) " EPS DT MX : ", EPS, DT, MX
WRITE (*,*) " ERROR CODE: : ", INDEX
WRITE(*,*)
WRITE(*,*) " U - NUMERICAL SOLUTION FOR DIFF STRIKES:"
WRITE(*,*) " X XI PRICE "
DO I = 1, NP+1
XI = XBYS(I)*EXP(-RRT) - (1.0-EXP(-RRT))/RRT
WRITE(*,9) XS*XBYS(I), SS*XI, SS*PRICE(I), SS*(PRICE(I)-XI)
ENDDO
ENDDO
9 FORMAT(F10.3, 4F15.7)
RETURN
END
C ------------------------------------------------------------------------------
SUBROUTINE ASIANVAL(
& ZZ, SS1, XS1, XSMIN, XSMAX, TIME1, RR1, SIGMA1,
& T0, TOUT, EPS, DT, PRICEBYS, NP, MODSEL,
& MF1, NPDE1, KORD1, MX1, NCC1, MAXDER1,
& XBYS, XBKPT, WORK, IWORK)
IMPLICIT REAL*8 (A-H, O-Z)
PARAMETER(MKORD=4, MDERV=0)
DIMENSION WORK
* (KORD1+NPDE1*(4+9*NPDE1)+(KORD1+(MX1-1)*(KORD1-NCC1))*
* (3*KORD1+2+NPDE1*(3*(KORD1-1)*NPDE1+MAXDER1+4)))
DIMENSION IWORK((NPDE1+1)*(KORD1+(MX1-1)*(KORD1-NCC1)))
DIMENSION XBKPT(MX1+1)
DIMENSION USOL(1,1,MDERV+1), SCRTCH(MKORD*(MDERV+1))
DIMENSION XBYS(NP), PRICEBYS(NP)
COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
COMMON /GEAR0/ HUSED, NQUSED, NS, NF, NJ
COMMON /GEAR1/ T,DTC,DTMN,DTMX,EPSC,UROUND,N,MFC,KFLAG,JSTART
COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W
COMMON /OPTION/ NOGAUS,MAXDER
COMMON /ISTART/ IW1, IW2, IW3, IDUM(15)
COMMON /PARAMS/ PI
COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL
COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR
C FOR COMMON BLOCKS:
SIGMA = SIGMA1
TIME = TIME1
RR = RR1
XS = XS1
SS = SS1
C FOR COMMON BLOCKS:
MF = MF1
NPDE = NPDE1
KORD = KORD1
MX = MX1
NCC = NCC1
MAXDER = MAXDER1
NINT = MX1
MODEL = MODSEL
PI = 4.0D0 * DATAN(1.0D0)
C CALCULATE FOR BOTH, FOR A CALL Z=+1 OR FOR A PUT Z=-1:
Z = ZZ
C WORKSPACE SETTINGS:
IWORK(1) = KORD+NPDE*(4+9*NPDE)+(KORD+(MX-1)*
* (KORD-NCC))*(3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4))
IWORK(2) = (NPDE+1)*(KORD+(NINT-1)*(KORD-NCC))
DO I = 1, IWORK(1)
WORK(I)=0.0
ENDDO
C OPTION SETTINGS:
SIGMAT = SIGMA * DSQRT(TIME)
RRT = RR*TIME
XM = 5.0D0 * SIGMAT
XL = -XM
XR = +XM
ETA = (SIGMA**2)*(TIME**3)/6.0D0
C SET SPACE POINTS:
NX = MX
DX = 2.0D0 * XM / NX
DO I = 1, NX + 1
XBKPT(I) = -XM + (I-1)*DX
ENDDO
C SOLVE PDE:
INDEX = 1
CALL PDECOL(T0, TOUT, DT, XBKPT, EPS,
& NX, KORD, NCC, NPDE, MF, INDEX, WORK, IWORK)
C OUTPUT U - NUMERICAL SOLUTION:
DO I = 1, NP+1
XI = XBYS(I)*DEXP(-RRT) - (1.0D0-DEXP(-RRT))/RRT
CALL VALUES(XI, USOL, SCRTCH, 1, 1, 1, 0, WORK)
PRICEBYS(I) = USOL(1,1,1)
ENDDO
RETURN
END
C ------------------------------------------------------------------------------
SUBROUTINE F(T, X, U, UX, UXX, FVAL, NPDE)
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE)
COMMON /GEAR0/ HUSED, NQUSED, NS, NF, NJ
COMMON /PARAMS/ PI
COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL
COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR
IF (MODEL.EQ.1) THEN
FR = (1.0D0-DEXP(-RR*T))/RRT
FVAL(1) = (0.5D0*SIGMAT*SIGMAT) * ((X+FR)**2) * UXX(1)
ENDIF
IF (MODEL.EQ.2) THEN
RT = (1.0D0-DEXP(-RR*T))/RR
PF = (X*SIGMA*SIGMA)/(4.0D0*DSQRT(PI*ETA))
FVAL(1) = (0.5D0*SIGMA*SIGMA) * ((X+RT)**2) * UXX(1)
FVAL(1) = FVAL(1) + PF * DEXP(-0.25D0*X*X/ETA) * (X+2.0D0*RT)
ENDIF
RETURN
END
C ------------------------------------------------------------------------------
SUBROUTINE BNDRY(T, X, U, UX, DBDU, DBDUX, DZDT, NPDE)
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE)
DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE)
COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL
COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR
C LEFT/RIGHT BOUNDARY MODEL 1:
IF (MODEL.EQ.1) THEN
IF (X.LE.-XM) THEN
DBDU (1,1) = (-Z*X + DABS(X) ) / 2.0D0
DBDUX(1,1) = 0.0D0
DZDT (1) = 0.0D0
RETURN
ENDIF
IF (X.LE.XM) THEN
DBDU (1,1) = (-Z*X + DABS(X) ) / 2.0D0
DBDUX(1,1) = 0.0D0
DZDT (1) = 0.0D0
RETURN
ENDIF
ENDIF
C LEFT/RIGHT BOUNDARY MODEL 2:
IF (MODEL.EQ.2) THEN
EPS = 1.0D-20
IF (X.LE.XL ) THEN
DBDU (1,1) = EPS
DBDUX(1,1) = 0.0D0
DZDT (1) = 0.0D0
RETURN
ENDIF
IF (X.GE.XR ) THEN
DBDU (1,1) = EPS
DBDUX(1,1) = 0.0D0
DZDT (1) = 0.0D0
RETURN
ENDIF
ENDIF
RETURN
END
C ------------------------------------------------------------------------------
SUBROUTINE UINIT(X, U, NPDE)
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION U(NPDE)
COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL
COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR
C NOTE : Z=+1 FOR A CALL AND Z-1 FOR A PUT
IF (MODEL.EQ.1) THEN
U(1) = ( (-Z*X) + DABS(-X) ) / 2.0D0
ENDIF
IF (MODEL.EQ.2) THEN
U(1) = 0.0D0
ENDIF
RETURN
END
C ------------------------------------------------------------------------------
SUBROUTINE DERIVF(T, X, U, UX, UXX, DFDU, DFDUX, DFDUXX, NPDE)
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION U(NPDE), UX(NPDE), UXX(NPDE)
DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE)
COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL
COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR
C
C IF THE USER DESIRES TO USE THE MF = 11 OR 21 OPTION IN ORDER TO SAVE
C ABOUT 10-20 PERCENT IN EXECUTION TIME (SEE BELOW), THEN THE USER MUST
C PROVIDE THE FOLLOWING SUBROUTINE WHICH PROVIDES INFORMATION ABOUT THE
C DERIVATIVES OF THE FUNCTION F ABOVE. THIS PROVIDES FOR MORE EFFICIENT
C JACOBIAN MATRIX GENERATION. ON MOST COMPUTER SYSTEMS, THE USER WILL
C BE REQUIRED TO SUPPLY THIS SUBROUTINE AS A DUMMY SUBROUTINE IF THE
C OPTIONS MF = 12 OR 22 ARE USED (SEE BELOW).
C
C THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES T, X, U, UX,
C AND UXX, AND THE USER SHOULD CONSTRUCT THIS ROUTINE TO PROVIDE
C THE FOLLOWING CORRESPONDING VALUES OF THE OUTPUT ARRAYS
C DFDU, DFDUX, AND DFDUXX FOR K,J = 1 TO NPDE...
C DFDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C PDE DEFINING FUNCTION F WITH RESPECT TO THE
C VARIABLE U(J).
C DFDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C PDE DEFINING FUNCTION F WITH RESPECT TO THE
C VARIABLE UX(J).
C DFDUXX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C PDE DEFINING FUNCTION F WITH RESPECT TO THE
C VARIABLE UXX(J).
C NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT
C VALUE.
PI = 4.0 * DATAN(1.0D0)
IF (MODEL.EQ.1) THEN
RT = (1.0D0-EXP(-RRT*T))/RRT
DFDU(1,1) = 0.0D0
DFDUX(1,1) = 0.0D0
DFDUXX(1,1) = (SIGMAT**2) * ( X + RT )
ENDIF
IF (MODEL.EQ.1) THEN
RT = (1.0D0-DEXP(-RR*T))/RR
F1 = (X*SIGMA*SIGMA)/(4.0D0*DSQRT(PI*ETA))
F2 = DEXP(-0.25D0*X*X/ETA)
F3 = (X+2.0D0*RT)
DF1 = F1 / X
DF2 = -2.0D0 * X * F2 / (4.0D0*ETA)
DF3 = 1.0D0
DFDUXX(1,1) = (SIGMA**2) * ( X + RT )
DFDUX(1,1) = 0.0D0
DFDU(1,1) = DF1*F2*F3 + F1*DF2*F3 + F1*F2*DF3
ENDIF
RETURN
END
C ##############################################################################
C ALGORITHM 540R (REMARK ON ALG.540), COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 343-344.
C
C
SUBROUTINE PDECOL(T0, TOUT, DT, XBKPT, EPS, NINT, KORD,
* NCC, NPDE, MF, INDEX, WORK, IWORK)
IMPLICIT REAL*8 (A-H, O-Z)
C
C
C-------------------------------------------------------------------------------
C
C THIS IS THE MARCH 24, 1978 VERSION OF PDECOL.
C
C THIS PACKAGE WAS CONSTRUCTED SO AS TO CONFORM TO AS MANY ANSI-FORTRAN
C RULES AS WAS CONVENIENTLY POSSIBLE. THE FORTRAN USED VIOLATES ANSI
C STANDARDS IN THE TWO WAYS LISTED BELOW....
C
C 1. SUBSCRIPTS OF THE GENERAL FORM C*V1 + V2 + V3 ARE USED
C (POSSIBLY IN A PERMUTED ORDER), WHERE C IS AN INTEGER CONSTANT
C AND V1, V2, AND V3 ARE INTEGER VARIABLES.
C
C 2. ARRAY NAMES APPEAR SINGLY IN DATA STATEMENTS IN THE ROUTINES
C BSPLVN AND COSET.
C
C MACHINE DEPENDENT FEATURES......
C
C THIS VERSION OF PDECOL WAS DESIGNED FOR USE ON CDC MACHINES WITH
C A WORD LENGTH OF 60 BITS. WE DO NOT RECOMMEND THE USE OF PDECOL WITH
C WORD LENGTHS OF LESS THAN 48 BITS. THE MOST IMPORTANT MACHINE
C AND WORD LENGTH DEPENDENT CONSTANTS ARE DEFINED IN THE BLOCK DATA
C AND IN SUBROUTINES COLPNT AND COSET. THE USER SHOULD CHECK THESE
C CAREFULLY FOR APPROPRIATENESS FOR HIS LOCAL SITUATION. THE FORTRAN
C FUNCTIONS USED BY EACH ROUTINE ARE LISTED IN THE COMMENTS TO
C FACILITATE CONVERSION TO DOUBLE PRECISION.
C
C-------------------------------------------------------------------------------
C
C PDECOL IS THE DRIVER ROUTINE FOR A SOPHISTICATED PACKAGE OF
C SUBROUTINES WHICH IS DESIGNED TO SOLVE THE GENERAL SYSTEM OF
C NPDE NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS OF AT MOST SECOND
C ORDER ON THE INTERVAL (XLEFT,XRIGHT) FOR T .GT. T0 WHICH IS OF THE
C FORM....
C
C DU/DT = F( T, X, U, UX, UXX )
C
C WHERE
C
C U = ( U(1), U(2), ... , U(NPDE) )
C UX = ( UX(1), UX(2), ... , UX(NPDE) )
C UXX = (UXX(1),UXX(2), ... ,UXX(NPDE) ) .
C
C EACH U(K) IS A FUNCTION OF THE SCALAR QUANTITIES T AND X.
C UX(K) REPRESENTS THE FIRST PARTIAL DERIVATIVE OF U(K) WITH RESPECT
C TO THE VARIABLE X, UXX(K) REPRESENTS THE SECOND PARTIAL DERIVATIVE
C OF U(K) WITH RESPECT TO THE VARIABLE X, AND DU/DT IS THE VECTOR OF
C PARTIAL DERIVATIVES OF U WITH RESPECT TO THE TIME VARIABLE T.
C F REPRESENTS AN ARBITRARY VECTOR VALUED FUNCTION WHOSE NPDE
C COMPONENTS DEFINE THE RESPECTIVE PARTIAL DIFFERENTIAL EQUATIONS OF
C THE PDE SYSTEM. SEE SUBROUTINE F DESCRIPTION BELOW.
C
C BOUNDARY CONDITIONS
C
C DEPENDING ON THE TYPE OF PDE(S), 0, 1, OR 2 BOUNDARY CONDITIONS
C ARE REQUIRED FOR EACH PDE IN THE SYSTEM. THESE ARE IMPOSED AT XLEFT
C AND/OR XRIGHT AND EACH MUST BE OF THE FORM....
C
C B(U,UX) = Z(T)
C
C WHERE B AND Z ARE ARBITRARY VECTOR VALUED FUNCTIONS WITH
C NPDE COMPONENTS AND U, UX, AND T ARE AS ABOVE. THESE BOUNDARY
C CONDITIONS MUST BE CONSISTENT WITH THE INITIAL CONDITIONS WHICH ARE
C DESCRIBED NEXT.
C
C INITIAL CONDITIONS
C
C EACH SOLUTION COMPONENT U(K) IS ASSUMED TO BE A KNOWN (USER
C PROVIDED) FUNCTION OF X AT THE INITIAL TIME T = T0. THE
C INITIAL CONDITION FUNCTIONS MUST BE CONSISTENT WITH THE BOUNDARY
C CONDITIONS ABOVE, I.E. THE INITIAL CONDITION FUNCTIONS MUST
C SATISFY THE BOUNDARY CONDITIONS FOR T = T0. SEE SUBROUTINE UINIT
C DESCRIPTION BELOW.
C
C-------------------------------------------------------------------------------
C
C REQUIRED USER SUPPLIED SUBROUTINES
C
C THE USER IS REQUIRED TO CONSTRUCT THREE SUBPROGRAMS AND A MAIN
C PROGRAM WHICH DEFINE THE PDE PROBLEM WHOSE SOLUTION IS TO BE
C ATTEMPTED. THE THREE SUBPROGRAMS ARE...
C
C 1) SUBROUTINE F( T, X, U, UX, UXX, FVAL, NPDE )
C DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE)
C THIS ROUTINE DEFINES THE DESIRED PARTIAL DIFFERENTIAL
C EQUATIONS TO BE SOLVED. THE PACKAGE PROVIDES VALUES OF THE
C INPUT SCALARS T AND X AND INPUT ARRAYS (LENGTH NPDE) U, UX,
C AND UXX, AND THE USER MUST CONSTRUCT THIS ROUTINE TO COMPUTE
C THE OUTPUT ARRAY FVAL (LENGTH NPDE) WHICH CONTAINS THE
C CORRESPONDING VALUES OF THE RIGHT HAND SIDES OF THE DESIRED
C PARTIAL DIFFERENTIAL EQUATIONS, I.E.
C
C FVAL(K) = THE VALUE OF THE RIGHT HAND SIDE OF THE K-TH PDE IN
C THE PDE SYSTEM ABOVE, FOR K = 1 TO NPDE.
C
C THE INCOMING VALUE OF THE SCALAR QUANTITY X WILL BE A
C COLLOCATION POINT VALUE (SEE INITAL AND COLPNT) AND THE
C INCOMING VALUES IN THE ARRAYS U, UX AND UXX CORRESPOND TO THIS
C POINT X AND TIME T.
C RETURN
C END
C
C 2) SUBROUTINE BNDRY( T, X, U, UX, DBDU, DBDUX, DZDT, NPDE )
C DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE)
C DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE)
C THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH NEEDED
C INFORMATION ABOUT THE BOUNDARY CONDITION FUNCTIONS B AND Z
C ABOVE. THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES
C T, X, U, AND UX, AND THE USER IS TO DEFINE THE CORRESPONDING
C OUTPUT VALUES OF THE DERIVATIVES OF THE FUNCTIONS B AND Z
C WHERE....
C DBDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO
C THE J-TH VARIABLE U(J).
C DBDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO
C THE J-TH VARIABLE UX(J).
C DZDT(K) = DERIVATIVE OF THE K-TH COMPONENT OF THE VECTOR
C FUNCTION Z(T) ABOVE WITH RESPECT TO THE
C VARIABLE T.
C NOTE... THE INCOMING VALUE OF X WILL BE EITHER XLEFT OR XRIGHT.
C IF NO BOUNDARY CONDITION IS DESIRED FOR SAY THE K-TH PDE AT
C ONE OR BOTH OF THE ENDPOINTS XLEFT OR XRIGHT, THEN DBDU(K,K)
C AND DBDUX(K,K) SHOULD BOTH BE SET TO ZERO WHEN BNDRY IS
C CALLED FOR THAT POINT. WE REFER TO THIS AS A NULL BOUNDARY
C CONDITION. THIS ROUTINE CAN BE STRUCTURED AS FOLLOWS...
C THE COMMON BLOCK /ENDPT/ IS NOT A PART OF PDECOL AND
C MUST BE SUPPLIED AND DEFINED BY THE USER.
C COMMON /ENDPT/ XLEFT
C IF( X .NE. XLEFT ) GO TO 10
C HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J),
C AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE LEFT BOUNDARY POINT
C X = XLEFT.
C RETURN
C 10 CONTINUE
C HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J),
C AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE RIGHT BOUNDARY POINT
C X = XRIGHT.
C RETURN
C END
C
C 3) SUBROUTINE UINIT( X, U, NPDE )
C DIMENSION U(NPDE)
C THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH THE
C NEEDED INITIAL CONDITION FUNCTION VALUES. THE PACKAGE
C PROVIDES A VALUE OF THE INPUT VARIABLE X, AND THE USER IS TO
C DEFINE THE PROPER INITIAL VALUES (AT T = T0) FOR ALL OF THE
C PDE COMPONENTS, I.E.
C U(K) = DESIRED INITIAL VALUE OF PDE COMPONENT U(K) AT
C X AND T = T0 FOR K = 1 TO NPDE.
C NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT
C VALUE. THE INITIAL CONDITIONS AND BOUNDARY CONDITIONS
C MUST BE CONSISTENT (SEE ABOVE).
C RETURN
C END
C
C-------------------------------------------------------------------------------
C
C OPTIONAL USER SUPPLIED SUBROUTINE
C
C IF THE USER DESIRES TO USE THE MF = 11 OR 21 OPTION IN ORDER TO SAVE
C ABOUT 10-20 PERCENT IN EXECUTION TIME (SEE BELOW), THEN THE USER MUST
C PROVIDE THE FOLLOWING SUBROUTINE WHICH PROVIDES INFORMATION ABOUT THE
C DERIVATIVES OF THE FUNCTION F ABOVE. THIS PROVIDES FOR MORE EFFICIENT
C JACOBIAN MATRIX GENERATION. ON MOST COMPUTER SYSTEMS, THE USER WILL
C BE REQUIRED TO SUPPLY THIS SUBROUTINE AS A DUMMY SUBROUTINE IF THE
C OPTIONS MF = 12 OR 22 ARE USED (SEE BELOW).
C
C 1) SUBROUTINE DERIVF( T, X, U, UX, UXX, DFDU, DFDUX, DFDUXX, NPDE )
C DIMENSION U(NPDE), UX(NPDE), UXX(NPDE)
C DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE)
C THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES T, X, U, UX,
C AND UXX, AND THE USER SHOULD CONSTRUCT THIS ROUTINE TO PROVIDE
C THE FOLLOWING CORRESPONDING VALUES OF THE OUTPUT ARRAYS
C DFDU, DFDUX, AND DFDUXX FOR K,J = 1 TO NPDE...
C DFDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C PDE DEFINING FUNCTION F WITH RESPECT TO THE
C VARIABLE U(J).
C DFDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C PDE DEFINING FUNCTION F WITH RESPECT TO THE
C VARIABLE UX(J).
C DFDUXX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C PDE DEFINING FUNCTION F WITH RESPECT TO THE
C VARIABLE UXX(J).
C NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT
C VALUE.
C RETURN
C END
C
C-------------------------------------------------------------------------------
C
C METHODS USED
C
C THE PACKAGE PDECOL IS BASED ON THE METHOD OF LINES AND USES A
C FINITE ELEMENT COLLOCATION PROCEDURE (WITH PIECEWISE POLYNOMIALS
C AS THE TRIAL SPACE) FOR THE DISCRETIZATION OF THE SPATIAL VARIABLE
C X. THE COLLOCATION PROCEDURE REDUCES THE PDE SYSTEM TO A SEMI-
C DISCRETE SYSTEM WHICH THEN DEPENDS ONLY ON THE TIME VARIABLE T.
C THE TIME INTEGRATION IS THEN ACCOMPLISHED BY USE OF SLIGHTLY
C MODIFIED STANDARD TECHNIQUES (SEE REFS. 1,2).
C
C PIECEWISE POLYNOMIALS
C
C THE USER IS REQUIRED TO SELECT THE PIECEWISE POLYNOMIAL SPACE
C WHICH IS TO BE USED TO COMPUTE HIS APPROXIMATE SOLUTION. FIRST, THE
C ORDER, KORD, OF THE POLYNOMIALS TO BE USED MUST BE SPECIFIED
C (KORD = POLYNOMIAL DEGREE + 1). NEXT, THE NUMBER OF PIECES
C (INTERVALS), NINT, INTO WHICH THE SPATIAL DOMAIN (XLEFT,XRIGHT) IS
C TO BE DIVIDED, IS CHOSEN. THE NINT + 1 DISTINCT BREAKPOINTS OF
C THE DOMAIN MUST BE DEFINED AND SET INTO THE ARRAY XBKPT IN
C STRICTLY INCREASING ORDER, I.E.
C XLEFT=XBKPT(1) .LT. XBKPT(2) .LT. ... .LT. XBKPT(NINT+1)=XRIGHT.
C THE APPROXIMATE SOLUTION AT ANY TIME T WILL BE A POLYNOMIAL OF
C ORDER KORD OVER EACH SUBINTERVAL (XBKPT(I),XBKPT(I+1)). THE
C NUMBER OF CONTINUITY CONDITIONS, NCC, TO BE IMPOSED ACROSS ALL OF
C THE BREAKPOINTS IS THE LAST PIECE OF USER SUPPLIED DATA WHICH IS
C REQUIRED TO UNIQUELY DETERMINE THE DESIRED PIECEWISE POLYNOMIAL
C SPACE. FOR EXAMPLE, NCC = 2 WOULD REQUIRE THAT THE APPROXIMATE
C SOLUTION (MADE UP OF THE SEPARATE POLYNOMIAL PIECES) AND ITS FIRST
C SPATIAL DERIVATIVE BE CONTINUOUS AT THE BREAKPOINTS AND HENCE ON
C THE ENTIRE DOMAIN (XLEFT,XRIGHT). NCC = 3 WOULD REQUIRE THAT THE
C APPROXIMATE SOLUTION AND ITS FIRST AND SECOND SPATIAL DERIVATIVES
C BE CONTINUOUS AT THE BREAKPOINTS, ETC. THE DIMENSION OF THIS LINEAR
C SPACE IS KNOWN AND FINITE AND IS NCPTS = KORD*NINT - NCC*(NINT-1).
C THE WELL-KNOWN B-SPLINE BASIS (SEE REF. 3) FOR THIS SPACE IS USED
C BY PDECOL AND IT CONSISTS OF NCPTS KNOWN PIECEWISE POLYNOMIAL
C FUNCTIONS BF(I,X), FOR I=1 TO NCPTS, WHICH DO NOT DEPEND ON THE
C TIME VARIABLE T. WE WISH TO EMPHASIZE THAT THE PIECEWISE POLYNOMIAL
C SPACE USED IN PDECOL (WHICH IS SELECTED BY THE USER) WILL DETERMINE
C THE MAGNITUDE OF THE SPATIAL DISCRETIZATION ERRORS IN THE COMPUTED
C APPROXIMATE SOLUTION. THE PACKAGE HAS NO CONTROL OVER ERRORS
C INTRODUCED BY THE USERS CHOICE OF THIS SPACE. SEE INPUT PARAMETERS
C BELOW.
C
C COLLOCATION OVER PIECEWISE POLYNOMIALS
C
C THE BASIC ASSUMPTION MADE IS THAT THE APPROXIMATE SOLUTION
C SATISFIES
C NCPTS
C U(T,X) = SUM C(I,T) * BF(I,X)
C I=1
C
C WHERE THE UNKNOWN COEFFICIENTS C DEPEND ONLY ON THE TIME T AND
C THE KNOWN BASIS FUNCTIONS DEPEND ONLY ON X (WE HAVE ASSUMED THAT
C NPDE = 1 FOR CONVENIENCE). SO, AT ANY GIVEN TIME T THE APPROX-
C IMATE SOLUTION IS A PIECEWISE POLYNOMIAL IN THE USER CHOSEN SPACE.
C THE SEMI-DISCRETE EQUATIONS (ACTUALLY ORDINARY DIFFERENTIAL
C EQUATIONS) WHICH DETERMINE THE COEFFICIENTS C ARE OBTAINED BY
C REQUIRING THAT THE ABOVE APPROXIMATE U(T,X) SATISFY THE PDE AND
C BOUNDARY CONDITIONS EXACTLY AT A SET OF NCPTS COLLOCATION POINTS
C (SEE COLPNT). THUS, PDECOL ACTUALLY COMPUTES THE BASIS FUNCTION
C COEFFICIENTS RATHER THAN SPECIFIC APPROXIMATE SOLUTION VALUES.
C
C REFERENCES
C
C 1. MADSEN, N.K. AND R.F. SINCOVEC, PDECOL - COLLOCATION SOFTWARE
C FOR PARTIAL DIFFERENTIAL EQUATIONS, ACM-TOMS, VOL. , NO. ,
C 2. SINCOVEC, R.F. AND N.K. MADSEN, SOFTWARE FOR NONLINEAR PARTIAL
C DIFFERENTIAL EQUATIONS, ACM-TOMS, VOL. 1, NO. 3,
C SEPTEMBER 1975, PP. 232-260.
C 3. HINDMARSH, A.C., PRELIMINARY DOCUMENTATION OF GEARIB.. SOLUTION
C OF IMPLICIT SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS WITH
C BANDED JACOBIANS, LAWRENCE LIVERMORE LAB, UCID-30130, FEBRUARY
C 1976.
C 4. DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J.
C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472.
C
C-------------------------------------------------------------------------------
C
C USE OF PDECOL
C
C PDECOL IS CALLED ONCE FOR EACH DESIRED OUTPUT VALUE (TOUT) OF THE
C TIME T, AND IT IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR,
C STIFIB, WHICH ADVANCES THE TIME BY TAKING SINGLE STEPS UNTIL
C T .GE. TOUT. INTERPOLATION TO THE EXACT TIME TOUT IS THEN DONE.
C SEE TOUT BELOW.
C
C
C SUMMARY OF SUGGESTED INPUT VALUES
C
C IT IS OF COURSE IMPOSSIBLE TO SUGGEST INPUT PARAMETER VALUES WHICH
C ARE APPROPRIATE FOR ALL PROBLEMS. THE FOLLOWING SUGGESTIONS ARE TO
C BE USED ONLY IF YOU HAVE NO IDEA OF BETTER VALUES FOR YOUR PROBLEM.
C
C DT = 1.E-10
C XBKPT = CHOOSE NINT+1 EQUALLY SPACED VALUES SUCH THAT XBKPT(1) =
C XLEFT AND XBKPT(NINT+1) = XRIGHT.
C EPS = 1.E-4
C NINT = ENOUGH SO THAT ANY FINE STRUCTURE OF THE PROBLEM MAY BE
C RESOLVED.
C KORD = 4
C NCC = 2
C MF = 22
C INDEX = 1 (ON FIRST CALL ONLY, THEN 0 THEREAFTER).
C
C
C THE INPUT PARAMETERS ARE..
C T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE
C (USED ONLY ON FIRST CALL).
C TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. SINCE
C THE PACKAGE CHOOSES ITS OWN TIME STEP SIZES, THE
C INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT
C AND THE PACKAGE WILL INTERPOLATE TO T = TOUT.
C DT = THE INITIAL STEP SIZE IN T, IF INDEX = 1, OR, THE
C MAXIMUM STEP SIZE ALLOWED (MUST BE .GT. 0), IF INDEX = 3.
C USED FOR INPUT ONLY WHEN INDEX = 1 OR 3. SEE BELOW.
C XBKPT = THE ARRAY OF PIECEWISE POLYNOMIAL BREAKPOINTS.
C THE NINT+1 VALUES MUST BE STRICTLY INCREASING WITH
C XBKPT(1) = XLEFT AND XBKPT(NINT+1) = XRIGHT (USED ONLY
C ON FIRST CALL).
C EPS = THE RELATIVE TIME ERROR BOUND (USED ONLY ON THE
C FIRST CALL, UNLESS INDEX = 4). SINGLE STEP ERROR
C ESTIMATES DIVIDED BY CMAX(I) WILL BE KEPT LESS THAN
C EPS IN ROOT-MEAN-SQUARE NORM. THE VECTOR CMAX OF WEIGHTS
C IS COMPUTED IN PDECOL. INITIALLY CMAX(I) IS SET TO
C DABS(C(I)), WITH A DEFAULT VALUE OF 1 IF DABS(C(I)) .LT. 1.
C THEREAFTER, CMAX(I) IS THE LARGEST VALUE
C OF DABS(C(I)) SEEN SO FAR, OR THE INITIAL CMAX(I) IF
C THAT IS LARGER. TO ALTER EITHER OF THESE, CHANGE THE
C APPROPRIATE STATEMENTS IN THE DO-LOOPS ENDING AT
C STATEMENTS 50 AND 130 BELOW. THE USER SHOULD EXERCISE
C SOME DISCRETION IN CHOOSING EPS. IN GENERAL, THE
C OVERALL RUNNING TIME FOR A PROBLEM WILL BE GREATER IF
C EPS IS CHOSEN SMALLER. THERE IS USUALLY LITTLE REASON TO
C CHOOSE EPS MUCH SMALLER THAN THE ERRORS WHICH ARE BEING
C INTRODUCED BY THE USERS CHOICE OF THE POLYNOMIAL SPACE.
C SEE RELATED COMMENTS CONCERNING CMAX BELOW STATEMENT 40.
C NINT = THE NUMBER OF SUBINTERVALS INTO WHICH THE SPATIAL DOMAIN
C (XLEFT,XRIGHT) IS TO BE DIVIDED (MUST BE .GE. 1)
C (USED ONLY ON FIRST CALL).
C KORD = THE ORDER OF THE PIECEWISE POLYNOMIAL SPACE TO BE USED.
C ITS VALUE MUST BE GREATER THAN 2 AND LESS THAN 21. FOR
C FIRST ATTEMPTS WE RECOMMEND KORD = 4. IF THE SOLUTION
C IS SMOOTH AND MUCH ACCURACY IS DESIRED, HIGHER VALUES
C MAY PROVE TO BE MORE EFFICIENT. WE HAVE SELDOM USED
C VALUES OF KORD IN EXCESS OF 8 OR 9, THOUGH THEY ARE
C AVAILABLE FOR USE IN PDECOL (USED ONLY ON FIRST CALL).
C NCC = THE NUMBER OF CONTINUITY CONDITIONS TO BE IMPOSED ON THE
C APPROXIMATE SOLUTION AT THE BREAKPOINTS IN XBKPT.
C NCC MUST BE GREATER THAN 1 AND LESS THAN KORD. WE
C RECOMMEND THE USE OF NCC = 2
C SINCE THEORY PREDICTS THAT DRAMATICALLY MORE
C ACCURATE RESULTS CAN OFTEN BE OBTAINED USING THIS CHOICE
C (USED ONLY ON FIRST CALL).
C NPDE = THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS IN THE SYSTEM
C TO BE SOLVED (USED ONLY ON FIRST CALL).
C MF = THE METHOD FLAG (USED ONLY ON FIRST CALL, UNLESS
C INDEX = 4). ALLOWED VALUES ARE 11, 12, 21, 22.
C FOR FIRST ATTEMPTS WE RECOMMEND THE USE OF MF = 22.
C MF HAS TWO DECIMAL DIGITS, METH AND MITER
C (MF = 10*METH + MITER).
C METH IS THE BASIC METHOD INDICATOR..
C METH = 1 MEANS THE ADAMS METHODS (GENERALIZATIONS OF
C CRANK-NICOLSON).
C METH = 2 MEANS THE BACKWARD DIFFERENTIATION
C FORMULAS (BDF), OR STIFF METHODS OF GEAR.
C MITER IS THE ITERATION METHOD INDICATOR
C AND DETERMINES HOW THE JACOBIAN MATRIX IS
C TO BE COMPUTED..
C MITER = 1 MEANS CHORD METHOD WITH ANALYTIC JACOBIAN.
C FOR THIS USER SUPPLIES SUBROUTINE DERIVF.
C SEE DESCRIPTION ABOVE.
C MITER = 2 MEANS CHORD METHOD WITH JACOBIAN CALCULATED
C INTERNALLY BY FINITE DIFFERENCES. SEE
C SUBROUTINES PSETIB AND DIFFF.
C INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL,
C WITH THE FOLLOWING VALUES AND MEANINGS..
C 1 THIS IS THE FIRST CALL FOR THIS PROBLEM.
C 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM,
C AND INTEGRATION IS TO CONTINUE.
C 2 SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT
C EXACTLY (NO INTERPOLATION IS DONE). SEE NOTE
C BELOW. ASSUMES TOUT .GE. THE CURRENT T.
C IF TOUT IS .LT. THE CURRENT TIME, THEN TOUT IS
C RESET TO THE CURRENT TIME AND CONTROL IS
C RETURNED TO THE USER. A CALL TO VALUES WILL
C PRODUCE ANSWERS FOR THE NEW VALUE OF TOUT.
C 3 SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING
C PROGRAM AFTER ONE STEP. TOUT IS IGNORED AND
C DT MUST BE SET .GT. 0 TO A MAXIMUM ALLOWED
C DT VALUE. SEE ABOVE.
C 4 THIS IS NOT THE FIRST CALL FOR THE PROBLEM,
C AND THE USER HAS RESET EPS AND/OR MF.
C SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0,
C IT NEED NOT BE RESET FOR NORMAL CONTINUATION.
C
C NOTE.. THE PACKAGE MUST HAVE TAKEN AT LEAST ONE SUCCESSFUL TIME
C STEP BEFORE A CALL WITH INDEX = 2 OR 4 IS ALLOWED.
C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL
C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN.
C ALL OTHER PARAMETERS WILL BE Y FOR THE NEXT CALL.
C A CHANGE OF PARAMETERS WITH INDEX = 4 CAN BE MADE AFTER
C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN PROVIDED AT LEAST
C ONE SUCCESSFUL TIME STEP HAS BEEN MADE.
C
C WORK = FLOATING POINT WORKING ARRAY FOR PDECOL. WE RECOMMEND
C THAT IT BE INITIALIZED TO ZERO BEFORE THE FIRST CALL
C TO PDECOL. ITS TOTAL LENGTH MUST BE AT LEAST
C
C KORD + 4*NPDE + 9*NPDE**2 + NCPTS*(3*KORD + 2) +
C NPDE*NCPTS*(3*ML + MAXDER + 7)
C
C WHERE ML AND MAXDER ARE DEFINED BELOW (SEE STORAGE
C ALLOCATION).
C
C IWORK = INTEGER WORKING ARRAY FOR PDECOL. THE FIRST TWO
C LOCATIONS MUST BE DEFINED AS FOLLOWS...
C IWORK(1) = LENGTH OF USERS ARRAY WORK
C IWORK(2) = LENGTH OF USERS ARRAY IWORK
C THE TOTAL LENGTH OF IWORK MUST BE AT LEAST
C NCPTS*(NPDE + 1).
C OUTPUT
C
C THE SOLUTION VALUES ARE NOT RETURNED DIRECTLY TO THE USER BY PDECOL.
C THE METHODS USED IN PDECOL COMPUTE BASIS FUNCTION COEFFICIENTS, SO
C THE USER (AFTER A RETURN FROM PDECOL) MUST CALL THE PACKAGE ROUTINE
C VALUES TO OBTAIN HIS APPROXIMATE SOLUTION VALUES AT ANY DESIRED SPACE
C POINTS X AT THE TIME T = TOUT. SEE THE COMMENTS IN SUBROUTINE VALUES
C FOR DETAILS ON HOW TO PROPERLY MAKE THE CALL.
C
C EXECUTION ERROR MESSAGES WILL BE PRINTED BY PDECOL ON LOGICAL UNIT
C LOUT WHICH IS THE ONLY VARIABLE IN THE COMMON BLOCK /IOUNIT/. A
C DEFAULT OF LOUT = 3 IS SET IN THE BLOCK DATA.
C
C THE COMMON BLOCK /GEAR0/ CONTAINS THE VARIABLES DTUSED, NQUSED,
C NSTEP, NFE, AND NJE AND CAN BE ACCESSED EXTERNALLY BY THE USER IF
C DESIRED. RESPECTIVELY, IT CONTAINS THE STEP SIZE LAST USED (SUCCESS-
C FULLY), THE ORDER LAST USED (SUCCESSFULLY), THE NUMBER OF STEPS TAKEN
C SO FAR, THE NUMBER OF RESIDUAL EVALUATIONS (RES CALLS) SO FAR,
C AND THE NUMBER OF MATRIX EVALUATIONS (PSETIB CALLS) SO FAR.
C DIFFUN CALLS ARE COUNTED IN WITH RESIDUAL EVALUATIONS.
C
C THE OUTPUT PARAMETERS ARE..
C DT = THE STEP SIZE USED LAST, WHETHER SUCCESSFULLY OR NOT.
C TOUT = THE OUTPUT VALUE OF T. IF INTEGRATION WAS SUCCESSFUL,
C AND THE INPUT VALUE OF INDEX WAS NOT 3, TOUT IS
C UNCHANGED FROM ITS INPUT VALUE. OTHERWISE, TOUT
C IS THE CURRENT VALUE OF T TO WHICH THE INTEGRATION
C HAS BEEN COMPLETED.
C INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS,
C WITH THE FOLLOWING VALUES AND MEANINGS..
C 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND.
C -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE
C ERROR TEST EVEN AFTER REDUCING DT BY A FACTOR OF
C 1.E10 FROM ITS INITIAL VALUE.
C -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS
C HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY
C A TEST ON EPS. TOO MUCH ACCURACY HAS BEEN REQUESTED.
C -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE
C CORRECTOR CONVERGENCE EVEN AFTER REDUCING DT BY A
C FACTOR OF 1.E10 FROM ITS INITIAL VALUE.
C -4 SINGULAR MATRIX ENCOUNTERED. PROBABLY DUE TO STORAGE
C OVERWRITES.
C -5 INDEX WAS 4 ON INPUT, BUT THE DESIRED CHANGES OF
C PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT
C WAS NOT BEYOND T. INTERPOLATION TO T = TOUT WAS
C PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN,
C SIMPLY CALL AGAIN WITH INDEX = 4 AND A NEW TOUT.
C -6 ILLEGAL INDEX VALUE.
C -7 ILLEGAL EPS VALUE.
C -8 AN ATTEMPT TO INTEGRATE IN THE WRONG DIRECTION. THE
C SIGN OF DT IS WRONG RELATIVE TO T0 AND TOUT.
C -9 DT .EQ. 0.0.
C -10 ILLEGAL NINT VALUE.
C -11 ILLEGAL KORD VALUE.
C -12 ILLEGAL NCC VALUE.
C -13 ILLEGAL NPDE VALUE.
C -14 ILLEGAL MF VALUE.
C -15 ILLEGAL BREAKPOINTS - NOT STRICTLY INCREASING.
C -16 INSUFFICIENT STORAGE FOR WORK OR IWORK.
C
C-------------------------------------------------------------------------------
C
C SUMMARY OF ALL PACKAGE ROUTINES
C
C PDECOL - STORAGE ALLOCATION, ERROR CHECKING, INITIALIZATION, REPEATED
C CALLS TO STIFIB TO ADVANCE THE TIME.
C
C INTERP - INTERPOLATES COMPUTED BASIS FUNCTION COEFFICIENTS TO THE
C DESIRED OUTPUT TIMES, TOUT, FOR USE BY VALUES.
C
C INITAL - INITIALIZATION, GENERATION AND STORAGE OF PIECEWISE POLY-
C NOMIAL SPACE BASIS FUNCTION VALUES AND DERIVATIVES, DET-
C ERMINES THE BASIS FUNCTION COEFFICINTS OF THE PIECEWISE
C POLYNOMIALS WHICH INTERPOLATE THE USERS INITIAL CONDITIONS.
C
C COLPNT - GENERATION OF REQUIRED COLLOCATION POINTS.
C
C BSPLVD - B-SPLINE PACKAGE ROUTINES WHICH ALLOW FOR EVALUATION OF
C BSPLVN ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE.
C INTERV
C
C VALUES - GENERATION AT ANY POINT(S) OF VALUES OF THE COMPUTED
C APPROXIMATE SOLUTION AND ITS DERIVATIVES WHICH ARE
C PIECEWISE POLYNOMIALS. THE SUBROUTINE IS CALLED ONLY BY
C THE USER.
C
C STIFIB - CORE INTEGRATOR, TAKES SINGLE TIME STEPS TO ADVANCE THE
C TIME. ASSEMBLES AND SOLVES THE PROPER NONLINEAR EQUATIONS
C WHICH ARE RELATED TO USE OF ADAMS OR GEAR TYPE INTEGRATION
C FORMULAS. CHOOSES PROPER STEP SIZE AND INTEGRATION FORMULA
C ORDER TO MAINTAIN A DESIRED ACCURACY. DESIGNED FOR ODE
C PROBLEMS OF THE FORM A * (DY/DT) = G(T,Y).
C
C COSET - GENERATES INTEGRATION FORMULA AND ERROR CONTROL COEFFICIENTS.
C
C RES - COMPUTES RESIDUAL VECTORS USED IN SOLVING THE NONLINEAR
C EQUATIONS BY A MODIFIED NEWTON METHOD.
C
C DIFFUN - COMPUTES A**-1 * G(T,Y) WHERE A AND G ARE AS ABOVE (STIFIB).
C
C ADDA - ADDS THE A MATRIX TO A GIVEN MATRIX IN BAND FORM.
C
C EVAL - EVALUATES THE COMPUTED PIECEWISE POLYNOMIAL SOLUTION AND
C DERIVATIVES AT COLLOCATION POINTS.
C
C GFUN - EVALUATES THE FUNCTION G(T,Y) BY CALLING EVAL AND THE USER
C SUBROUTINES F AND BNDRY.
C
C PSETIB - GENERATES PROPER JACOBIAN MATRICES REQUIRED BY THE MODIFIED
C NEWTON METHOD.
C
C DIFFF - PERFORMS SAME ROLE AS THE USER ROUTINE DERIVF. COMPUTES
C DERIVATIVE APPROXIMATIONS BY USE OF FINITE DIFFERENCES.
C
C DECB - PERFORM AN LU DECOMPOSTION AND FORWARD AND BACKWARD
C SOLB SUBSTITUTION FOR SOLVING BANDED SYSTEMS OF LINEAR EQUATIONS.
C
C-----------------------------------------------------------------------
C
C
C STORAGE ALLOCATION
C
C SINCE PDECOL IS A DYNAMICALLY DIMENSIONED PROGRAM, MOST OF ITS
C WORKING STORAGE IS PROVIDED BY THE USER IN THE ARRAYS WORK AND IWORK.
C THE FOLLOWING GIVES A LIST OF THE ARRAYS WHICH MAKE UP THE CONTENTS
C WORK AND IWORK, THEIR LENGTHS, AND THEIR USES. WHEN MORE THAN ONE
C NAME IS GIVEN, IT INDICATES THAT DIFFERENT NAMES ARE USED FOR THE
C SAME ARRAY IN DIFFERENT PARTS OF THE PROGRAM. THE DIFFERENT NAMES
C OCCUR BECAUSE PDECOL IS AN AMALGAMATION OF SEVERAL OTHER CODES
C WRITTEN BY DIFFERENT PEOPLE AND WE HAVE TRIED TO LEAVE THE SEPARATE
C PARTS AS UNCHANGED FROM THEIR ORIGINAL VERSIONS AS POSSIBLE.
C
C
C NAMES LENGTH USE
C --------- ------------ -------------------------------------
C
C BC 4*NPDE**2 BOUNDARY CONDITION INFORMATION.
C WORK
C
C A 3*KORD*NCPTS BASIS FUNCTION VALUES AT COLLOCATION POINT
C WORK(IW1)
C
C XT NCPTS + KORD BREAKPOINT SEQUENCE FOR GENERATION OF BASI
C WORK(IW2) FUNCTION VALUES.
C
C XC NCPTS COLLOCATION POINTS.
C WORK(IW3)
C
C CMAX NPDE*NCPTS VALUES USED IN ESTIMATING TIME
C YMAX INTEGRATION ERRORS.
C WORK(IW4)
C
C ERROR NPDE*NCPTS TIME INTEGRATION ERRORS.
C WORK(IW5)
C
C SAVE1 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION
C WORK(IW6) METHOD.
C
C SAVE2 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION
C WORK(IW7) METHOD.
C
C SAVE3 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION
C WORK(IW8) METHOD.
C
C UVAL 3*NPDE WORKING STORAGE FOR VALUES OF U, UX, AND
C WORK(IW9) UXX AT ONE POINT.
C
C C NPDE*NCPTS* CURRENT BASIS FUNCTION COEFFICIENT VALUES
C Y (MAXDER+1) AND THEIR SCALED TIME DERIVATIVES.
C WORK(IW10)
C
C DFDU NPDE**2 WORKING STORAGE USED TO COMPUTE THE
C WORK(IW11) JACOBIAN MATRIX.
C
C DFDUX NPDE**2 WORKING STORAGE USED TO COMPUTE THE
C WORK(IW12) JACOBIAN MATRIX.
C
C DFDUXX NPDE**2 WORKING STORAGE USED TO COMPUTE THE
C WORK(IW13) JACOBIAN MATRIX.
C
C DBDU NPDE**2 BOUNDARY CONDITION INFORMATION.
C WORK(IW14)
C
C DBDUX NPDE**2 BOUNDARY CONDITION INFORMATION.
C WORK(IW15)
C
C DZDT NPDE BOUNDARY CONDITION INFORMATION.
C WORK(IW16)
C
C PW NPDE*NCPTS* STORAGE AND PROCESSING OF THE JACOBIAN
C WORK(IW17) (3*ML+1) MATRIX.
C
C ILEFT NCPTS POINTERS TO BREAKPOINT SEQUENCE FOR
C IWORK GENERATION OF BASIS FUNCTION VALUES.
C
C IPIV NPDE*NCPTS PIVOT INFORMATION FOR THE LU DECOMPOSED
C IWORK(IW18) JACOBIAN MATRIX PW.
C
C WHERE...
C
C NCPTS = KORD*NINT - NCC*(NINT-1)
C ML = NPDE*(KORD+IQUAD-1) - 1
C IQUAD = 1 IF KORD = 3 AND A NULL BOUNDARY CONDITION EXISTS
C IQUAD = 0 OTHERWISE
C MAXDER = 5 UNLESS OTHERWISE SET BY THE USER INTO /OPTION/.
C
C THE COMMON BLOCK /OPTION/ CONTAINS THE TWO VARIABLES NOGAUS AND
C MAXDER. NOGAUS IS SET .EQ. 0 IN THE BLOCK DATA. IT CAN BE CHANGED
C TO BE SET .EQ. 1 IF THE GAUSS-LEGENDRE COLLOCATION POINTS ARE NOT
C DESIRED WHEN NCC = 2 (SEE ABOVE AND COLPNT). MAXDER IS SET
C .EQ. 5 IN THE BLOCK DATA AND ITS VALUE REPRESENTS THE
C MAXIMUM ORDER OF TIME INTEGRATION FORMULA ALLOWED. ITS VALUE
C AFFECTS THE STORAGE REQUIRED IN WORK AND MAY BE CHANGED IF
C DESIRED. SEE COSET FOR RESTRICTIONS. THESE CHANGES MAY BE MADE BY
C THE USER BY ACCESSING /OPTION/ IN HIS CALLING PROGRAM (BEFORE THE
C FIRST CALL TO PDECOL) OR BY CHANGING THE DATA STATEMENT IN
C THE BLOCK DATA.
C
C-----------------------------------------------------------------------
C
C
C COMMUNICATION
C
C EACH SUBROUTINE IN THE PACKAGE CONTAINS A COMMUNICATION SUMMARY
C AS INDICATED BELOW.
C
C PACKAGE ROUTINES CALLED.. EVAL,INITAL,INTERP,STIFIB
C USER ROUTINES CALLED.. BNDRY
C CALLED BY.. USERS MAIN PROGRAM
C FORTRAN FUNCTIONS USED.. ABS,DMAX1,FLOAT,DSQRT
C-----------------------------------------------------------------------
SAVE
COMMON /GEAR0/ DTUSED,NQUSED,NSTEP,NFE,NJE
COMMON /GEAR1/ T,DTC,DTMN,DTMX,EPSC,UROUND,N,MFC,KFLAG,JSTART
COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W
COMMON /OPTION/ NOGAUS,MAXDER
COMMON /SIZES/ NIN,KOR,NC,NPD,NCPTS,NEQN,IQUAD
COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,
* IW11,IW12,IW13,IW14,IW15,IW16,IW17,IW18
COMMON /IOUNIT/ LOUT
DIMENSION WORK(KORD+NPDE*(4+9*NPDE)+(KORD+(NINT-1)*(KORD-NCC))*
* (3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4))),
* IWORK((NPDE+1)*(KORD+(NINT-1)*(KORD-NCC))), XBKPT(NINT+1)
IF (INDEX .EQ. 0) GO TO 60
IF (INDEX .EQ. 2) GO TO 70
IF (INDEX .EQ. 4) GO TO 80
IF (INDEX .EQ. 3) GO TO 90
C-----------------------------------------------------------------------
C SEVERAL CHECKS ARE MADE HERE TO DETERMINE IF THE INPUT PARAMETERS
C HAVE LEGAL VALUES. ERROR CHECKS ARE MADE ON INDEX, EPS, (T0-TOUT)*DT,
C DT, NINT, KORD, NCC, NPDE, MF, WHETHER THE BREAKPOINT SEQUENCE IS
C STRICTLY INCREASING, AND WHETHER THERE IS SUFFICIENT STORAGE
C PROVIDED FOR WORK AND IWORK. PROBLEM DEPENDENT PARAMETERS ARE
C CALCULATED AND PLACED IN COMMON.
C-----------------------------------------------------------------------
IERID = -6
IF (INDEX .NE. 1) GO TO 320
IERID = IERID - 1
IF (EPS .LE. 0.) GO TO 320
IERID = IERID - 1
IF ((T0-TOUT)*DT .GT. 0.) GO TO 320
IERID = IERID - 1
IF (DT .EQ. 0.0) GO TO 320
IERID = IERID - 1
NIN = NINT
IF (NIN .LT. 1) GO TO 320
IERID = IERID - 1
KOR = KORD
IF (KOR .LT. 3 .OR. KOR .GT. 20) GO TO 320
IERID = IERID - 1
NC = NCC
IF (NCC .LT. 2 .OR. NCC .GE. KOR) GO TO 320
IERID = IERID - 1
NPD = NPDE
NPDE2 = NPD*NPD
IF (NPDE .LT. 1) GO TO 320
IERID = IERID - 1
IF (MF.NE.22.AND.MF.NE.21.AND.MF.NE.12.AND.MF.NE.11) GO TO 320
IERID = IERID - 1
DO 10 K=1,NIN
IF(XBKPT(K) .GE. XBKPT(K+1)) GO TO 320
10 CONTINUE
NCPTS = KOR + (NIN - 1) * (KOR - NCC)
NEQN = NPDE * NCPTS
ML = (KOR-1)*NPDE - 1
MU = ML
MW = ML + ML + 1
N0W = NEQN*MW
IWSAVE = IWORK(1)
IISAVE = IWORK(2)
IW1 = 4*NPDE2 + 1
IW2 = IW1 + 3*KORD*NCPTS
IW3 = IW2 + NCPTS + KORD
IW4 = IW3 + NCPTS
IW5 = IW4 + NEQN
IW6 = IW5 + NEQN
IW7 = IW6 + NEQN
IW8 = IW7 + NEQN
IW9 = IW8 + NEQN
IW10 = IW9 + 3*NPDE
IW11 = IW10 + NEQN*(MAXDER+1)
IW12 = IW11 + NPDE2
IW13 = IW12 + NPDE2
IW14 = IW13 + NPDE2
IW15 = IW14 + NPDE2
IW16 = IW15 + NPDE2
IW17 = IW16 + NPDE
IW18 = NCPTS + 1
IERID = IERID - 1
IWSTOR = IW17 + NEQN*(3*ML+1) - 1
IISTOR = IW18 + NEQN - 1
IF ( IWSAVE .LT. IWSTOR .OR. IISAVE .LT. IISTOR ) GO TO 335
C-----------------------------------------------------------------------
C PERFORM INITIALIZATION TASKS. IF KORD .EQ. 3 THEN CALCULATE THE BAND-
C WIDTH OF THE ASSOCIATED MATRIX PROBLEM BY DETERMINING THE TYPE OF
C BOUNDARY CONDITIONS, THEN CHECK FOR SUFFICIENT STORAGE AGAIN.
C-----------------------------------------------------------------------
CALL INITAL(KOR,WORK(IW1),WORK(IW6),XBKPT,WORK(IW2),WORK(IW3),
* WORK(IW17),IWORK(IW18),IWORK)
IF(IQUAD .NE. 0) GO TO 280
IF( KOR .NE. 3 ) GO TO 40
CALL EVAL(1,NPDE,WORK(IW6),WORK(IW9),WORK(IW1),IWORK)
CALL BNDRY(T0,WORK(IW3),WORK(IW9),WORK(IW9+NPDE),WORK(IW14),
* WORK(IW15),WORK(IW16),NPDE)
DO 20 K=1,NPDE
I = K + NPDE*(K-1) - 1
IF(WORK(IW14+I) .EQ. 0.0 .AND. WORK(IW15+I) .EQ. 0.0)
* IQUAD = 1
20 CONTINUE
CALL EVAL(NCPTS,NPDE,WORK(IW6),WORK(IW9),WORK(IW1),IWORK)
CALL BNDRY(T0,WORK(IW3+NCPTS-1),WORK(IW9),WORK(IW9+NPDE),
* WORK(IW14),WORK(IW15),WORK(IW16),NPDE)
DO 30 K=1,NPDE
I = K + NPDE*(K-1) - 1
IF(WORK(IW14+I) .EQ. 0.0 .AND. WORK(IW15+I) .EQ. 0.0)
* IQUAD = 1
30 CONTINUE
ML = ML + IQUAD*NPDE
MU = ML
MW = ML + ML + 1
N0W = NEQN*MW
40 CONTINUE
IWSTOR = IW17 + NEQN*(3*ML+1) - 1
IF ( IWSAVE .LT. IWSTOR ) GO TO 335
C-----------------------------------------------------------------------
C IF INITIAL VALUES OF CMAX OTHER THAN THOSE SET BELOW ARE DESIRED,
C THEY SHOULD BE SET HERE. ALL CMAX(I) MUST BE POSITIVE.
C HAVING PROPER VALUES OF CMAX FOR THE PROBLEM BEING SOLVED IS AS
C IMPORTANT AS CHOOSING EPS (SEE ABOVE), SINCE ERRORS ARE
C MEASURED RELATIVE TO CMAX. IF VALUES FOR DTMN OR DTMX, THE
C BOUNDS ON DABS(DT), OTHER THAN THOSE BELOW ARE DESIRED, THEY
C SHOULD BE SET BELOW.
C-----------------------------------------------------------------------
DO 50 I = 1,NEQN
I1 = I - 1
WORK(IW4+I1) = DABS(WORK(IW6+I1))
IF (WORK(IW4+I1) .LT. 1.) WORK(IW4+I1) = 1.
50 WORK(IW10+I1) = WORK(IW6+I1)
N = NEQN
T = T0
DTC = DT
DTMN = DABS(DT)
DTUSED = 0.
EPSC = EPS
MFC = MF
JSTART = 0
EPSJ = DSQRT(UROUND)
NM1 = NEQN - 1
N0ML = NEQN*ML
NHCUT = 0
KFLAG = 0
TOUTP = T0
IF ( T0 .EQ. TOUT ) GO TO 360
60 DTMX = DABS(TOUT-TOUTP)*10.
GO TO 140
C
70 DTMX = DABS(TOUT-TOUTP)*10.
IF ((T-TOUT)*DTC .GE. 0.) GO TO 340
GO TO 150
C
80 IF ((T-TOUT)*DTC .GE. 0.) GO TO 300
JSTART = -1
EPSC = EPS
MFC = MF
GO TO 100
C
90 DTMX = DT
100 IF ((T+DTC) .EQ. T) WRITE(LOUT,110)
110 FORMAT(36H WARNING.. T + DT = T ON NEXT STEP.)
C-----------------------------------------------------------------------
C TAKE A TIME STEP BY CALLING THE INTEGRATOR.
C-----------------------------------------------------------------------
CALL STIFIB (NEQN,WORK(IW10),WORK(IW4),WORK(IW5),WORK(IW6),
* WORK(IW7),WORK(IW8),WORK(IW17),IWORK(IW18),WORK,IWORK)
C
KGO = 1 - KFLAG
GO TO (120, 160, 220, 260, 280), KGO
C KFLAG = 0, -1, -2, -3 -4
C
120 CONTINUE
C-----------------------------------------------------------------------
C NORMAL RETURN FROM INTEGRATOR.
C
C THE WEIGHTS CMAX(I) ARE UPDATED. IF DIFFERENT VALUES ARE DESIRED,
C THEY SHOULD BE SET HERE. A TEST IS MADE FOR EPS BEING TOO SMALL
C FOR THE MACHINE PRECISION.
C
C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY
C STEP SHOULD BE INSERTED HERE.
C
C IF INDEX = 3, SAVE1 IS SET TO THE CURRENT C VALUES ON RETURN.
C IF INDEX = 2, DT IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF
C ERROR), AND THEN THE CURRENT C VALUES ARE PUT IN SAVE1 ON RETURN.
C FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR
C UNLESS TOUT HAS BEEN REACHED. THEN INTERPOLATED VALUES OF C ARE
C COMPUTED AND STORED IN SAVE1 ON RETURN.
C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE
C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 340 INSTEAD OF 360.
C-----------------------------------------------------------------------
D = 0.
DO 130 I = 1,NEQN
I1 = I - 1
AYI = DABS(WORK(IW10+I1))
WORK(IW4+I1) = DMAX1(WORK(IW4+I1), AYI)
130 D = D + (AYI/WORK(IW4+I1))**2
D = D*(UROUND/EPS)**2
IF (D .GT. FLOAT(NEQN)) GO TO 240
IF (INDEX .EQ. 3) GO TO 340
IF (INDEX .EQ. 2) GO TO 150
140 IF ((T-TOUT)*DTC .LT. 0.) GO TO 100
CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6))
GO TO 360
C
150 IF (((T+DTC)-TOUT)*DTC .LE. 0.) GO TO 100
IF (DABS(T-TOUT) .LE. 100.*UROUND*DTMX) GO TO 340
IF ((T-TOUT)*DTC .GE. 0.) GO TO 340
DTC = (TOUT - T)*(1. - 4.*UROUND)
JSTART = -1
GO TO 100
C-----------------------------------------------------------------------
C ON AN ERROR RETURN FROM INTEGRATOR, AN IMMEDIATE RETURN OCCURS IF
C KFLAG = -2 OR -4, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE.
C TO RECOVER, DT AND DTMN ARE REDUCED BY A FACTOR OF .1 UP TO 10
C TIMES BEFORE GIVING UP.
C-----------------------------------------------------------------------
160 WRITE (LOUT,170) T
170 FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/
* 41H ERROR TEST FAILED WITH DABS(DT) = DTMIN/)
180 IF (NHCUT .EQ. 10) GO TO 200
NHCUT = NHCUT + 1
DTMN = .1*DTMN
DTC = .1*DTC
WRITE (LOUT,190) DTC
190 FORMAT(25H DT HAS BEEN REDUCED TO ,E16.8,
* 26H AND STEP WILL BE RETRIED//)
JSTART = -1
GO TO 100
C
200 WRITE (LOUT,210)
210 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//)
GO TO 340
C
220 WRITE (LOUT,230) T,DTC
230 FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,6H DT =,
* E16.8/52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//)
GO TO 340
C
240 WRITE (LOUT,250) T
250 FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/
* 56H EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/)
KFLAG = -2
GO TO 340
C
260 WRITE (LOUT,270) T
270 FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/
* 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/)
GO TO 180
C
280 WRITE (LOUT,290)
290 FORMAT(//28H SINGULAR MATRIX ENCOUNTERED,
* 35H PROBABLY DUE TO STORAGE OVERWRITES//)
KFLAG = -4
GO TO 340
C
300 WRITE(LOUT,310) T,TOUT,DTC
310 FORMAT(//45H INDEX = -1 ON INPUT WITH (T-TOUT)*DT .GE. 0./
* 4H T =,E16.8,9H TOUT =,E16.8,8H DTC =,E16.8/
* 44H INTERPOLATION WAS DONE AS ON NORMAL RETURN./
* 41H DESIRED PARAMETER CHANGES WERE NOT MADE.)
CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6))
INDEX = -5
RETURN
C
320 WRITE(LOUT,330) IERID
330 FORMAT(//24H ILLEGAL INPUT...INDEX= ,I3//)
INDEX = IERID
RETURN
C
335 WRITE(LOUT,336) IWSTOR,IWSAVE,IISTOR,IISAVE
336 FORMAT(//21H INSUFFICIENT STORAGE/24H WORK MUST BE OF LENGTH,
* I10,5X,12HYOU PROVIDED,I10/24H IWORK MUST BE OF LENGTH,I10,5X,
* 12HYOU PROVIDED,I10//)
INDEX = IERID
RETURN
C
340 TOUT = T
DO 350 I = 1,NEQN
I1 = I - 1
350 WORK(IW6+I1) = WORK(IW10+I1)
360 INDEX = KFLAG
TOUTP = TOUT
DT = DTUSED
IF (KFLAG .NE. 0) DT = DTC
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE VALUES(X,USOL,SCTCH,NDIM1,NDIM2,NPTS,NDERV,WORK)
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C SUBROUTINE VALUES COMPUTES THE SOLUTION U AND THE FIRST NDERV
C DERIVATIVES OF U AT THE NPTS POINTS X AND AT TIME TOUT AND RETURNS
C THEM IN THE ARRAY USOL. THIS ROUTINE MUST BE USED TO OBTAIN
C SOLUTION VALUES SINCE PDECOL DOES NOT RETURN ANY SOLUTION VALUES
C TO THE USER. SEE PDECOL.
C
C THE CALLING PARAMETERS ARE...
C X = AN ARBITRARY VECTOR OF SPATIAL POINTS OF LENGTH NPTS AT
C WHICH THE SOLUTION AND THE FIRST NDERV DERIVATIVE VALUES
C ARE TO BE CALCULATED. IF X .LT. XLEFT ( X .GT. XRIGHT )
C THEN THE PIECEWISE POLYNOMIAL OVER THE LEFTMOST ( RIGHT-
C MOST ) INTERVAL IS EVALUATED TO CALCULATE THE SOLUTION
C VALUES AT THIS UNUSUAL VALUE OF X. SEE PDECOL.
C
C USOL = AN ARRAY WHICH CONTAINS THE SOLUTION AND THE FIRST
C NDERV DERIVATIVES OF THE SOLUTION AT ALL THE POINTS IN
C THE INPUT VECTOR X. IN PARTICULAR, USOL(J,I,K) CONTAINS
C THE VALUE OF THE (K-1)-ST DERIVATIVE OF THE J-TH PDE
C COMPONENT AT THE I-TH POINT OF THE X VECTOR FOR
C J = 1 TO NPDE, I = 1 TO NPTS, AND K = 1 TO NDERV+1.
C
C SCTCH = A USER SUPPLIED WORKING STORAGE ARRAY OF LENGTH AT LEAST
C KORD*(NDERV+1). SEE BELOW AND PDECOL FOR DEFINITIONS OF
C THESE PARAMETERS.
C
C NDIM1 = THE FIRST DIMENSION OF THE OUTPUT ARRAY USOL IN THE CALLING
C PROGRAM. NDIM1 MUST BE .GE. NPDE.
C
C NDIM2 = THE SECOND DIMENSION OF THE OUTPUT ARRAY USOL IN THE
C CALLING PROGRAM. NDIM2 MUST BE .GE. NPTS.
C
C NPTS = THE NUMBER OF POINTS IN THE X VECTOR.
C
C NDERV = THE NUMBER OF DERIVATIVE VALUES OF THE SOLUTION THAT ARE
C TO BE CALCULATED. NDERV SHOULD BE LESS THAN KORD SINCE
C THE KORD-TH DERIVATIVE OF A POLYNOMIAL OF DEGREE KORD-1
C IS EQUAL TO ZERO. SEE PDECOL.
C
C WORK = THE USERS WORKING STORAGE ARRAY WHICH IS USED IN THIS CASE
C TO PROVIDE THE CURRENT BASIS FUNCTION COEFFICIENTS AND THE
C PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE.
C
C PACKAGE ROUTINES CALLED.. BSPLVD,INTERV
C USER ROUTINES CALLED.. NONE
C CALLED BY.. USERS MAIN PROGRAM
C FORTRAN FUNCTIONS USED.. NONE
C
C-----------------------------------------------------------------------
SAVE ILEFT, MFLAG
COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IDUM(12)
COMMON /OPTION/ NOGAUS,MAXDER
DIMENSION USOL(NDIM1,NDIM2,NDERV+1),X(NPTS),SCTCH(KORD*(NDERV+1)),
* WORK(KORD+NPDE*(4+9*NPDE)+(KORD+(NINT-1)*(KORD-NCC))*
* (3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4)))
DATA ILEFT/0/, MFLAG/0/
NDERV1 = NDERV + 1
DO 20 IPTS=1,NPTS
CALL INTERV(WORK(IW2),NCPTS,X(IPTS),ILEFT,MFLAG)
CALL BSPLVD(WORK(IW2),KORD,X(IPTS),ILEFT,SCTCH,NDERV1)
IK = ILEFT - KORD
DO 10 M=1,NDERV1
I1 = (M-1)*KORD
DO 10 K=1,NPDE
USOL(K,IPTS,M) = 0.
DO 10 I=1,KORD
I2 = (I+IK-1)*NPDE + IW6 - 1
USOL(K,IPTS,M) = USOL(K,IPTS,M) + WORK(I2+K) * SCTCH(I+I1)
10 CONTINUE
20 CONTINUE
RETURN
END
C
C
C ##############################################################################
C
C
BLOCK DATA
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C IN THE FOLLOWING DATA STATEMENT, SET..
C LOUT = THE LOGICAL UNIT NUMBER FOR THE OUTPUT OF MESSAGES DURING
C THE INTEGRATION.
C NOGAUS = SET .EQ. 1 IF THE GAUSS-LEGENDRE COLLOCATION POINTS ARE
C NOT DESIRED WHEN NCC = 2 (SEE PDECOL AND COLPNT).
C MAXDER = SET .EQ. 5. ITS VALUE REPRESENTS THE MAXIMUM ORDER OF
C THE TIME INTEGRATION ALLOWED. ITS VALUE AFFECTS THE STOR-
C AGE REQUIRED IN WORK AND MAY BE CHANGED IF DESIRED
C (SEE COSET FOR RESTRICTIONS).
C UROUND = THE UNIT ROUNDOFF OF THE MACHINE, I.E. THE SMALLEST
C POSITIVE U SUCH THAT 1. + U .NE. 1. ON THE MACHINE.
C-------------------------------------------------------------------------------
COMMON /GEAR1/ DUM(5),UROUND,IDUM(4)
COMMON /OPTION/ NOGAUS,MAXDER
COMMON /IOUNIT/ LOUT
C***
C*** UROUND SET TO SINGLE PRECISION FOR A SUN SPARC2
C***
C*** DATA LOUT,NOGAUS,MAXDER,UROUND/6, 0, 5, 5.960464D-08/
C
DATA LOUT,NOGAUS,MAXDER,UROUND/66, 0, 5, 2.22D-16/
END
C
C
C ##############################################################################
C
C
SUBROUTINE INITAL(K,A,RHS,X,XT,XC,PW,IPIV,ILEFT)
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C INITAL IS CALLED ONLY ONCE BY PDECOL TO PERFORM INITIALIZATION TASKS.
C THESE TASKS INCLUDE - 1) DEFINING THE PIECEWISE POLYNOMIAL SPACE
C BREAKPOINT SEQUENCE, 2) CALLING THE SUBROUTINE COLPNT TO DEFINE THE
C REQUIRED COLLOCATION POINTS, 3) DEFING THE PIECEWISE POLYNOMIAL SPACE
C BASIS FUNCTION VALUES (PLUS FIRST AND SECOND DERIVATIVE VALUES) AT
C THE COLLOCATION POINTS, AND 4) DEFINING THE INITIAL BASIS FUNCTION
C COEFFICIENTS WHICH DETERMINE THE PIECEWISE POLYNOMIAL WHICH
C INTERPOLATES THE USER SUPPLIED (UINIT) INITIAL CONDITION FUNCTION(S)
C AT THE COLLOCATION POINTS.
C
C K = ORDER OF PIECEWISE POLYNOMIAL SPACE.
C A = BASIS FUNCTION VALUES GENERATED BY INITAL.
C RHS = TEMPORARY STORAGE USED TO RETURN INITIAL CONDITION COEFFICIENT
C VALUES.
C X = USER DEFINED PIECEWISE POLYNOMIAL BREAKPOINTS.
C XT = PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE GENERATED BY INITAL.
C XC = COLLOCATION POINTS GENERATED BY INITAL.
C PW = STORAGE FOR BAND MATRIX USED TO GENERATE INITIAL
C COEFFICIENT VALUES.
C IPIV = PIVOT INFORMATION FOR LINEAR EQUATION SOLVER DECB-SOLB.
C ILEFT = POINTERS TO BREAKPOINT SEQUENCE GENERATED BY INITAL.
C
C PACKAGE ROUTINES CALLED.. BSPLVD,COLPNT,DECB,INTERV,SOLB
C USER ROUTINES CALLED.. UINIT
C CALLED BY.. PDECOL
C FORTRAN FUNCTIONS USED.. MAX0,MIN0
C-------------------------------------------------------------------------------
COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IER
COMMON /GEAR9/ EPSJ,R0,ML,MU,IDUM(3),N0W
DIMENSION A(K,3,NCPTS),RHS(NEQN),X(NINT+1),XT(NCPTS+KORD),
* XC(NCPTS),PW(NEQN*(3*ML+1)),
* IPIV(NEQN),ILEFT(NCPTS)
MFLAG = -2
IER = 0
C-----------------------------------------------------------------------
C SET UP THE PIECEWISE POLYNOMIAL SPACE BREAKPOINT SEQUENCE.
C-----------------------------------------------------------------------
KRPT = KORD - NCC
DO 10 I=1,KORD
XT(NCPTS+I) = X(NINT+1)
10 XT(I) = X(1)
DO 20 I=2,NINT
I1 = (I-2)*KRPT + KORD
DO 20 J=1,KRPT
20 XT(I1+J) = X(I)
C-----------------------------------------------------------------------
C SET UP COLLOCATION POINTS ARRAY XC.
C-----------------------------------------------------------------------
CALL COLPNT(X, XC, XT)
C-----------------------------------------------------------------------
C GENERATE THE ILEFT ARRAY. STORE THE BASIS FUNCTION VALUES IN THE
C ARRAY A. THE ARRAY A IS DIMENSIONED A(KORD,3,NCPTS) AND A(K,J,I)
C CONTAINS THE VALUE OF THE (J-1)-ST DERIVATIVE (J = 1,2,3) OF THE K-TH
C NONZERO BASIS FUNCTION (K = 1, ... ,KORD) AT THE I-TH COLLOCATION
C POINT (I = 1, ... ,NCPTS). SET UP RHS FOR INTERPOLATING THE INITIAL
C CONDITIONS AT THE COLLOCATION POINTS. SET THE INTERPOLATION MATRIX
C INTO THE BANDED MATRIX PW.
C-----------------------------------------------------------------------
DO 30 I=1,N0W
30 PW(I) = 0.
DO 40 I=1,NCPTS
CALL INTERV(XT,NCPTS,XC(I),ILEFT(I),MFLAG)
CALL BSPLVD(XT,KORD,XC(I),ILEFT(I),A(1,1,I),3)
I1 = NPDE * (I-1)
CALL UINIT(XC(I),RHS(I1+1),NPDE)
ICOL = ILEFT(I) - I - 1
JL = MAX0(1,I+2-NCPTS)
JU = MIN0(KORD,KORD+I-2)
DO 40 J=JL,JU
J1 = I1 + NEQN * (NPDE * (ICOL + J) - 1)
DO 40 JJ=1,NPDE
40 PW(JJ+J1) = A(J,1,I)
C-----------------------------------------------------------------------
C LU DECOMPOSE THE MATRIX PW.
C-----------------------------------------------------------------------
CALL DECB (NEQN,NEQN,ML,MU,PW,IPIV,IER)
IF ( IER .NE. 0 ) RETURN
C-----------------------------------------------------------------------
C SOLVE THE LINEAR SYSTEM PW*Z = RHS. THIS GIVES THE BASIS FUNCTION
C COEFFICIENTS FOR THE INITIAL CONDITIONS.
C-----------------------------------------------------------------------
CALL SOLB (NEQN,NEQN,ML,MU,PW,RHS,IPIV)
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE COLPNT(X, XC, XT)
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C COLPNT IS CALLED ONLY ONCE BY INITAL TO DEFINE THE REQUIRED COLLOCA-
C TION POINTS WHICH ARE TO BE USED WITH THE USER SELECTED PIECEWISE
C POLYNOMIAL SPACE. THE COLLOCATION POINTS ARE CHOSEN SUCH THAT THEY
C ARE EITHER THE POINTS AT WHICH THE PIECEWISE POLYNOMIAL SPACE BASIS
C FUNCTIONS ATTAIN THEIR UNIQUE MAXIMUM VALUES, OR, THE GAUSS-LEGENDRE
C QUADRATURE POINTS WITHIN EACH PIECEWISE POLYNOMIAL SPACE SUBINTERVAL,
C DEPENDING UPON THE SPACE BEING USED AND THE DESIRE OF THE USER.
C
C X = USER DEFINED PIECEWISE POLYNOMIAL BREAKPOINTS.
C XC = COLLOCATION POINTS DEFINED BY COLPNT.
C XT = PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE.
C
C PACKAGE ROUTINES CALLED.. BSPLVD,INTERV
C USER ROUTINES CALLED.. NONE
C CALLED BY.. INITAL
C FORTRAN FUNCTIONS USED.. NONE
C-------------------------------------------------------------------------------
SAVE ILEFT
COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
COMMON /OPTION/ NOGAUS,MAXDER
DIMENSION RHO(40),X(NINT+1),XC(NCPTS),XT(NCPTS+KORD)
DATA ILEFT/0/
C-----------------------------------------------------------------------
C IF THE VARIABLE NOGAUS IN THE COMMON BLOCK /OPTION/ IS SET .EQ. 1,
C THE USE OF THE GAUSS-LEGENDRE POINTS IS PROHIBITED FOR ALL CASES.
C NOGAUS IS CURRENTLY SET .EQ. 0 BY A DATA STATEMENT IN THE BLOCK DATA.
C THE USER MAY CHANGE THIS AS DESIRED.
C-----------------------------------------------------------------------
IF ( NCC .NE. 2 .OR. NOGAUS .EQ. 1 ) GO TO 200
C-----------------------------------------------------------------------
C COMPUTE THE COLLOCATION POINTS TO BE AT THE GAUSS-LEGENDRE POINTS IN
C EACH PIECEWISE POLYNOMIAL SPACE SUBINTERVAL. THE ARRAY RHO IS SET TO
C CONTAIN THE GAUSS-LEGENDRE POINTS FOR THE STANDARD INTERVAL (-1,1).
C-----------------------------------------------------------------------
IPTS = KORD - 2
GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,
* 180),IPTS
10 RHO(1) = 0.
GO TO 190
20 RHO(2) = .577350269189626D-00
RHO(1) = - RHO(2)
GO TO 190
30 RHO(3) = .774596669241483D-00
RHO(1) = - RHO(3)
RHO(2) = 0.
GO TO 190
40 RHO(3) = .339981043584856D-00
RHO(2) = - RHO(3)
RHO(4) = .861136311594053D-00
RHO(1) = - RHO(4)
GO TO 190
50 RHO(4) = .538469310105683D-00
RHO(2) = - RHO(4)
RHO(5) = .906179845938664D-00
RHO(1) = - RHO(5)
RHO(3) = 0.
GO TO 190
60 RHO(4) = .238619186083197D-00
RHO(3) = - RHO(4)
RHO(5) = .661209386466265D-00
RHO(2) = - RHO(5)
RHO(6) = .932469514203152D-00
RHO(1) = - RHO(6)
GO TO 190
70 RHO(5) = .405845151377397D-00
RHO(3) = - RHO(5)
RHO(6) = .741531185599394D-00
RHO(2) = - RHO(6)
RHO(7) = .949107912342759D-00
RHO(1) = - RHO(7)
RHO(4) = 0.
GO TO 190
80 RHO(5) = .183434642495650D-00
RHO(4) = - RHO(5)
RHO(6) = .525532409916329D-00
RHO(3) = - RHO(6)
RHO(7) = .796666477413627D-00
RHO(2) = - RHO(7)
RHO(8) = .960289856497536D-00
RHO(1) = - RHO(8)
GO TO 190
90 RHO( 5) = .0
RHO( 6) = .324253423403809D-00
RHO( 7) = .613371432700590D-00
RHO( 8) = .836031107326636D-00
RHO( 9) = .968160239507626D-00
DO 95 I=1,4
95 RHO(I) = -RHO(10-I)
GO TO 190
100 RHO( 6) = .148874338981631D-00
RHO( 7) = .433395394129247D-00
RHO( 8) = .679409568299024D-00
RHO( 9) = .865063366688984D-00
RHO(10) = .973906528517172D-00
DO 105 I=1,5
105 RHO(I) = -RHO(11-I)
GO TO 190
110 RHO( 6) = .0
RHO( 7) = .269543155952345D-00
RHO( 8) = .519096129206812D-00
RHO( 9) = .730152005574049D-00
RHO(10) = .887062599768095D-00
RHO(11) = .978228658146057D-00
DO 115 I=1,5
115 RHO(I) = -RHO(12-I)
GO TO 190
120 RHO( 7) = .125233408511469D-00
RHO( 8) = .367831498998180D-00
RHO( 9) = .587317954286617D-00
RHO(10) = .769902674194305D-00
RHO(11) = .904117256370475D-00
RHO(12) = .981560634246719D-00
DO 125 I=1,6
125 RHO(I) = -RHO(13-I)
GO TO 190
130 RHO( 7) = .0
RHO( 8) = .230458315955135D-00
RHO( 9) = .448492751036447D-00
RHO(10) = .642349339440340D-00
RHO(11) = .801578090733310D-00
RHO(12) = .917598399222978D-00
RHO(13) = .984183054718588D-00
DO 135 I=1,6
135 RHO(I) = -RHO(14-I)
GO TO 190
140 RHO( 8) = .108054948707344D-00
RHO( 9) = .319112368927890D-00
RHO(10) = .515248636358154D-00
RHO(11) = .687292904811685D-00
RHO(12) = .827201315069765D-00
RHO(13) = .928434883663574D-00
RHO(14) = .986283808696812D-00
DO 145 I=1,7
145 RHO(I) = -RHO(15-I)
GO TO 190
150 RHO( 8) = .0
RHO( 9) = .201194093997435D-00
RHO(10) = .394151347077563D-00
RHO(11) = .570972172608539D-00
RHO(12) = .724417731360170D-00
RHO(13) = .848206583410427D-00
RHO(14) = .937273392400706D-00
RHO(15) = .987992518020485D-00
DO 155 I = 1,7
155 RHO(I) = -RHO(16-I)
GO TO 190
160 RHO( 9) = .950125098376374D-01
RHO(10) = .281603550779259D-00
RHO(11) = .458016777657227D-00
RHO(12) = .617876244402644D-00
RHO(13) = .755404408355003D-00
RHO(14) = .865631202387832D-00
RHO(15) = .944575023073233D-00
RHO(16) = .989400934991650D-00
DO 165 I=1,8
165 RHO(I) = -RHO(17-I)
GO TO 190
170 RHO( 9) = .0
RHO(10) = .178484181495848D-00
RHO(11) = .351231763453876D-00
RHO(12) = .512690537086477D-00
RHO(13) = .657671159216691D-00
RHO(14) = .781514003896801D-00
RHO(15) = .880239153726986D-00
RHO(16) = .950675521768768D-00
RHO(17) = .990575475314417D-00
DO 175 I=1,8
175 RHO(I) = -RHO(18-I)
GO TO 190
180 RHO(10) = .847750130417353D-01
RHO(11) = .251886225691506D-00
RHO(12) = .411751161462843D-00
RHO(13) = .559770831073948D-00
RHO(14) = .691687043060353D-00
RHO(15) = .803704958972523D-00
RHO(16) = .892602466497556D-00
RHO(17) = .955823949571398D-00
RHO(18) = .991565168420931D-00
DO 185 I=1,9
185 RHO(I) = -RHO(19-I)
C-----------------------------------------------------------------------
C COMPUTE THE GAUSS-LEGENDRE COLLOCATION POINTS IN EACH SUBINTERVAL.
C-----------------------------------------------------------------------
190 DO 195 I=1,NINT
FAC = ( X(I+1) - X(I) ) * .5
DO 195 J = 1,IPTS
KNOT = IPTS * (I-1) + J + 1
195 XC(KNOT) = X(I) + FAC * ( RHO(J) + 1. )
XC(1) = X(1)
XC(NCPTS) = X(NINT+1)
RETURN
C-----------------------------------------------------------------------
C COMPUTE THE COLLOCATION POINTS TO BE AT THE POINTS WHERE THE BASIS
C FUNCTIONS ATTAIN THEIR MAXIMA. A BISECTION METHOD IS USED TO FIND
C THE POINTS TO MACHINE PRECISION. THIS PROCESS COULD BE SPEEDED UP
C BY USING A SECANT METHOD IF DESIRED.
C-----------------------------------------------------------------------
200 ITOP = NCPTS - 1
MFLAG = -2
XC(1) = X(1)
XC(NCPTS) = X(NINT+1)
DO 240 I=2,ITOP
XOLD = 1.E+20
XL = XT(I)
XR = XT(I+KORD)
210 XNEW = .5 * (XL + XR)
IF( XOLD .EQ. XNEW ) GO TO 240
CALL INTERV(XT,NCPTS,XNEW,ILEFT,MFLAG)
CALL BSPLVD(XT,KORD,XNEW,ILEFT,RHO,2)
DO 220 J=1,KORD
IF( I .EQ. J + ILEFT - KORD ) GO TO 230
220 CONTINUE
230 XVAL = RHO(KORD+J)
IF( XVAL .EQ. 0.0 ) XR = XNEW
IF( XVAL .GT. 0.0 ) XL = XNEW
IF( XVAL .LT. 0.0 ) XR = XNEW
XOLD = XNEW
GO TO 210
240 XC(I) = XR
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE BSPLVD ( XT, K, X, ILEFT, VNIKX, NDERIV )
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE
C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE.
C SEE REFERENCE BELOW.
C
C CALCULATES THE VALUE AND THE FIRST NDERIV-1 DERIVATIVES OF ALL
C B-SPLINES WHICH DO NOT VANISH AT X. THE ROUTINE FILLS THE TWO-
C DIMENSIONAL ARRAY VNIKX(J,IDERIV), J=IDERIV, ... ,K WITH NONZERO
C VALUES OF B-SPLINES OF ORDER K+1-IDERIV, IDERIV=NDERIV, ... ,1, BY
C REPEATED CALLS TO BSPLVN.
C
C XT = PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE.
C K = ORDER OF THE PIECEWISE POLYNOMIAL SPACE.
C X = POINT AT WHICH THE B-SPLINE IS TO BE EVALUATED.
C ILEFT = POINTER TO THE BREAKPOINT SEQUENCE.
C VNIKX = TABLE OF B-SPLINE VALUES AND DERIVATIVES.
C NDERIV = DETERMINES NUMBER OF DERIVATIVES TO BE GENERATED.
C
C REFERENCE
C
C DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J.
C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472.
C
C PACKAGE ROUTINES CALLED.. BSPLVN
C USER ROUTINES CALLED.. NONE
C CALLED BY.. COLPNT,INITAL,VALUES
C FORTRAN FUNCTIONS USED.. FLOAT,MAX0
C-------------------------------------------------------------------------------
COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
DIMENSION XT(NCPTS+KORD),VNIKX(K,NDERIV)
DIMENSION A(20,20)
KO = K + 1 - NDERIV
CALL BSPLVN(XT,KO,1,X,ILEFT,VNIKX(NDERIV,NDERIV))
IF (NDERIV .LE. 1) GO TO 120
IDERIV = NDERIV
DO 20 I=2,NDERIV
IDERVM = IDERIV-1
DO 10 J=IDERIV,K
10 VNIKX(J-1,IDERVM) = VNIKX(J,IDERIV)
IDERIV = IDERVM
CALL BSPLVN(XT,0,2,X,ILEFT,VNIKX(IDERIV,IDERIV))
20 CONTINUE
DO 40 I=1,K
DO 30 J=1,K
30 A(I,J) = 0.
40 A(I,I) = 1.
KMD = K
DO 110 M=2,NDERIV
KMD = KMD - 1
FKMD = FLOAT(KMD)
I = ILEFT
J = K
50 JM1 = J-1
IPKMD = I + KMD
DIFF = XT(IPKMD) -XT(I)
IF (JM1 .EQ. 0) GO TO 80
IF (DIFF .EQ. 0.) GO TO 70
DO 60 L=1,J
60 A(L,J) = (A(L,J) - A(L,J-1))/DIFF*FKMD
70 J = JM1
I = I - 1
GO TO 50
80 IF (DIFF .EQ. 0.) GO TO 90
A(1,1) = A(1,1)/DIFF*FKMD
90 DO 110 I=1,K
V = 0.
JLOW = MAX0(I,M)
DO 100 J=JLOW,K
100 V = A(I,J)*VNIKX(J,M) + V
110 VNIKX(I,M) = V
120 RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE BSPLVN ( XT, JHIGH, INDEX, X, ILEFT, VNIKX )
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE
C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE.
C SEE REFERENCE BELOW.
C
C CALCULATES THE VALUE OF ALL POSSIBLY NONZERO B-SPLINES AT THE
C POINT X OF ORDER MAX(JHIGH,(J+1)(INDEX-1)) FOR THE BREAKPOINT SEQ-
C UENCE XT. ASSUMING THAT XT(ILEFT) .LE. X .LE. XT(ILEFT+1), THE ROUT-
C INE RETURNS THE B-SPLINE VALUES IN THE ONE DIMENSIONAL ARRAY VNIKX.
C
C FOR DEFINITIONS OF CALLING ARGUMENTS SEE ABOVE AND BSPLVD.
C
C REFERENCE
C
C DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J.
C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472.
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. NONE
C CALLED BY.. BSPLVD
C FORTRAN FUNCTIONS USED.. NONE
C-------------------------------------------------------------------------------
SAVE J,DELTAM,DELTAP
COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
DIMENSION DELTAM(20),DELTAP(20)
DIMENSION XT(NCPTS+KORD),VNIKX(*)
DATA J/1/,DELTAM/20*0.D-00/,DELTAP/20*0.D-00/
GO TO (10,20),INDEX
10 J = 1
VNIKX(1) = 1.
IF (J .GE. JHIGH) GO TO 40
20 IPJ = ILEFT+J
DELTAP(J) = XT(IPJ) - X
IMJP1 = ILEFT-J+1
DELTAM(J) = X - XT(IMJP1)
VMPREV = 0.
JP1 = J+1
DO 30 L=1,J
JP1ML = JP1-L
VM = VNIKX(L)/(DELTAP(L) + DELTAM(JP1ML))
VNIKX(L) = VM*DELTAP(L) + VMPREV
30 VMPREV = VM*DELTAM(JP1ML)
VNIKX(JP1) = VMPREV
J = JP1
IF (J .LT. JHIGH) GO TO 20
40 RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE INTERV ( XT, LXT, X, ILEFT, MFLAG )
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE
C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE.
C SEE REFERENCE BELOW.
C
C COMPUTES LARGEST ILEFT IN (1,LXT) SUCH THAT XT(ILEFT) .LE. X. THE
C PROGRAM STARTS THE SEARCH FOR ILEFT WITH THE VALUE OF ILEFT THAT WAS
C RETURNED AT THE PREVIOUS CALL (AND WAS SAVED IN THE LOCAL VARIABLE
C ILO) TO MINIMIZE THE WORK IN THE COMMON CASE THAT THE VALUE OF X ON
C THIS CALL IS CLOSE TO THE VALUE OF X ON THE PREVIOUS CALL. SHOULD
C THIS ASSUMPTION NOT BE VALID, THEN THE PROGRAM LOCATES ILO AND IHI
C SUCH THAT XT(ILO) .LE. X .LT. XT(IHI) AND, ONCE THEY ARE FOUND USES
C BISECTION TO FIND THE CORRECT VALUE FOR ILEFT. MFLAG IS AN ERROR FLAG.
C
C FOR DEFINITIONS OF CALLING ARGUMENTS SEE ABOVE AND BSPLVD.
C
C REFERENCE
C
C DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J.
C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472.
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. NONE
C CALLED BY.. COLPNT,INITAL,VALUES
C FORTRAN FUNCTIONS USED.. NONE
C-------------------------------------------------------------------------------
SAVE ILO
DIMENSION XT(LXT)
IF(MFLAG.EQ.-2) ILO = 1
IHI = ILO + 1
IF (IHI .LT. LXT) GO TO 20
IF (X .GE. XT(LXT)) GO TO 110
IF (LXT .LE. 1) GO TO 90
ILO = LXT - 1
GO TO 21
20 IF (X .GE. XT(IHI)) GO TO 40
21 IF (X .GE. XT(ILO)) GO TO 100
C-----------------------------------------------------------------------
C NOW X .LT. XT(IHI). FIND LOWER BOUND.
C-----------------------------------------------------------------------
30 ISTEP = 1
31 IHI = ILO
ILO = IHI - ISTEP
IF (ILO .LE. 1) GO TO 35
IF (X .GE. XT(ILO)) GO TO 50
ISTEP = ISTEP*2
GO TO 31
35 ILO = 1
IF (X .LT. XT(1)) GO TO 90
GO TO 50
C-----------------------------------------------------------------------
C NOW X .GE. XT(ILO). FIND UPPER BOUND.
C-----------------------------------------------------------------------
40 ISTEP = 1
41 ILO = IHI
IHI = ILO + ISTEP
IF (IHI .GE. LXT) GO TO 45
IF (X .LT. XT(IHI)) GO TO 50
ISTEP = ISTEP*2
GO TO 41
45 IF (X .GE. XT(LXT)) GO TO 110
IHI = LXT
C-----------------------------------------------------------------------
C NOW XT(ILO) .LE. X .LT. XT(IHI). NARROW THE INTERVAL.
C-----------------------------------------------------------------------
50 MIDDLE = (ILO + IHI)/2
IF (MIDDLE .EQ. ILO) GO TO 100
C-----------------------------------------------------------------------
C NOTE.. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1.
C-----------------------------------------------------------------------
IF (X .LT. XT(MIDDLE)) GO TO 53
ILO = MIDDLE
GO TO 50
53 IHI = MIDDLE
GO TO 50
C-----------------------------------------------------------------------
C SET OUTPUT AND RETURN.
C-----------------------------------------------------------------------
90 MFLAG = -1
ILEFT = 1
RETURN
100 MFLAG = 0
ILEFT = ILO
RETURN
110 MFLAG = 1
ILEFT = LXT
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE STIFIB (N0,Y,YMAX,ERROR,SAVE1,SAVE2,SAVE3,
* PW,IPIV,WORK,IWORK)
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C STIFIB PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM,
C A(Y,T)*(DY/DT) = G(Y,T), WHERE Y = (Y(1),Y(2), ... ,Y(N)).
C STIFIB IS FOR USE WHEN THE MATRICES A AND DG/DY HAVE BANDED OR NEARLY
C BANDED FORM. THE DEPENDENCE OF A(Y,T) ON Y IS ASSUMED TO BE WEAK.
C
C REFERENCE
C
C HINDMARSH, A.C., PRELIMINARY DOCUMENTATION OF GEARIB.. SOLUTION
C OF IMPLICIT SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS WITH
C BANDED JACOBIANS, LAWRENCE LIVERMORE LAB, UCID-30130, FEBRUARY
C 1976.
C
C COMMUNICATION WITH STIFIB IS DONE WITH THE FOLLOWING VARIABLES..
C
C Y AN N0 BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
C AND THEIR SCALED DERIVATIVES. LMAX IS 13 FOR THE ADAMS
C METHODS AND 6 FOR THE GEAR METHODS. LMAX - 1 = MAXDER
C IS THE MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET.
C Y(I,J+1) CONTAINS THE J-TH DERIVATIVE OF Y(I), SCALED BY
C H**J/FACTORIAL(J) (J = 0,1,...,NQ).
C N0 A CONSTANT INTEGER .GE. N, USED FOR DIMENSIONING PURPOSES.
C T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN.
C H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
C HMIN, THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEP SIZE
C HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY
C TIME, BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE.
C EPS THE RELATIVE ERROR BOUND. SEE DESCRIPTION IN PDECOL.
C UROUND THE UNIT ROUNDOFF OF THE MACHINE.
C N THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
C MF THE METHOD FLAG. SEE DESCRIPTION IN PDECOL.
C KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
C 0 THE STEP WAS SUCCESFUL.
C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED
C WITH DABS(H) = HMIN.
C -2 THE REQUESTED ERROR IS SMALLER THAN CAN
C BE HANDLED FOR THIS PROBLEM.
C -3 CORRECTOR CONVERGENCE COULD NOT BE
C ACHIEVED FOR DABS(H) = HMIN.
C -4 SINGULAR A-MATRIX ENCOUNTERED.
C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND
C THE Y ARRAY ARE AS OF THE BEGINNING OF THE LAST
C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
C JSTART AN INTEGER USED ON INPUT AND OUTPUT.
C ON INPUT, IT HAS THE FOLLOWING VALUES AND MEANINGS..
C 0 PERFORM THE FIRST STEP.
C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST.
C .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF
C H, EPS, N, AND/OR MF.
C ON EXIT, JSTART IS NQ, THE CURRENT ORDER OF THE METHOD.
C YMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL
C ERRORS IN Y ARE COMPARED.
C ERROR AN ARRAY OF N ELEMENTS. ERROR(I)/TQ(2) IS THE ESTIMATED
C ONE-STEP ERROR IN Y(I).
C SAVE1,SAVE2,SAVE3 THREE WORKING STORAGE ARRAYS, EACH OF LENGTH N.
C PW A BLOCK OF LOCATIONS USED FOR THE CHORD ITERATION
C MATRIX. SEE DESCRIPTION IN PDECOL.
C IPIV AN INTEGER ARRAY OF LENGTH N FOR PIVOT INFORMATION.
C ML,MU THE LOWER AND UPPER HALF BANDWIDTHS, RESPECTIVELY, OF
C THE CHORD ITERATION MATRIX. SEE DESCRIPTION IN PDECOL.
C WORK,IWORK WORKING ARRAYS WHICH ARE USED TO PASS APPROPRIATE
C ARRAYS TO OTHER SUBROUTINES.
C
C PACKAGE ROUTINES CALLED.. COSET,DIFFUN,PSETIB,RES,SOLB
C USER ROUTINES CALLED.. NONE
C CALLED BY.. PDECOL
C FORTRAN FUNCTIONS USED.. ABS,DMAX1,DMIN1,FLOAT
C-------------------------------------------------------------------------------
C SAVE EL,OLDL0,TQ,IER,NQ,L,METH,MITER
SAVE BND,CON,CRATE,D,D1,E,EDN,ENQ1,ENQ2,ENQ3,EPSOLD,
* EUP,FN,HOLD,OLDL0,PR1,PR2,PR3,R1,RC,RH,RMAX,TOLD,
* I,IDOUB,IER,IREDO,IRET,IWEVAL,J,J1,J2,L,LMAX,M,MEO,METH,
* MFOLD,MIO,MITER,NEWQ,NOLD,NQ,NSTEPJ
COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
* IW12,IW13,IW14,IW15,IW16,IW17,IW18
COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,N,MF,KFLAG,JSTART
COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W
COMMON /GEAR0/ HUSED,NQUSED,NSTEP,NFE,NJE
COMMON /OPTION/ NOGAUS,MAXDER
DIMENSION Y(NEQN,MAXDER+1),YMAX(NEQN),ERROR(NEQN),SAVE1(NEQN),
* SAVE2(NEQN),
* SAVE3(NEQN),PW(NEQN*(3*ML+1)),IPIV(NEQN),
* IWORK((NPDE+1)*NCPTS),
* WORK(KORD+NPDE*(4+9*NPDE)+(KORD+(NINT-1)*(KORD-NCC))*
* (3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4)))
DIMENSION EL(13),TQ(4)
DATA EL(2)/1./, OLDL0/1./, TQ(1)/0./, IER/0/
KFLAG = 0
TOLD = T
IF (JSTART .GT. 0) GO TO 200
IF (JSTART .NE. 0) GO TO 120
C-----------------------------------------------------------------------
C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT IS
C CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE
C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
C FOR THE NEXT INCREASE.
C-----------------------------------------------------------------------
NQ = 1
IER = 0
CALL DIFFUN (N, T, Y, SAVE1, IER, PW, IPIV, WORK, IWORK)
IF ( IER .NE. 0 ) GO TO 685
DO 110 I = 1,N
110 Y(I,2) = H*SAVE1(I)
METH = MF/10
MITER = MF - 10*METH
L = 2
IDOUB = 3
RMAX = 1.E+04
RC = 0.
CRATE = 1.
EPSOLD = EPS
HOLD = H
MFOLD = MF
NOLD = N
NSTEP = 0
NSTEPJ = 0
NFE = 0
NJE = 1
IRET = 3
GO TO 130
C-----------------------------------------------------------------------
C IF THE CALLER HAS CHANGED METH, COSET IS CALLED TO SET
C THE COEFFICIENTS OF THE METHOD. IF THE CALLER HAS CHANGED
C N, EPS, OR METH, THE CONSTANTS E, EDN, EUP, AND BND MUST BE RESET.
C E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS
C TO TEST FOR INCREASING THE ORDER, EDN FOR DECREASING THE ORDER.
C BND IS USED TO TEST FOR CONVERGENCE OF THE CORRECTOR ITERATES.
C IF THE CALLER HAS CHANGED H, Y MUST BE RESCALED.
C IF H OR METH HAS BEEN CHANGED, IDOUB IS RESET TO L + 1 TO PREVENT
C FURTHER CHANGES IN H FOR THAT MANY STEPS.
C-----------------------------------------------------------------------
120 IF (MF .EQ. MFOLD) GO TO 150
MEO = METH
MIO = MITER
METH = MF/10
MITER = MF - 10*METH
MFOLD = MF
IF (MITER .NE. MIO) IWEVAL = MITER
IF (METH .EQ. MEO) GO TO 150
IDOUB = L + 1
IRET = 1
130 CALL COSET (METH, NQ, EL, TQ)
LMAX = MAXDER + 1
RC = RC*EL(1)/OLDL0
OLDL0 = EL(1)
140 FN = FLOAT(N)
EDN = FN*(TQ(1)*EPS)**2
E = FN*(TQ(2)*EPS)**2
EUP = FN*(TQ(3)*EPS)**2
BND = FN*(TQ(4)*EPS)**2
GO TO (160, 170, 200), IRET
150 IF ((EPS .EQ. EPSOLD) .AND. (N .EQ. NOLD)) GO TO 160
EPSOLD = EPS
NOLD = N
IRET = 1
GO TO 140
160 IF (H .EQ. HOLD) GO TO 200
RH = H/HOLD
H = HOLD
IREDO = 3
GO TO 175
170 RH = DMAX1(RH,HMIN/ DABS(H))
175 RH = DMIN1(RH,HMAX/ DABS(H),RMAX)
R1 = 1.
DO 180 J = 2,L
R1 = R1*RH
DO 180 I = 1,N
180 Y(I,J) = Y(I,J)*R1
H = H*RH
RC = RC*RH
IDOUB = L + 1
IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
C MULTIPLYING THE Y ARRAY BY THE PASCAL TRIANGLE MATRIX.
C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1).
C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, OR THE CALLER HAS
C CHANGED MITER, IWEVAL IS SET TO MITER TO FORCE PW TO BE UPDATED.
C IN ANY CASE, PW IS UPDATED AT LEAST EVERY 40-TH STEP.
C PW IS THE CHORD ITERATION MATRIX A - H*EL(1)*(DG/DY).
C-----------------------------------------------------------------------
200 IF ( DABS(RC-1.) .GT. 0.3) IWEVAL = MITER
IF (NSTEP .GE. NSTEPJ+40) IWEVAL = MITER
T = T + H
DO 210 J1 = 1,NQ
DO 210 J2 = J1,NQ
J = (NQ + J1) - J2
DO 210 I = 1,N
210 Y(I,J) = Y(I,J) + Y(I,J+1)
C-----------------------------------------------------------------------
C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS
C MADE ON THE R.M.S. NORM OF EACH CORRECTION, USING BND, WHICH
C IS DEPENDENT ON EPS. THE SUM OF THE CORRECTIONS IS ACCUMULATED
C IN THE VECTOR ERROR(I). THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR
C LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1.
C THE UPDATED H*YDOT IS STORED IN SAVE2.
C-----------------------------------------------------------------------
220 DO 230 I = 1,N
SAVE2(I) = Y(I,2)
230 ERROR(I) = 0.
M = 0
CALL RES (T, H, Y, SAVE2, SAVE3, NPDE, NCPTS, WORK(IW1), IWORK,
* WORK, WORK(IW14), WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9))
NFE = NFE + 1
IF (IWEVAL .LE. 0) GO TO 350
C-----------------------------------------------------------------------
C IF INDICATED, THE MATRIX PW IS REEVALUATED BEFORE STARTING THE
C CORRECTOR ITERATION. IWEVAL IS SET TO 0 AS AN INDICATOR
C THAT THIS HAS BEEN DONE. PW IS COMPUTED AND PROCESSED IN PSETIB.
C-----------------------------------------------------------------------
IWEVAL = 0
RC = 1.0D0
NJE = NJE + 1
NSTEPJ = NSTEP
CON = -H*EL(1)
CALL PSETIB (Y, PW, N0, CON, MITER, IER, WORK(IW1), IWORK,
* WORK(IW3),WORK(IW9),SAVE2,IPIV,YMAX,WORK(IW11),WORK(IW12),
* WORK(IW13),WORK(IW16),WORK(IW14),WORK(IW15),WORK,NPDE)
IF (IER .NE. 0) GO TO 420
C-----------------------------------------------------------------------
C COMPUTE THE CORRECTOR ERROR, R SUB M, AND SOLVE THE LINEAR SYSTEM
C WITH THAT AS RIGHT-HAND SIDE AND PW AS COEFFICIENT MATRIX,
C USING THE LU DECOMPOSITION OF PW.
C-----------------------------------------------------------------------
350 CALL SOLB (N0, N, ML, MU, PW, SAVE3, IPIV)
370 D = 0.0D0
DO 380 I = 1,N
ERROR(I) = ERROR(I) + SAVE3(I)
D = D + (SAVE3(I)/YMAX(I))**2
SAVE1(I) = Y(I,1) + EL(1)*ERROR(I)
380 SAVE2(I) = Y(I,2) + ERROR(I)
C-----------------------------------------------------------------------
C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
C-----------------------------------------------------------------------
400 IF (M .NE. 0) CRATE = DMAX1(.9*CRATE,D/D1)
IF ((D*DMIN1(1.D0,2.0D0*CRATE)) .LE. BND) GO TO 450
D1 = D
M = M + 1
IF (M .EQ. 3) GO TO 410
CALL RES(T, H, SAVE1, SAVE2, SAVE3, NPDE, NCPTS, WORK(IW1), IWORK,
* WORK, WORK(IW14), WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9))
GO TO 350
C-----------------------------------------------------------------------
C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES.
C IF THE MATRIX PW IS NOT UP TO DATE, IT IS REEVALUATED FOR THE
C NEXT TRY. OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES
C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF NOT, A
C NO-CONVERGENCE EXIT IS TAKEN.
C-----------------------------------------------------------------------
410 NFE = NFE + 2
IF (IWEVAL .EQ. -1) GO TO 440
420 T = TOLD
RMAX = 2.
DO 430 J1 = 1,NQ
DO 430 J2 = J1,NQ
J = (NQ + J1) - J2
DO 430 I = 1,N
430 Y(I,J) = Y(I,J) - Y(I,J+1)
IF ( DABS(H) .LE. HMIN*1.00001) GO TO 680
RH = .25
IREDO = 1
GO TO 170
440 IWEVAL = MITER
GO TO 220
C-----------------------------------------------------------------------
C THE CORRECTOR HAS CONVERGED. IWEVAL IS SET TO -1 TO SIGNAL
C THAT PW MAY NEED UPDATING ON SUBSEQUENT STEPS. THE ERROR TEST
C IS MADE AND CONTROL PASSES TO STATEMENT 500 IF IT FAILS.
C-----------------------------------------------------------------------
450 IWEVAL = -1
NFE = NFE + M
D = 0.
DO 460 I = 1,N
460 D = D + (ERROR(I)/YMAX(I))**2
IF (D .GT. E) GO TO 500
C-----------------------------------------------------------------------
C AFTER A SUCCESSFUL STEP, UPDATE THE Y ARRAY.
C CONSIDER CHANGING H IF IDOUB = 1. OTHERWISE DECREASE IDOUB BY 1.
C IF IDOUB IS THEN 1 AND NQ .LT. MAXDER, THEN ERROR IS SAVED FOR
C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A
C FACTOR OF AT LEAST 1.1. IF NOT, IDOUB IS SET TO 10 TO PREVENT
C TESTING FOR THAT MANY STEPS.
C-----------------------------------------------------------------------
KFLAG = 0
IREDO = 0
NSTEP = NSTEP + 1
HUSED = H
NQUSED = NQ
DO 470 J = 1,L
DO 470 I = 1,N
470 Y(I,J) = Y(I,J) + EL(J)*ERROR(I)
IF (IDOUB .EQ. 1) GO TO 520
IDOUB = IDOUB - 1
IF (IDOUB .GT. 1) GO TO 700
IF (L .EQ. LMAX) GO TO 700
DO 490 I = 1,N
490 Y(I,LMAX) = ERROR(I)
GO TO 700
C-----------------------------------------------------------------------
C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
C RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
C ONE LOWER ORDER.
C-----------------------------------------------------------------------
500 KFLAG = KFLAG - 1
T = TOLD
DO 510 J1 = 1,NQ
DO 510 J2 = J1,NQ
J = (NQ + J1) - J2
DO 510 I = 1,N
510 Y(I,J) = Y(I,J) - Y(I,J+1)
RMAX = 2.
IF ( DABS(H) .LE. HMIN*1.00001) GO TO 660
IF (KFLAG .LE. -3) GO TO 640
IREDO = 2
PR3 = 1.E+20
GO TO 540
C-----------------------------------------------------------------------
C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
C PR1, PR2, AND PR3 ARE COMPUTED, BY WHICH H COULD BE DIVIDED
C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
C IN THE CASE OF FAILURE, PR3 = 1.E20 TO AVOID AN ORDER INCREASE.
C THE SMALLEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
C ADDITIONAL SCALED DERIVATIVE.
C-----------------------------------------------------------------------
520 PR3 = 1.E+20
IF (L .EQ. LMAX) GO TO 540
D1 = 0.
DO 530 I = 1,N
530 D1 = D1 + ((ERROR(I) - Y(I,LMAX))/YMAX(I))**2
ENQ3 = .5/ FLOAT(L+1)
PR3 = ((D1/EUP)**ENQ3)*1.4 + 1.4D-06
540 ENQ2 = .5/ FLOAT(L)
PR2 = ((D/E)**ENQ2)*1.2 + 1.2D-06
PR1 = 1.E+20
IF (NQ .EQ. 1) GO TO 560
D = 0.
DO 550 I = 1,N
550 D = D + (Y(I,L)/YMAX(I))**2
ENQ1 = .5/ FLOAT(NQ)
PR1 = ((D/EDN)**ENQ1)*1.3 + 1.3D-06
560 IF (PR2 .LE. PR3) GO TO 570
IF (PR3 .LT. PR1) GO TO 590
GO TO 580
570 IF (PR2 .GT. PR1) GO TO 580
NEWQ = NQ
RH = 1./PR2
GO TO 620
580 NEWQ = NQ - 1
RH = 1./PR1
GO TO 620
590 NEWQ = L
RH = 1./PR3
IF (RH .LT. 1.1) GO TO 610
DO 600 I = 1,N
600 Y(I,NEWQ+1) = ERROR(I)*EL(L)/ FLOAT(L)
GO TO 630
610 IDOUB = 10
GO TO 700
620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1)) GO TO 610
C-----------------------------------------------------------------------
C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
C IN ANY CASE H IS RESET ACCORDING TO RH AND THE Y ARRAY IS RESCALED.
C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
C-----------------------------------------------------------------------
IF (NEWQ .EQ. NQ) GO TO 170
630 NQ = NEWQ
L = NQ + 1
IRET = 2
GO TO 130
C-----------------------------------------------------------------------
C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED.
C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
C Y ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN
C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED.
C AFTER A TOTAL OF 7 FAILURES, AN EXIT IS TAKEN WITH KFLAG = -2.
C-----------------------------------------------------------------------
640 IF (KFLAG .EQ. -7) GO TO 670
RH = .1
RH = DMAX1(HMIN/ DABS(H),RH)
H = H*RH
IER = 0
CALL DIFFUN (N, T, Y, SAVE1, IER, PW, IPIV, WORK, IWORK)
IF ( IER .NE. 0 ) GO TO 685
NJE = NJE + 1
DO 650 I = 1,N
650 Y(I,2) = H*SAVE1(I)
IWEVAL = MITER
IDOUB = 10
IF (NQ .EQ. 1) GO TO 200
NQ = 1
L = 2
IRET = 3
GO TO 130
C-----------------------------------------------------------------------
C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD
C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
C-----------------------------------------------------------------------
660 KFLAG = -1
GO TO 700
670 KFLAG = -2
GO TO 700
680 KFLAG = -3
GO TO 700
685 KFLAG = -4
GO TO 700
690 RMAX = 10.
700 HOLD = H
JSTART = NQ
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE GFUN ( T,C,UDOT,NPDE,NCPTS,A,BC,DBDU,DBDUX,DZDT,
* XC,UVAL,ILEFT )
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL.
C
C SUBROUTINE GFUN COMPUTES THE FUNCTION UDOT=G(C,T), THE RIGHT-
C HAND SIDE OF THE SEMI-DISCRETE APPROXIMATION TO THE ORIGINAL
C SYSTEM OF PARTIAL DIFFERENTIAL EQUATIONS AND UPDATES THE BOUNDARY
C CONDITION INFORMATION.
C
C PACKAGE ROUTINES CALLED.. EVAL
C USER ROUTINES CALLED.. BNDRY,F
C CALLED BY.. DIFFUN,PSETIB,RES
C FORTRAN FUNCTIONS USED.. NONE
C-------------------------------------------------------------------------------
COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPT,NEQN,IQUAD
DIMENSION C(NPDE,NCPTS),UDOT(NPDE,NCPTS)
DIMENSION A(NCPTS*3*KORD),BC(NPDE,NPDE,4),
* XC(NCPTS),UVAL(NPDE,3),ILEFT(NCPTS)
DIMENSION DZDT(NPDE),DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE)
DO 10 K=1,4
DO 10 J=1,NPDE
DO 10 I=1,NPDE
BC(I,J,K) = 0.0
10 CONTINUE
C-----------------------------------------------------------------------
C UPDATE THE LEFT BOUNDARY VALUES. SAVE LEFT BOUNDARY CONDITION
C INFORMATION IN THE FIRST 2*NPDE*NPDE LOCATIONS OF BC.
C
C NOTE.. UVAL(K,1) = U(K), UVAL(K,2) = UX(K), AND UVAL(K,3) = UXX(K).
C-----------------------------------------------------------------------
CALL EVAL(1,NPDE,C,UVAL,A,ILEFT)
CALL BNDRY(T,XC(1),UVAL,UVAL(1,2),DBDU,DBDUX,DZDT,NPDE)
CALL F(T,XC(1),UVAL,UVAL(1,2),UVAL(1,3),UDOT,NPDE)
ILIM = KORD + 2
DO 30 K=1,NPDE
BC(K,K,1) = 1.
IF( DBDU(K,K) .EQ. 0.0 .AND. DBDUX(K,K) .EQ. 0.0 ) GO TO 30
UDOT(K,1) = DZDT(K)
DO 20 J=1,NPDE
BC(K,J,2) = A(ILIM) * DBDUX(K,J)
BC(K,J,1) = DBDU(K,J) - BC(K,J,2)
20 CONTINUE
30 CONTINUE
C-----------------------------------------------------------------------
C MAIN LOOP TO FORM RIGHT SIDE OF ODES AT THE COLLOCATION POINTS.
C-----------------------------------------------------------------------
ILIM = NCPTS - 1
DO 40 I=2,ILIM
CALL EVAL(I,NPDE,C,UVAL,A,ILEFT)
CALL F(T,XC(I),UVAL,UVAL(1,2),UVAL(1,3),UDOT(1,I),NPDE)
40 CONTINUE
C-----------------------------------------------------------------------
C UPDATE THE RIGHT BOUNDARY VALUES. SAVE THE RIGHT BOUNDARY CONDITION
C INFORMATION IN THE LAST 2*NPDE*NPDE LOCATIONS IN BC.
C-----------------------------------------------------------------------
CALL EVAL(NCPTS,NPDE,C,UVAL,A,ILEFT)
CALL F(T,XC(NCPTS),UVAL,UVAL(1,2),UVAL(1,3),UDOT(1,NCPTS),NPDE)
CALL BNDRY(T,XC(NCPTS),UVAL,UVAL(1,2),DBDU,DBDUX,DZDT,NPDE)
ILIM = NCPTS * 3 * KORD - KORD - 1
DO 60 K=1,NPDE
BC(K,K,4) = 1.
IF( DBDU(K,K) .EQ. 0.0 .AND. DBDUX(K,K) .EQ. 0.0 ) GO TO 60
UDOT(K,NCPTS) = DZDT(K)
DO 50 J=1,NPDE
BC(K,J,3) = A(ILIM) * DBDUX(K,J)
BC(K,J,4) = DBDU(K,J) - BC(K,J,3)
50 CONTINUE
60 CONTINUE
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE EVAL( ICPT,NPDE,C,UVAL,A,ILEFT )
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL.
C
C SUBROUTINE EVAL EVALUATES U(K), UX(K), AND UXX(K), K=1 TO NPDE,
C AT THE COLLOCATION POINT WITH INDEX ICPT USING THE VALUES OF
C THE BASIS FUNCTION COEFFICIENTS IN C AND THE BASIS FUNCTION VALUES
C STORED IN A. THE RESULTS ARE STORED IN UVAL AS FOLLOWS..
C UVAL(K,1) = U(K), UVAL(K,2) = UX(K), AND UVAL(K,3) = UXX(K).
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. NONE
C CALLED BY.. GFUN,PDECOL,PSETIB
C FORTRAN FUNCTIONS USED.. NONE
C-------------------------------------------------------------------------------
COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD
DIMENSION C(NPDE,NCPTS),UVAL(NPDE,3),A(NCPTS*3*KORD),ILEFT(NCPTS)
IK = ILEFT(ICPT) - KORD
IC = 3*KORD*(ICPT-1)
DO 10 M=1,3
ICC = IC + KORD*(M-1)
DO 10 J=1,NPDE
UVAL(J,M) = 0.
DO 10 I=1,KORD
UVAL(J,M) = UVAL(J,M) + C(J,I+IK)*A(I+ICC)
10 CONTINUE
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE DIFFUN( N, T, Y, YDOT, IER, PW, IPIV, WORK, IWORK )
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL.
C
C THIS ROUTINE COMPUTES YDOT = A(Y,T)**-1 * G(Y,T) BY USE OF
C THE ROUTINES GFUN, ADDA, DECB, AND SOLB.
C
C PACKAGE ROUTINES CALLED.. ADDA,DECB,GFUN,SOLB
C USER ROUTINES CALLED.. NONE
C CALLED BY.. STIFIB
C FORTRAN FUNCTIONS USED.. NONE
C-------------------------------------------------------------------------------
COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W
COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
COMMON /ISTART/ IW1,IW2,IW3,IDUM(5),IW9,IW10,IW11,IW12,IW13,IW14,
* IW15,IW16,IW17,IW18
COMMON /OPTION/ NOGAUS,MAXDER
DIMENSION Y(NEQN),YDOT(NEQN),PW(NEQN*(3*ML+1)),
* IPIV(NEQN),IWORK((NPDE+1)*NCPTS),
* WORK(KORD+NPDE*(4+9*NPDE)+(KORD+(NINT-1)*(KORD-NCC))*
* (3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4)))
CALL GFUN (T, Y, YDOT, NPDE, NCPTS, WORK(IW1), WORK, WORK(IW14),
* WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9), IWORK)
DO 10 I = 1,N0W
10 PW(I) = 0.
N0 = NM1 + 1
CALL ADDA (PW, N0, WORK(IW1), IWORK, WORK, NPDE)
CALL DECB (N0, N, ML, MU, PW, IPIV, IER)
IF ( IER .NE. 0 ) RETURN
CALL SOLB (N0, N, ML, MU, PW, YDOT, IPIV)
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE ADDA( PW,N0,A,ILEFT,BC,NPDE )
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL AND STIFIB.
C
C SUBROUTINE ADDA ADDS THE MATRIX A TO THE MATRIX STORED IN PW IN
C BAND FORM. PW IS STORED BY DIAGONALS WITH THE LOWERMOST DIAGONAL
C STORED IN THE FIRST COLUMN OF THE ARRAY.
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. NONE
C CALLED BY.. DIFFUN,PSETIB
C FORTRAN FUNCTIONS USED.. NONE
C-------------------------------------------------------------------------------
COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD
COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W
DIMENSION PW(NEQN,3*ML+1),A(3*KORD*NCPTS),
* ILEFT(NCPTS),BC(NPDE,NPDE,4)
C-----------------------------------------------------------------------
C ADD THE BOUNDARY CONDITION PORTIONS OF THE A MATRIX TO PW ( THE FIRST
C AND LAST BLOCK ROWS).
C-----------------------------------------------------------------------
ICOL = (ILEFT(1) + IQUAD - 1) * NPDE
DO 10 I=1,NPDE
IBOT = NEQN - NPDE + I
DO 10 J=1,NPDE
IND = ICOL + J - I
PW(I,IND) = PW(I,IND) + BC(I,J,1)
PW(I,IND+NPDE) = PW(I,IND+NPDE) + BC(I,J,2)
PW(IBOT,IND-NPDE) = PW(IBOT,IND-NPDE) + BC(I,J,3)
PW(IBOT,IND) = PW(IBOT,IND) + BC(I,J,4)
10 CONTINUE
C-----------------------------------------------------------------------
C UPDATE THE REMAINING ROWS OF PW BY ADDING THE APPROPRIATE VALUES
C IN A TO PW.
C-----------------------------------------------------------------------
IND = NCPTS - 1
DO 20 I=2,IND
I1 = (I-1) * NPDE
I2 = (I-1) * KORD * 3
ICOL = ILEFT(I) - I + IQUAD - 1
DO 20 J=1,KORD
J1 = (ICOL+J) * NPDE
J2 = I2 + J
DO 20 JJ=1,NPDE
20 PW(I1+JJ,J1) = PW(I1+JJ,J1) + A(J2)
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE RES( T,H,C,V,R,NPDE,NCPTS,A,ILEFT,BC,DBDU,DBDUX,DZDT,
* XC,UVAL )
IMPLICIT REAL*8 (A-H, O-Z)
C-------------------------------------------------------------------------------
C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL.
C
C SUBROUTINE RES COMPUTES THE RESIDUAL VECTOR R = H*G(C,T) - A(C,T)*V
C WHERE H IS THE CURRENT TIME STEP SIZE, G IS A VECTOR, A IS A
C MATRIX, V IS A VECTOR, AND T IS THE CURRENT TIME.
C
C PACKAGE ROUTINES CALLED.. GFUN
C USER ROUTINES CALLED.. NONE
C CALLED BY.. STIFIB
C FORTRAN FUNCTIONS USED.. NONE
C-------------------------------------------------------------------------------
SAVE
COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPT,NEQN,IQUAD
DIMENSION C(NPDE,NCPTS),R(NPDE,NCPTS),V(NPDE,NCPTS)
DIMENSION A(3*KORD*NCPTS),ILEFT(NCPTS),BC(NPDE,NPDE,4),XC(NCPTS),
* UVAL(3*NPDE)
DIMENSION DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE),DZDT(NPDE)
C-----------------------------------------------------------------------
C FORM G(C,T) AND STORE IN R.
C-----------------------------------------------------------------------
CALL GFUN(T,C,R,NPDE,NCPTS,A,BC,DBDU,DBDUX,DZDT,XC,UVAL,ILEFT)
C-----------------------------------------------------------------------
C FORM THE FIRST AND LAST BLOCK ROWS OF THE RESIDUAL VECTOR
C WHICH ARE DEPENDENT ON THE BOUNDARY CONDITIONS.
C-----------------------------------------------------------------------
ILIM = NCPTS - 1
DO 20 I=1,NPDE
SUM1 = 0.0
SUM2 = 0.0
DO 10 J=1,NPDE
SUM1 = SUM1 + BC(I,J,1) * V(J,1) + BC(I,J,2) * V(J,2)
SUM2 = SUM2 + BC(I,J,3) * V(J,ILIM) + BC(I,J,4) * V(J,NCPTS)
10 CONTINUE
R(I,1) = H * R(I,1) - SUM1
R(I,NCPTS) = H * R(I,NCPTS) - SUM2
20 CONTINUE
C-----------------------------------------------------------------------
C FORM THE REMAINING COMPONENTS OF THE RESIDUAL VECTOR.
C-----------------------------------------------------------------------
DO 50 ICPTS=2,ILIM
I2 = (ICPTS-1) * KORD * 3
ICOL = ILEFT(ICPTS) - KORD
DO 40 JJ=1,NPDE
SUM1 = 0.
DO 30 J=1,KORD
SUM1 = SUM1 + A(I2+J) * V(JJ,ICOL+J)
30 CONTINUE
R(JJ,ICPTS) = H*R(JJ,ICPTS) - SUM1
40 CONTINUE
50 CONTINUE
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE PSETIB( C,PW,N0,CON,MITER,IER,A,ILEFT,XC,UVAL,
* SAVE2,IPIV,CMAX,DFDU,DFDUX,DFDUXX,DZDT,DBDU,DBDUX,BC,NPDE )
IMPLICIT REAL*8 (A-H, O-Z)
C-----------------------------------------------------------------------
C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL AND STIFIB.
C
C PSETIB IS CALLED BY STIFIB TO COMPUTE AND PROCESS THE MATRIX
C PW = A - H*EL(1)*(DG/DC), WHERE A AND DG/DC ARE TREATED IN BAND
C FORM. DG/DC IS COMPUTED, EITHER WITH THE AID OF THE USER-SUPPLIED
C ROUTINE DERIVF IF MITER = 1, OR BY FINITE DIFFERENCING WITH THE AID
C OF THE PACKAGE-SUPPLIED ROUTINE DIFFF IF MITER = 2. FINALLY,
C PW IS SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER
C SOLUTION OF LINEAR SYSTEMS WITH PW AS COEFFICIENT MATRIX.
C SEE SUBROUTINES DECB AND SOLB.
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C WITH PSETIB USES THE FOLLOWING..
C EPSJ = DSQRT(UROUND), USED IN THE NUMERICAL JACOBIAN INCREMENTS.
C MW = ML + MU + 1.
C NM1 = N0 - 1.
C N0ML = N0*ML.
C N0W = N0*MW.
C
C PACKAGE ROUTINES CALLED.. ADDA,DECB,DIFFF,EVAL,GFUN
C USER ROUTINES CALLED.. BNDRY,DERIVF
C CALLED BY.. STIFIB
C FORTRAN FUNCTIONS USED.. ABS,FLOAT,MAX0,MIN0,DSQRT
C-----------------------------------------------------------------------
COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD
COMMON /GEAR1/ T,H,DUMMY(3),UROUND,N,IDUMMY(3)
COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W
DIMENSION PW(NEQN,3*ML+1),C(NEQN),CMAX(NEQN)
DIMENSION A(3*KORD*NCPTS),ILEFT(NCPTS),BC(4*NPDE*NPDE),
* XC(NCPTS),UVAL(NPDE,3),SAVE2(NEQN),IPIV(NEQN)
DIMENSION DFDU(NPDE,NPDE),DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE)
DIMENSION DZDT(NPDE),DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE)
DO 10 I=1,NEQN
DO 5 J=1,MW
5 PW(I,J)=0.0E0
10 CONTINUE
IF ( MITER .EQ. 1 ) GO TO 25
CALL GFUN (T, C, SAVE2, NPDE, NCPTS,A,BC,DBDU,DBDUX,DZDT,XC,
* UVAL,ILEFT)
D = 0.
DO 20 I = 1,N
20 D = D + SAVE2(I)**2
R0 = DABS(H)* DSQRT(D/FLOAT(N0))*1.E+03*UROUND
C-----------------------------------------------------------------------
C COMPUTE BLOCK ROWS OF JACOBIAN.
C-----------------------------------------------------------------------
25 DO 30 I=1,NCPTS
I1 = (I-1)*NPDE
I2 = (I-1)*KORD*3
CALL EVAL(I,NPDE,C,UVAL,A,ILEFT)
IF ( MITER .EQ. 1 )
* CALL DERIVF(T,XC(I),UVAL,UVAL(1,2),UVAL(1,3),
* DFDU,DFDUX,DFDUXX,NPDE)
IF ( MITER .EQ. 2 )
* CALL DIFFF(T,XC(I),I,UVAL,UVAL(1,2),UVAL(1,3),
* DFDU,DFDUX,DFDUXX,NPDE,CMAX,SAVE2)
ICOL = ILEFT(I) - I + IQUAD - 1
KLOW = MAX0(1,I+2-NCPTS)
KUP = MIN0(KORD,KORD+I-2)
DO 30 KBLK=KLOW,KUP
J1 = (ICOL+KBLK)*NPDE
J2 = I2 + KBLK
J3 = J2 + KORD
J4 = J3 + KORD
DO 30 L=1,NPDE
DO 30 K=1,NPDE
PW(I1+K,J1-K+L) = DFDU(K,L)*A(J2) + DFDUX(K,L)*A(J3)
* + DFDUXX(K,L)*A(J4)
30 CONTINUE
C-----------------------------------------------------------------------
C MODIFY THE LAST AND THE FIRST BLOCK ROWS FOR THE BOUNDARY CONDITIONS.
C CURRENT INFORMATION FOR THE RIGHT BOUNDARY CONDITION IS ALREADY IN
C THE ARRAYS DBDU, DBDUX AS A RESULT OF A PREVIOUS CALL TO GFUN.
C-----------------------------------------------------------------------
IROW = NEQN - NPDE
DO 50 K=1,NPDE
IROW = IROW + 1
IF(DBDU(K,K) .EQ. 0.0 .AND. DBDUX(K,K) .EQ. 0.0) GO TO 50
DO 40 J=1,MW
PW(IROW,J) = 0.0
40 CONTINUE
50 CONTINUE
CALL EVAL(1,NPDE,C,UVAL,A,ILEFT)
CALL BNDRY(T,XC(1),UVAL,UVAL(1,2),DBDU,DBDUX,DZDT,NPDE)
DO 70 K=1,NPDE
IF(DBDU(K,K) .EQ. 0.0 .AND. DBDUX(K,K) .EQ. 0.0) GO TO 70
DO 60 J=1,MW
PW(K,J) = 0.0
60 CONTINUE
70 CONTINUE
DO 80 I = 1,N0
DO 85 J=1,MW
85 PW(I,J)=PW(I,J)*CON
80 CONTINUE
C-----------------------------------------------------------------------
C ADD MATRIX A(C,T) TO PW.
C-----------------------------------------------------------------------
CALL ADDA (PW, N0, A, ILEFT, BC, NPDE)
C-----------------------------------------------------------------------
C DO LU DECOMPOSITION ON PW.
C-----------------------------------------------------------------------
CALL DECB (N0, N, ML, MU, PW, IPIV, IER)
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE DIFFF( T,X,IPT,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE,CMAX,
* SAVE2)
IMPLICIT REAL*8 (A-H, O-Z)
C-----------------------------------------------------------------------
C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL.
C
C SUBROUTINE DIFFF IS USED IF MITER=2 TO PROVIDE FINITE DIFFERENCE
C APPROXIMATIONS FOR THE PARTIAL DERIVATIVES OF THE K-TH USER DEFINED
C FUNCTION IN THE F ROUTINE WITH RESPECT TO THE VARIABLES U, UX, AND
C UXX. THESE PARTIALS WITH RESPECT TO U, UX, AND UXX ARE COMPUTED,
C STORED, AND RETURNED IN THE NPDE BY NPDE ARRAYS DFDU, DFDUX, AND
C DFDUXX, RESPECTIVELY, AT COLLOCATION POINT NUMBER IPT.
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. F
C CALLED BY.. PSETIB
C FORTRAN FUNCTIONS USED.. DMAX1
C-----------------------------------------------------------------------
COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W
COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD
DIMENSION U(NPDE),UX(NPDE),UXX(NPDE),DFDU(NPDE,NPDE),
* DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE),CMAX(NEQN),SAVE2(NEQN)
ID = (IPT-1) * NPDE
DO 40 J=1,NPDE
UJ = U(J)
R = EPSJ * CMAX(J)
R = DMAX1(R,R0)
U(J) = U(J) + R
RINV = 1. / R
CALL F(T,X,U,UX,UXX,DFDU(1,J),NPDE)
DO 10 I=1,NPDE
10 DFDU(I,J) = ( DFDU(I,J) - SAVE2(I+ID) ) * RINV
U(J) = UJ
UJ = UX(J)
UX(J) = UX(J) + R
CALL F(T,X,U,UX,UXX,DFDUX(1,J),NPDE)
DO 20 I=1,NPDE
20 DFDUX(I,J) = ( DFDUX(I,J) - SAVE2(I+ID) ) * RINV
UX(J) = UJ
UJ = UXX(J)
UXX(J) = UXX(J) + R
CALL F(T,X,U,UX,UXX,DFDUXX(1,J),NPDE)
DO 30 I=1,NPDE
30 DFDUXX(I,J) = ( DFDUXX(I,J) - SAVE2(I+ID) ) * RINV
UXX(J) = UJ
40 CONTINUE
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE INTERP ( TOUT, Y, N0, Y0 )
IMPLICIT REAL*8 (A-H, O-Z)
C-----------------------------------------------------------------------
C CALLING ARGUMENTS ARE DEFINED BELOW AND IN STIFIB
C
C SUBROUTINE INTERP COMPUTES INTERPOLATED VALUES OF THE DEPENDENT
C VARIABLE Y AND STORES THEM IN Y0. THE INTERPOLATION IS TO THE
C POINT T = TOUT, AND USES THE NORDSIECK HISTORY ARRAY Y, AS FOLLOWS..
C NQ
C Y0(I) = SUM Y(I,J+1)*S**J ,
C J=0
C WHERE S = -(T-TOUT)/H.
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. NONE
C CALLED BY.. PDECOL
C FORTRAN FUNCTIONS USED.. NONE
C-----------------------------------------------------------------------
COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD
COMMON /OPTION/ NOGAUS,MAXDER
COMMON /GEAR1/ T,H,DUMMY(4),N,IDUMMY(2),JSTART
DIMENSION Y0(NEQN),Y(NEQN,MAXDER+1)
DO 10 I = 1,N
10 Y0(I) = Y(I,1)
L = JSTART + 1
S = (TOUT - T)/H
S1 = 1.
DO 30 J = 2,L
S1 = S1*S
DO 20 I = 1,N
20 Y0(I) = Y0(I) + S1*Y(I,J)
30 CONTINUE
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE COSET ( METH, NQ, EL, TQ )
IMPLICIT REAL*8 (A-H, O-Z)
C-----------------------------------------------------------------------
C COSET IS CALLED BY THE INTEGRATOR AND SETS COEFFICIENTS USED THERE.
C THE VECTOR EL, OF LENGTH NQ + 1, DETERMINES THE BASIC METHOD.
C THE VECTOR TQ, OF LENGTH 4, IS INVOLVED IN ADJUSTING THE STEP SIZE
C IN RELATION TO TRUNCATION ERROR. ITS VALUES ARE GIVEN BY THE
C PERTST ARRAY.
C
C THE VECTORS EL AND TQ DEPEND ON METH AND NQ.
C THE MAXIMUM ORDER, MAXDER, OF THE METHODS AVAILABLE IS CURRENTLY
C 12 FOR THE ADAMS METHODS AND 5 FOR THE BDF METHODS. MAXDER DEFAULTS
C TO 5 UNLESS THE USER SETS MAXDER TO SOME OTHER LEGITIMATE VALUE
C THROUGH THE COMMON BLOCK /OPTION/. SEE PDECOL FOR ADDITIONAL DETAILS.
C LMAX = MAXDER + 1 IS THE NUMBER OF COLUMNS IN THE Y ARRAY (SEE STIFIB
C AND THE VARIABLE C, Y, OR WORK(IW10) IN PDECOL.
C
C THE COEFFICIENTS IN PERTST NEED BE GIVEN TO ONLY ABOUT
C ONE PERCENT ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW
C IS.. COEFFICIENTS FOR ORDER NQ - 1, COEFFICIENTS FOR ORDER NQ,
C COEFFICIENTS FOR ORDER NQ + 1. WITHIN EACH GROUP ARE THE
C COEFFICIENTS FOR THE ADAMS METHODS, FOLLOWED BY THOSE FOR THE
C BDF METHODS.
C
C REFERENCE
C
C GEAR, C.W., NUMERICAL INITIAL VALUE PROBLEMS IN ORDINARY
C DIFFERENTIAL EQUATIONS, PRENTICE-HALL, ENGLEWOOD CLIFFS,
C N. J., 1971.
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. NONE
C CALLED BY.. STIFIB
C FORTRAN FUNCTIONS USED.. FLOAT
C-----------------------------------------------------------------------
DIMENSION PERTST(12,2,3),EL(13),TQ(4)
DATA PERTST / 1.,1.,2.,1.,.3158,.07407,.01391,.002182,
* .0002945,.00003492,.000003692,.0000003524,
* 1.,1.,.5,.1667,.04167,1.,1.,1.,1.,1.,1.,1.,
* 2.,12.,24.,37.89,53.33,70.08,87.97,106.9,
* 126.7,147.4,168.8,191.0,
* 2.0,4.5,7.333,10.42,13.7,1.,1.,1.,1.,1.,1.,1.,
* 12.0,24.0,37.89,53.33,70.08,87.97,106.9,
* 126.7,147.4,168.8,191.0,1.,
* 3.0,6.0,9.167,12.5,1.,1.,1.,1.,1.,1.,1.,1. /
C
GO TO (1,2),METH
1 GO TO (101,102,103,104,105,106,107,108,109,110,111,112),NQ
2 GO TO (201,202,203,204,205),NQ
C-----------------------------------------------------------------------
C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY.
C FOR A GIVEN ORDER NQ, THEY CAN BE CALCULATED BY USE OF THE
C GENERATING POLYNOMIAL L(T), WHOSE COEFFICIENTS ARE EL(I)..
C L(T) = EL(1) + EL(2)*T + ... + EL(NQ+1)*T**NQ.
C FOR THE IMPLICIT ADAMS METHODS, L(T) IS GIVEN BY
C DL/DT = (T+1)*(T+2)* ... *(T+NQ-1)/K, L(-1) = 0,
C WHERE K = FACTORIAL(NQ-1).
C FOR THE BDF METHODS,
C L(T) = (T+1)*(T+2)* ... *(T+NQ)/K,
C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
C
C THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS..
C IMPLICIT ADAMS METHODS OF ORDERS 1 TO 12,
C BDF METHODS OF ORDERS 1 TO 5.
C-----------------------------------------------------------------------
101 EL(1) = 1.0D-00
GO TO 900
102 EL(1) = 0.5D-00
EL(3) = 0.5D-00
GO TO 900
103 EL(1) = 4.1666666666667D-01
EL(3) = 0.75D-00
EL(4) = 1.6666666666667D-01
GO TO 900
104 EL(1) = 0.375D-00
EL(3) = 9.1666666666667D-01
EL(4) = 3.3333333333333D-01
EL(5) = 4.1666666666667D-02
GO TO 900
105 EL(1) = 3.4861111111111D-01
EL(3) = 1.0416666666667D-00
EL(4) = 4.8611111111111D-01
EL(5) = 1.0416666666667D-01
EL(6) = 8.3333333333333D-03
GO TO 900
106 EL(1) = 3.2986111111111D-01
EL(3) = 1.1416666666667D-00
EL(4) = 0.625D-00
EL(5) = 1.7708333333333D-01
EL(6) = 0.025D-00
EL(7) = 1.3888888888889D-03
GO TO 900
107 EL(1) = 3.1559193121693D-01
EL(3) = 1.225D-00
EL(4) = 7.5185185185185D-01
EL(5) = 2.5520833333333D-01
EL(6) = 4.8611111111111D-02
EL(7) = 4.8611111111111D-03
EL(8) = 1.9841269841270D-04
GO TO 900
108 EL(1) = 3.0422453703704D-01
EL(3) = 1.2964285714286D-00
EL(4) = 8.6851851851852D-01
EL(5) = 3.3576388888889D-01
EL(6) = 7.7777777777778D-02
EL(7) = 1.0648148148148D-02
EL(8) = 7.9365079365079D-04
EL(9) = 2.4801587301587D-05
GO TO 900
109 EL(1) = 2.9486800044092D-01
EL(3) = 1.3589285714286D-00
EL(4) = 9.7655423280423D-01
EL(5) = 0.4171875D-00
EL(6) = 1.1135416666667D-01
EL(7) = 0.01875D-00
EL(8) = 1.9345238095238D-03
EL(9) = 1.1160714285714D-04
EL(10)= 2.7557319223986D-06
GO TO 900
110 EL(1) = 2.8697544642857D-01
EL(3) = 1.4144841269841D-00
EL(4) = 1.0772156084656D-00
EL(5) = 4.9856701940035D-01
EL(6) = 0.1484375D-00
EL(7) = 2.9060570987654D-02
EL(8) = 3.7202380952381D-03
EL(9) = 2.9968584656085D-04
EL(10)= 1.3778659611993D-05
EL(11)= 2.7557319223986D-07
GO TO 900
111 EL(1) = 2.8018959644394D-01
EL(3) = 1.4644841269841D-00
EL(4) = 1.1715145502646D-00
EL(5) = 5.7935819003527D-01
EL(6) = 1.8832286155203D-01
EL(7) = 4.1430362654321D-02
EL(8) = 6.2111441798942D-03
EL(9) = 6.2520667989418D-04
EL(10)= 4.0417401528513D-05
EL(11)= 1.5156525573192D-06
EL(12)= 2.5052108385442D-08
GO TO 900
112 EL(1) = 2.7426554003160D-01
EL(3) = 1.5099386724387D-00
EL(4) = 1.2602711640212D-00
EL(5) = 6.5923418209877D-01
EL(6) = 2.3045800264550D-01
EL(7) = 5.5697246105232D-02
EL(8) = 9.4394841269841D-03
EL(9) = 1.1192749669312D-03
EL(10)= 9.0939153439153D-05
EL(11)= 4.8225308641975D-06
EL(12)= 1.5031265031265D-07
EL(13)= 2.0876756987868D-09
GO TO 900
201 EL(1) = 1.0D-00
GO TO 900
202 EL(1) = 6.6666666666667D-01
EL(3) = 3.3333333333333D-01
GO TO 900
203 EL(1) = 5.4545454545455D-01
EL(3) = EL(1)
EL(4) = 9.0909090909091D-02
GO TO 900
204 EL(1) = 0.48D-00
EL(3) = 0.7D-00
EL(4) = 0.2D-00
EL(5) = 0.02D-00
GO TO 900
205 EL(1) = 4.3795620437956D-01
EL(3) = 8.2116788321168D-01
EL(4) = 3.1021897810219D-01
EL(5) = 5.4744525547445D-02
EL(6) = 3.6496350364964D-03
C
900 DO 910 K = 1,3
910 TQ(K) = PERTST(NQ,METH,K)
TQ(4) = .5D-00*TQ(2)/ FLOAT(NQ+2)
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE DECB ( NDIM, N, ML, MU, B, IPIV, IER )
IMPLICIT REAL*8 (A-H, O-Z)
C-----------------------------------------------------------------------
C SUBROUTINES DECB AND SOLB FORM A TWO SUBROUTINE PACKAGE FOR THE
C DIRECT SOLUTION OF A SYSTEM OF LINEAR EQUATIONS IN WHICH THE
C COEFFICIENT MATRIX IS REAL AND BANDED.
C
C LU DECOMPOSITION OF BAND MATRIX A.. L*U = P*A , WHERE P IS A
C PERMUTATION MATRIX, L IS A UNIT LOWER TRIANGULAR MATRIX,
C AND U IS AN UPPER TRIANGULAR MATRIX.
C N = ORDER OF MATRIX.
C B = N BY (2*ML+MU+1) ARRAY CONTAINING THE MATRIX A ON INPUT
C AND ITS FACTORED FORM ON OUTPUT.
C ON INPUT, B(I,K) (1.LE.I.LE.N) CONTAINS THE K-TH
C DIAGONAL OF A, OR A(I,J) IS STORED IN B(I,J-I+ML+1).
C ON OUTPUT, B CONTAINS THE L AND U FACTORS, WITH
C U IN COLUMNS 1 TO ML+MU+1, AND L IN COLUMNS
C ML+MU+2 TO 2*ML+MU+1.
C ML,MU = WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, NOT
C COUNTING THE MAIN DIAGONAL. TOTAL BANDWIDTH IS ML+MU+1.
C NDIM = THE FIRST DIMENSION (COLUMN LENGTH) OF THE ARRAY B.
C NDIM MUST BE .GE. N.
C IPIV = ARRAY OF LENGTH N CONTAINING PIVOT INFORMATION.
C IER = ERROR INDICATOR..
C = 0 IF NO ERROR,
C = K IF THE K-TH PIVOT CHOSEN WAS ZERO (A IS SINGULAR).
C THE INPUT ARGUMENTS ARE NDIM, N, ML, MU, B.
C THE OUTPUT ARGUMENTS ARE B, IPIV, IER.
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. NONE
C CALLED BY.. DIFFUN,INITAL,PSETIB
C FORTRAN FUNCTIONS USED.. ABS,MIN0
C-----------------------------------------------------------------------
DIMENSION B(NDIM,2*ML+MU+1),IPIV(N)
IER = 0
IF (N .EQ. 1) GO TO 92
LL = ML + MU + 1
N1 = N - 1
IF (ML .EQ. 0) GO TO 32
DO 30 I = 1,ML
II = MU + I
K = ML + 1 - I
DO 10 J = 1,II
10 B(I,J) = B(I,J+K)
K = II + 1
DO 20 J = K,LL
20 B(I,J) = 0.
30 CONTINUE
32 LR = ML
DO 90 NR = 1,N1
NP = NR + 1
IF (LR .NE. N) LR = LR + 1
MX = NR
XM = DABS(B(NR,1))
IF (ML .EQ. 0) GO TO 42
DO 40 I = NP,LR
IF ( DABS(B(I,1)) .LE. XM) GO TO 40
MX = I
XM = DABS(B(I,1))
40 CONTINUE
42 IPIV(NR) = MX
IF (MX .EQ. NR) GO TO 60
DO 50 I = 1,LL
XX = B(NR,I)
B(NR,I) = B(MX,I)
50 B(MX,I) = XX
60 XM = B(NR,1)
IF (XM .EQ. 0.) GO TO 100
B(NR,1) = 1./XM
IF (ML .EQ. 0) GO TO 90
XM = -B(NR,1)
KK = MIN0(N-NR,LL-1)
DO 80 I = NP,LR
J = LL + I - NR
XX = B(I,1)*XM
B(NR,J) = XX
DO 70 II = 1,KK
70 B(I,II) = B(I,II+1) + XX*B(NR,II+1)
80 B(I,LL) = 0.
90 CONTINUE
92 NR = N
IF (B(N,1) .EQ. 0.) GO TO 100
B(N,1) = 1./B(N,1)
RETURN
100 IER = NR
RETURN
END
C
C
C ##############################################################################
C
C
SUBROUTINE SOLB ( NDIM, N, ML, MU, B, Y, IPIV )
IMPLICIT REAL*8 (A-H, O-Z)
C-----------------------------------------------------------------------
C SUBROUTINES DECB AND SOLB FORM A TWO SUBROUTINE PACKAGE FOR THE
C DIRECT SOLUTION OF A SYSTEM OF LINEAR EQUATIONS IN WHICH THE
C COEFFICIENT MATRIX IS REAL AND BANDED.
C
C SOLUTION OF A*X = C GIVEN LU DECOMPOSITION OF A FROM DECB.
C Y = RIGHT-HAND VECTOR C, OF LENGTH N, ON INPUT,
C = SOLUTION VECTOR X ON OUTPUT.
C ALL THE ARGUMENTS ARE INPUT ARGUMENTS.
C THE OUTPUT ARGUMENT IS Y.
C
C PACKAGE ROUTINES CALLED.. NONE
C USER ROUTINES CALLED.. NONE
C CALLED BY.. DIFFUN,INITAL,STIFIB
C FORTRAN FUNCTIONS USED.. MIN0
C-----------------------------------------------------------------------
DIMENSION B(NDIM,2*ML+MU+1),Y(N),IPIV(N)
IF (N .EQ. 1) GO TO 60
N1 = N - 1
LL = ML + MU + 1
IF (ML .EQ. 0) GO TO 32
DO 30 NR = 1,N1
IF (IPIV(NR) .EQ. NR) GO TO 10
J = IPIV(NR)
XX = Y(NR)
Y(NR) = Y(J)
Y(J) = XX
10 KK = MIN0(N-NR,ML)
DO 20 I = 1,KK
20 Y(NR+I) = Y(NR+I) + Y(NR)*B(NR,LL+I)
30 CONTINUE
32 LL = LL - 1
Y(N) = Y(N)*B(N,1)
KK = 0
DO 50 NB = 1,N1
NR = N - NB
IF (KK .NE. LL) KK = KK + 1
DP = 0.
IF (LL .EQ. 0) GO TO 50
DO 40 I = 1,KK
40 DP = DP + B(NR,I+1)*Y(NR+I)
50 Y(NR) = (Y(NR) - DP)*B(NR,1)
RETURN
60 Y(1) = Y(1)*B(1,1)
RETURN
END
C ------------------------------------------------------------------------------
Computing file changes ...
