SUBROUTINE RADAU5(N,FCN,X,Y,XEND,H,
& RTOL,ATOL,ITOL,
& JAC ,IJAC,MLJAC,MUJAC,
& MAS ,IMAS,MLMAS,MUMAS,
& SOLOUT,IOUT,
& WORK,LWORK,IWORK,LIWORK,LRCONT,IDID)
C ----------------------------------------------------------
C NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
C SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS
C M*Y'=F(X,Y).
C THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M .NE. I)
C OR EXPLICIT (M=I).
C THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA)
C OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT.
C C.F. SECTION IV.8
C
C AUTHORS: E. HAIRER AND G. WANNER
C UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
C CH-1211 GENEVE 24, SWITZERLAND
C E-MAIL: HAIRER@CGEUGE51.BITNET, WANNER@CGEUGE51.BITNET
C
C THIS CODE IS PART OF THE BOOK:
C E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
C EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
C SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
C SPRINGER-VERLAG (1990)
C
C VERSION OF NOVEMBER 14, 1989
C
C INPUT PARAMETERS
C ----------------
C N DIMENSION OF THE SYSTEM
C
C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
C VALUE OF F(X,Y):
C SUBROUTINE FCN(N,X,Y,F)
C REAL*8 X,Y(N),F(N)
C F(1)=... ETC.
C
C X INITIAL X-VALUE
C
C Y(N) INITIAL VALUES FOR Y
C
C XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
C
C H INITIAL STEP SIZE GUESS;
C FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
C H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD.
C THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
C ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW
C STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE.
C (IF H=0.D0, THE CODE PUTS H=1.D-6).
C
C RTOL,ATOL RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
C CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
C
C ITOL SWITCH FOR RTOL AND ATOL:
C ITOL=0: BOTH RTOL AND ATOL ARE SCALARS.
C THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
C Y(I) BELOW RTOL*ABS(Y(I))+ATOL
C ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
C THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
C RTOL(I)*ABS(Y(I))+ATOL(I).
C
C JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
C THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
C (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY
C A DUMMY SUBROUTINE IN THE CASE IJAC=0).
C FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM
C SUBROUTINE JAC(N,X,Y,DFY,LDFY)
C REAL*8 X,Y(N),DFY(LDFY,N)
C DFY(1,1)= ...
C LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS
C FURNISHED BY THE CALLING PROGRAM.
C IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO
C BE FULL AND THE PARTIAL DERIVATIVES ARE
C STORED IN DFY AS
C DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
C ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
C THE PARTIAL DERIVATIVES ARE STORED
C DIAGONAL-WISE AS
C DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
C
C IJAC SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
C IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
C DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
C IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
C
C MLJAC SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
C MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
C 0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN
C MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
C THE MAIN DIAGONAL).
C
C MUJAC UPPER BANDWITH OF JACOBIAN MATRIX (>= NUMBER OF NON-
C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
C NEED NOT BE DEFINED IF MLJAC=N.
C
C ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS -----
C ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
C
C MAS NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
C MATRIX M.
C IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
C MATRIX AND NEEDS NOT TO BE DEFINED;
C SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
C IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
C SUBROUTINE MAS(N,AM,LMAS)
C REAL*8 AM(LMAS,N)
C AM(1,1)= ....
C IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
C AS FULL MATRIX LIKE
C AM(I,J) = M(I,J)
C ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
C DIAGONAL-WISE AS
C AM(I-J+MUMAS+1,J) = M(I,J).
C
C IMAS GIVES INFORMATION ON THE MASS-MATRIX:
C IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
C MATRIX, MAS IS NEVER CALLED.
C IMAS=1: MASS-MATRIX IS SUPPLIED.
C
C MLMAS SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
C MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
C 0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
C MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
C THE MAIN DIAGONAL).
C MLMAS IS SUPPOSED TO BE .LE. MLJAC.
C
C MUMAS UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
C NEED NOT BE DEFINED IF MLMAS=N.
C MUMAS IS SUPPOSED TO BE .LE. MUJAC.
C
C SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
C NUMERICAL SOLUTION DURING INTEGRATION.
C IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
C SUPPLY A DUMMY SUBROUTINE IF IOUT=0.
C IT MUST HAVE THE FORM
C SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN)
C REAL*8 X,Y(N)
C ....
C SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
C GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
C THE FIRST GRID-POINT).
C "XOLD" IS THE PRECEEDING GRID-POINT.
C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
C IS SET <0, RADAU5 RETURNS TO THE CALLING PROGRAM.
C
C ----- CONTINUOUS OUTPUT: -----
C DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
C FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
C THE REAL*8 FUNCTION
C >>> CONTR5(I,S) <<<
C WHICH PROVIDES AN APPROXIMATION TO THE I-TH
C COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
C S SHOULD LIE IN THE INTERVAL [XOLD,X].
C
C IOUT SWITCH FOR CALLING THE SUBROUTINE SOLOUT:
C IOUT=0: SUBROUTINE IS NEVER CALLED
C IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT.
C
C WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK".
C WORK(1), WORK(2),.., WORK(7) SERVE AS PARAMETERS
C FOR THE CODE. FOR STANDARD USE OF THE CODE
C WORK(1),..,WORK(7) MUST BE SET TO ZERO BEFORE
C CALLING. SEE BELOW FOR A MORE SOPHISTICATED USE.
C WORK(8),..,WORK(LWORK) SERVE AS WORKING SPACE
C FOR ALL VECTORS AND MATRICES.
C "LWORK" MUST BE AT LEAST
C N*(LJAC+LMAS+3*LE+8)+7
C WHERE
C LJAC=N IF MLJAC=N (FULL JACOBIAN)
C LJAC=MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.)
C AND
C LMAS=0 IF IMAS=0
C LMAS=N IF IMAS=1 AND MLMAS=N (FULL)
C LMAS=MLMAS+MUMAS+1 IF MLMAS<N (BANDED MASS-M.)
C AND
C LE=N IF MLJAC=N (FULL JACOBIAN)
C LE=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.)
C
C IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
C MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
C STORAGE REQUIREMENT IS
C LWORK = 4*N*N+8*N+7.
C
C LWORK DECLARED LENGHT OF ARRAY "WORK".
C
C IWORK INTEGER WORKING SPACE OF LENGHT "LIWORK".
C IWORK(1),IWORK(2),...,IWORK(7) SERVE AS PARAMETERS
C FOR THE CODE. FOR STANDARD USE, SET IWORK(1),..,
C IWORK(7) TO ZERO BEFORE CALLING.
C IWORK(8),...,IWORK(LIWORK) SERVE AS WORKING AREA.
C "LIWORK" MUST BE AT LEAST 3*N+7.
C
C LIWORK DECLARED LENGHT OF ARRAY "IWORK".
C
C LRCONT DECLARED LENGTH OF COMMON BLOCK
C >>> COMMON /CONT/ICONT(3),RCONT(LRCONT) <<<
C WHICH MUST BE DECLARED IN THE CALLING PROGRAM.
C "LRCONT" MUST BE AT LEAST
C 4*N+4 .
C THIS IS USED FOR STORING THE COEFFICIENTS OF THE
C CONTINUOUS SOLUTION AND MAKES THE CALLING LIST FOR THE
C FUNCTION "CONTR5" AS SIMPLE AS POSSIBLE.
C
C ----------------------------------------------------------------------
C
C SOPHISTICATED SETTING OF PARAMETERS
C -----------------------------------
C SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK
C WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(7)
C AS WELL AS IWORK(1),..,IWORK(7) DIFFERENT FROM ZERO.
C FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
C
C IWORK(1) IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN
C MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY
C ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN.
C IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N)
C AND NOT FOR IMPLICIT SYSTEMS (IMAS=1). IT IS
C ALSO NOT GOOD FOR SPARSE JACOBIANS.
C
C IWORK(2) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
C THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000.
C
C IWORK(3) THE MAXIMUM NUMBER OF NEWTON ITERATIONS FOR THE
C SOLUTION OF THE IMPLICIT SYSTEM IN EACH STEP.
C THE DEFAULT VALUE (FOR IWORK(3)=0) IS 7.
C
C IWORK(4) IF IWORK(4).EQ.0 THE EXTRAPOLATED COLLOCATION SOLUTION
C IS TAKEN AS STARTING VALUE FOR NEWTON'S METHOD.
C IF IWORK(4).NE.0 ZERO STARTING VALUES ARE USED.
C THE LATTER IS RECOMMENDED IF NEWTON'S METHOD HAS
C DIFFICULTIES WITH CONVERGENCE (THIS IS THE CASE WHEN
C NSTEP IS LARGER THAN NACCPT + NREJCT).
C DEFAULT IS IWORK(4)=0.
C
C THE FOLLOWING 3 PARAMETERS ARE IMPORTANT FOR
C DIFFERENTIAL-ALGEBRAIC SYSTEMS OF INDEX > 1.
C THE FUNCTION-SUBROUTINE SHOULD BE WRITTEN SUCH THAT
C THE INDEX 1,2,3 VARIABLES APPEAR IN THIS ORDER.
C IN ESTIMATING THE ERROR THE INDEX 2 VARIABLES ARE
C MULTIPLIED BY H, THE INDEX 3 VARIABLES BY H**2.
C
C IWORK(5) DIMENSION OF THE INDEX 1 VARIABLES (MUST BE > 0). FOR
C ODE'S THIS EQUALS THE DIMENSION OF THE SYSTEM.
C DEFAULT IWORK(5)=N.
C
C IWORK(6) DIMENSION OF THE INDEX 2 VARIABLES. DEFAULT IWORK(6)=0.
C
C IWORK(7) DIMENSION OF THE INDEX 3 VARIABLES. DEFAULT IWORK(7)=0.
C
C
C WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
C
C WORK(2) THE SAFETY FACTOR IN STEP SIZE PREDICTION,
C DEFAULT 0.9D0.
C
C WORK(3) DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
C INCREASE WORK(3), TO 0.1 SAY, WHEN JACOBIAN EVALUATIONS
C ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER
C (0.001D0, SAY). NEGATIV WORK(3) FORCES THE CODE THE
C COMPUTE THE JACOBIAN AFTER EVERY ACCEPTED STEP.
C DEFAULT 0.001D0.
C
C WORK(4) STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
C SMALLER VALUES OF WORK(4) MAKE THE CODE SLOWER, BUT SAFER.
C DEFAULT 0.03D0.
C
C WORK(5) AND WORK(6) : IF WORK(5) < HNEW/HOLD < WORK(6), THEN THE
C STEP SIZE IS NOT CHANGED. THIS SAVES, TOGETHER WITH A
C LARGE WORK(3), LU-DECOMPOSITIONS AND COMPUTING TIME FOR
C LARGE SYSTEMS. FOR SMALL SYSTEMS ONE MAY HAVE
C WORK(5)=1.D0, WORK(6)=1.2D0, FOR LARGE FULL SYSTEMS
C WORK(5)=0.99D0, WORK(6)=2.D0 MIGHT BE GOOD.
C DEFAULTS WORK(5)=1.D0, WORK(6)=1.2D0 .
C
C WORK(7) MAXIMAL STEP SIZE, DEFAULT XEND-X.
C
C-----------------------------------------------------------------------
C
C OUTPUT PARAMETERS
C -----------------
C X X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
C (AFTER SUCCESSFUL RETURN X=XEND).
C
C Y(N) NUMERICAL SOLUTION AT X
C
C H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
C
C IDID REPORTS ON SUCCESSFULNESS UPON RETURN:
C IDID=1 COMPUTATION SUCCESSFUL,
C IDID=-1 COMPUTATION UNSUCCESSFUL.
C
C-----------------------------------------------------------------------
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C DECLARATIONS
C *** *** *** *** *** *** *** *** *** *** *** *** ***
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION Y(N),ATOL(1),RTOL(1),WORK(LWORK),IWORK(LIWORK)
LOGICAL IMPLCT,JBAND,ARRET,STARTN
EXTERNAL FCN,JAC,MAS,SOLOUT
COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
C --- COMMON STAT CAN BE INSPECTED FOR STATISTICAL PURPOSES:
C --- NFCN NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
C EVALUATION OF THE JACOBIAN ARE NOT COUNTED)
C --- NJAC NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
C OR NUMERICALLY)
C --- NSTEP NUMBER OF COMPUTED STEPS
C --- NACCPT NUMBER OF ACCEPTED STEPS
C --- NREJCT NUMBER OF REJECTED STEPS (DUE TO ERROR TEST),
C (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED)
C --- NDEC NUMBER OF LU-DECOMPOSITIONS OF BOTH MATRICES
C --- NSOL NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS, OF BOTH
C SYSTEMS; THE NSTEP FORWARD-BACKWARD SUBSTITUTIONS,
C NEEDED FOR STEP SIZE SELECTION, ARE NOT COUNTED
C *** *** *** *** *** *** ***
C SETTING THE PARAMETERS
C *** *** *** *** *** *** ***
NFCN=0
NJAC=0
NSTEP=0
NACCPT=0
NREJCT=0
NDEC=0
NSOL=0
ARRET=.FALSE.
C -------- SWITCH FOR TRANSFORMATION OF JACOBIAN TO HESSIAN FORM ---
NHESS=IWORK(1)
IF (N.LE.2) NHESS=0
C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
IF(IWORK(2).EQ.0)THEN
NMAX=100000
ELSE
NMAX=IWORK(2)
IF(NMAX.LE.0)THEN
WRITE(4,*)' WRONG INPUT IWORK(2)=',IWORK(2)
ARRET=.TRUE.
END IF
END IF
C -------- NIT MAXIMAL NUMBER OF NEWTON ITERATIONS
IF(IWORK(3).EQ.0)THEN
NIT=7
ELSE
NIT=IWORK(3)
IF(NIT.LE.0)THEN
WRITE(4,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
ARRET=.TRUE.
END IF
END IF
C -------- STARTN SWITCH FOR STARTING VALUES OF NEWTON ITERATIONS
IF(IWORK(4).EQ.0)THEN
STARTN=.FALSE.
ELSE
STARTN=.TRUE.
END IF
C -------- PARAMETER FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS
NIND1=IWORK(5)
NIND2=IWORK(6)
NIND3=IWORK(7)
IF (NIND1.EQ.0) NIND1=N
IF (NIND1+NIND2+NIND3.NE.N)THEN
WRITE(4,*)' CURIOUS INPUT FOR WORK(5,6,7)=',NIND1,NIND2,NIND3
ARRET=.TRUE.
END IF
C -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0
IF(WORK(1).EQ.0.D0)THEN
UROUND=1.D-16
ELSE
UROUND=WORK(1)
IF(UROUND.LE.1.D-19.OR.UROUND.GE.1.D0)THEN
WRITE(4,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1)
ARRET=.TRUE.
END IF
END IF
C --------- SAFE SAFETY FACTOR IN STEP SIZE PREDICTION
IF(WORK(2).EQ.0.D0)THEN
SAFE=0.9D0
ELSE
SAFE=WORK(2)
IF(SAFE.LE..001D0.OR.SAFE.GE.1.D0)THEN
WRITE(4,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
ARRET=.TRUE.
END IF
END IF
C ------ THET DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
IF(WORK(3).EQ.0.D0)THEN
THET=0.001D0
ELSE
THET=WORK(3)
END IF
C --- FNEWT STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
IF(WORK(4).EQ.0.D0)THEN
FNEWT=0.03D0
ELSE
FNEWT=WORK(4)
END IF
C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST.
IF(WORK(5).EQ.0.D0)THEN
QUOT1=1.D0
ELSE
QUOT1=WORK(5)
END IF
IF(WORK(6).EQ.0.D0)THEN
QUOT2=1.2D0
ELSE
QUOT2=WORK(6)
END IF
C -------- MAXIMAL STEP SIZE
IF(WORK(7).EQ.0.D0)THEN
HMAX=XEND-X
ELSE
HMAX=WORK(7)
END IF
C --------- CHECK IF TOLERANCES ARE O.K.
IF (ITOL.EQ.0) THEN
IF (ATOL(1).LE.0.D0.OR.RTOL(1).LE.10.D0*UROUND) THEN
WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
ARRET=.TRUE.
END IF
ELSE
DO 15 I=1,N
IF (ATOL(I).LE.0.D0.OR.RTOL(I).LE.10.D0*UROUND) THEN
WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
ARRET=.TRUE.
END IF
15 CONTINUE
END IF
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C COMPUTATION OF ARRAY ENTRIES
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C ---- IMPLICIT, BANDED OR NOT ?
IMPLCT=IMAS.NE.0
JBAND=MLJAC.NE.N
C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
C -- JACOBIAN AND MATRICES E1, E2
IF(JBAND)THEN
LDJAC=MLJAC+MUJAC+1
LDE1=2*MLJAC+MUJAC+1
ELSE
LDJAC=N
LDE1=N
END IF
C -- MASS MATRIX
IF (IMPLCT) THEN
IF (MLMAS.NE.N) THEN
LDMAS=MLMAS+MUMAS+1
ELSE
LDMAS=N
END IF
C ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC"
IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
WRITE (6,*) 'BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF
& "JAC"'
ARRET=.TRUE.
END IF
ELSE
LDMAS=0
END IF
LDMAS2=MAX(1,LDMAS)
C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN
IF ((IMPLCT.OR.JBAND).AND.NHESS.NE.0) THEN
WRITE(4,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH
&FULL JACOBIAN'
ARRET=.TRUE.
END IF
C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
IEZ1=8
IEZ2=IEZ1+N
IEZ3=IEZ2+N
IEY0=IEZ3+N
IESCAL=IEY0+N
IEF1=IESCAL+N
IEF2=IEF1+N
IEF3=IEF2+N
IEJAC=IEF3+N
IEMAS=IEJAC+N*LDJAC
IEE1=IEMAS+N*LDMAS
IEE2R=IEE1+N*LDE1
IEE2I=IEE2R+N*LDE1
C ------ TOTAL STORAGE REQUIREMENT -----------
ISTORE=IEE2I+N*LDE1-1
IF(ISTORE.GT.LWORK)THEN
WRITE(4,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
ARRET=.TRUE.
END IF
C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
IEIP1=8
IEIP2=IEIP1+N
IEIPH=IEIP2+N
C --------- TOTAL REQUIREMENT ---------------
ISTORE=IEIPH+N-1
IF(ISTORE.GT.LIWORK)THEN
WRITE(4,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
ARRET=.TRUE.
END IF
C --------- CONTROL OF LENGTH OF COMMON BLOCK "CONT" -------
IF(LRCONT.LT.(4*N+4))THEN
WRITE(4,*)' INSUFF. STORAGE FOR RCONT, MIN. LRCONT=',4*N+4
ARRET=.TRUE.
END IF
C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
IF (ARRET) THEN
IDID=-1
RETURN
END IF
C -------- CALL TO CORE INTEGRATOR ------------
CALL RADCOR(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,
& JAC,IJAC,MLJAC,MUJAC,MAS,IMAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
& NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,NHESS,STARTN,
& NIND1,NIND2,NIND3,
& IMPLCT,JBAND,LDJAC,LDE1,LDMAS2,WORK(IEZ1),WORK(IEZ2),
& WORK(IEZ3),WORK(IEY0),WORK(IESCAL),WORK(IEF1),WORK(IEF2),
& WORK(IEF3),WORK(IEJAC),WORK(IEE1),WORK(IEE2R),WORK(IEE2I),
& WORK(IEMAS),IWORK(IEIP1),IWORK(IEIP2),IWORK(IEIPH))
C ----------- RETURN -----------
RETURN
END
C
C
C
C ----- ... AND HERE IS THE CORE INTEGRATOR ----------
C
SUBROUTINE RADCOR(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,
& JAC,IJAC,MLJAC,MUJAC,MAS,IMAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
& NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,NHESS,STARTN,
& NIND1,NIND2,NIND3,IMPLCT,BANDED,LDJAC,LDE1,LDMAS,Z1,Z2,Z3,
& Y0,SCAL,F1,F2,F3,FJAC,E1,E2R,E2I,FMAS,IP1,IP2,IPHES)
C ----------------------------------------------------------
C CORE INTEGRATOR FOR RADAU5
C PARAMETERS SAME AS IN RADAU5 WITH WORKSPACE ADDED
C ----------------------------------------------------------
C DECLARATIONS
C ----------------------------------------------------------
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 Y(N),Z1(N),Z2(N),Z3(N),Y0(N),SCAL(N),F1(N),F2(N),F3(N)
REAL*8 FJAC(LDJAC,N),FMAS(LDMAS,N)
REAL*8 E1(LDE1,N),E2R(LDE1,N),E2I(LDE1,N)
REAL*8 ATOL(1),RTOL(1)
INTEGER IP1(N),IP2(N),IPHES(N)
COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
COMMON /CONT/NN,NN2,NN3,XSOL,HSOL,C2M1,C1M1,CONT(1)
LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,STARTN
LOGICAL INDEX1,INDEX2,INDEX3
C *** *** *** *** *** *** ***
C INITIALISATIONS
C *** *** *** *** *** *** ***
C --------- DUPLIFY N FOR COMMON BLOCK CONT -----
NN=N
NN2=2*N
NN3=3*N
C -------- CHECK THE INDEX OF THE PROBLEM -----
INDEX1=NIND1.NE.0
INDEX2=NIND2.NE.0
INDEX3=NIND3.NE.0
C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
IF(IMPLCT)CALL MAS(N,FMAS,LDMAS)
C ---------- CONSTANTS ---------
SQ6=DSQRT(6.D0)
C1=(4.D0-SQ6)/10.D0
C2=(4.D0+SQ6)/10.D0
C1M1=C1-1.D0
C2M1=C2-1.D0
C1MC2=C1-C2
DD1=-(13.D0+7.D0*SQ6)/3.D0
DD2=(-13.D0+7.D0*SQ6)/3.D0
DD3=-1.D0/3.D0
U1=(6.D0+81.D0**(1.D0/3.D0)-9.D0**(1.D0/3.D0))/30.D0
ALPH=(12.D0-81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))/60.D0
BETA=(81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))*DSQRT(3.D0)/60.D0
CNO=ALPH**2+BETA**2
T11=9.1232394870892942792D-02
T12=-0.14125529502095420843D0
T13=-3.0029194105147424492D-02
T21=0.24171793270710701896D0
T22=0.20412935229379993199D0
T23=0.38294211275726193779D0
T31=0.96604818261509293619D0
TI11=4.3255798900631553510D0
TI12=0.33919925181580986954D0
TI13=0.54177053993587487119D0
TI21=-4.1787185915519047273D0
TI22=-0.32768282076106238708D0
TI23=0.47662355450055045196D0
TI31=-0.50287263494578687595D0
TI32=2.5719269498556054292D0
TI33=-0.59603920482822492497D0
POSNEG=SIGN(1.D0,XEND-X)
HMAXN=MIN(ABS(HMAX),ABS(XEND-X))
IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
H=MIN(ABS(H),HMAXN)
H=SIGN(H,POSNEG)
HOLD=H
REJECT=.FALSE.
FIRST=.TRUE.
FACCON=1.D0
CFAC=SAFE*(1+2*NIT)
NSING=0
XOLD=X
IF (IOUT.NE.0) THEN
IRTRN=1
NRSOL=1
XOSOL=XOLD
XSOL=X
DO 7 I=1,N
7 CONT(I)=Y(I)
NSOLU=N
HSOL=HOLD
CALL SOLOUT(NRSOL,XOSOL,XSOL,CONT,NSOLU,IRTRN)
IF (IRTRN.LT.0) GOTO 79
END IF
MLE=MLJAC
MUE=MUJAC
MBJAC=MLJAC+MUJAC+1
MBB=MLMAS+MUMAS+1
MDIAG=MLE+MUE+1
MDIFF=MLE+MUE-MUMAS
MBDIAG=MUMAS+1
N2=2*N
N3=3*N
IF (ITOL.EQ.0) THEN
DO 8 I=1,N
8 SCAL(I)=ATOL(1)+RTOL(1)*DABS(Y(I))
ELSE
DO 9 I=1,N
9 SCAL(I)=ATOL(I)+RTOL(I)*DABS(Y(I))
END IF
HHFAC=H
CALL FCN(N,X,Y,Y0)
NFCN=NFCN+1
C --- BASIC INTEGRATION STEP
10 CONTINUE
C *** *** *** *** *** *** ***
C COMPUTATION OF THE JACOBIAN
C *** *** *** *** *** *** ***
NJAC=NJAC+1
IF (IJAC.EQ.0) THEN
C --- COMPUTE JACOBIAN MATRIX NUMERICALLY
IF (BANDED) THEN
C --- JACOBIAN IS BANDED
MUJACP=MUJAC+1
MD=MIN(MBJAC,N)
DO 16 K=1,MD
J=K
12 F1(J)=Y(J)
F2(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
Y(J)=Y(J)+F2(J)
J=J+MD
IF (J.LE.N) GOTO 12
CALL FCN(N,X,Y,CONT)
J=K
LBEG=MAX(1,J-MUJAC)
14 LEND=MIN(N,J+MLJAC)
Y(J)=F1(J)
MUJACJ=MUJACP-J
DO 15 L=LBEG,LEND
15 FJAC(L+MUJACJ,J)=(CONT(L)-Y0(L))/F2(J)
J=J+MD
LBEG=LEND+1
IF (J.LE.N) GOTO 14
16 CONTINUE
ELSE
C --- JACOBIAN IS FULL
DO 18 I=1,N
YSAFE=Y(I)
DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
Y(I)=YSAFE+DELT
CALL FCN(N,X,Y,CONT)
DO 17 J=1,N
17 FJAC(J,I)=(CONT(J)-Y0(J))/DELT
18 Y(I)=YSAFE
END IF
ELSE
C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
CALL JAC(N,X,Y,FJAC,LDJAC)
END IF
CALJAC=.TRUE.
IF (NHESS.NE.0) CALL ELMHES (LDJAC,N,1,N,FJAC,IPHES)
20 CONTINUE
C *** *** *** *** *** *** ***
C COMPUTE THE MATRICES E1 AND E2 AND THEIR DECOMPOSITIONS
C *** *** *** *** *** *** ***
FAC1=1.D0/(H*U1)
ALPHN=ALPH/(H*CNO)
BETAN=BETA/(H*CNO)
IF (IMPLCT) THEN
IF (BANDED) THEN
C --- THE MATRIX E (B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX)
DO 127 J=1,N
I1J=MAX0(1,MUJAC+2-J)
I2J=MIN(MBJAC,MUJAC+1-J+N)
DO 125 I=I1J,I2J
FJ=-FJAC(I,J)
IMLE=I+MLE
E1(IMLE,J)=FJ
E2R(IMLE,J)=FJ
125 E2I(IMLE,J)=0.D0
I1B=MAX0(1,MUMAS+2-J)
I2B=MIN0(MBB,MUMAS+1-J+N)
DO 126 I=I1B,I2B
IB=I+MDIFF
BB=FMAS(I,J)
E1(IB,J)=E1(IB,J)+FAC1*BB
E2R(IB,J)=E2R(IB,J)+ALPHN*BB
126 E2I(IB,J)=BETAN*BB
127 CONTINUE
CALL DECB(N,LDE1,E1,MLE,MUE,IP1,IER)
IF (IER.NE.0) GOTO 78
CALL DECBC(N,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
IF (IER.NE.0) GOTO 78
ELSE
IF (MLMAS.NE.N) THEN
C --- THE MATRIX E (B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX)
DO 225 J=1,N
DO 225 I=1,N
FJ=-FJAC(I,J)
E1(I,J)=FJ
E2R(I,J)=FJ
225 E2I(I,J)=0.D0
DO 226 J=1,N
I1=MAX0(1,J-MUMAS)
I2=MIN0(N,J+MLMAS)
DO 226 I=I1,I2
BB=FMAS(I-J+MBDIAG,J)
E1(I,J)=E1(I,J)+FAC1*BB
E2R(I,J)=E2R(I,J)+ALPHN*BB
226 E2I(I,J)=BETAN*BB
CALL DEC(N,LDE1,E1,IP1,IER)
IF (IER.NE.0) GOTO 78
CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
IF (IER.NE.0) GOTO 78
ELSE
C --- THE MATRIX E (B IS A FULL MATRIX, JACOBIAN A FULL MATRIX)
IF (MLJAC.EQ.N) THEN
DO 324 J=1,N
DO 324 I=1,N
FJ=-FJAC(I,J)
BB=FMAS(I,J)
E1(I,J)=BB*FAC1+FJ
E2R(I,J)=BB*ALPHN+FJ
324 E2I(I,J)=BB*BETAN
CALL DEC(N,LDE1,E1,IP1,IER)
IF (IER.NE.0) GOTO 78
CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
IF (IER.NE.0) GOTO 78
ELSE
C --- THE MATRIX E (B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX)
END IF
END IF
END IF
ELSE
IF (BANDED) THEN
C --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX)
DO 427 J=1,N
I1J=MAX0(1,MUJAC+2-J)
I2J=MIN(MBJAC,MUJAC+1-J+N)
DO 425 I=I1J,I2J
FJ=-FJAC(I,J)
IMLE=I+MLE
E1(IMLE,J)=FJ
E2R(IMLE,J)=FJ
425 E2I(IMLE,J)=0.D0
E1(MDIAG,J)=E1(MDIAG,J)+FAC1
E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN
E2I(MDIAG,J)=BETAN
427 CONTINUE
CALL DECB(N,LDE1,E1,MLE,MUE,IP1,IER)
IF (IER.NE.0) GOTO 78
CALL DECBC(N,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
IF (IER.NE.0) GOTO 78
ELSE
C --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
IF (NHESS.EQ.0) THEN
DO 526 J=1,N
DO 525 I=1,N
FJ=-FJAC(I,J)
E1(I,J)=FJ
E2R(I,J)=FJ
525 E2I(I,J)=0.D0
E1(J,J)=E1(J,J)+FAC1
E2R(J,J)=E2R(J,J)+ALPHN
526 E2I(J,J)=BETAN
CALL DEC(N,LDE1,E1,IP1,IER)
IF (IER.NE.0) GOTO 78
CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
IF (IER.NE.0) GOTO 78
ELSE
DO 624 J=1,N-1
J1=J+1
FJ=-FJAC(J1,J)
E1(J1,J)=FJ
E2R(J1,J)=FJ
624 E2I(J1,J)=0.D0
DO 626 J=1,N
DO 625 I=1,J
FJ=-FJAC(I,J)
E1(I,J)=FJ
E2I(I,J)=0.D0
625 E2R(I,J)=FJ
E1(J,J)=E1(J,J)+FAC1
E2R(J,J)=E2R(J,J)+ALPHN
626 E2I(J,J)=BETAN
CALL DECH(N,LDE1,E1,1,IP1,IER)
IF (IER.NE.0) GOTO 78
CALL DECHC(N,LDE1,E2R,E2I,1,IP2,IER)
IF (IER.NE.0) GOTO 78
END IF
END IF
END IF
NDEC=NDEC+1
30 CONTINUE
NSTEP=NSTEP+1
IF (NSTEP.GT.NMAX.OR.X+1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
IF (INDEX2) THEN
DO 465 I=NIND1+1,NIND1+NIND2
465 SCAL(I)=SCAL(I)/HHFAC
END IF
IF (INDEX3) THEN
DO 466 I=NIND1+NIND2+1,NIND1+NIND2+NIND3
466 SCAL(I)=SCAL(I)/(HHFAC*HHFAC)
END IF
XPH=X+H
C *** *** *** *** *** *** ***
C STARTING VALUES FOR NEWTON ITERATION
C *** *** *** *** *** *** ***
IF (FIRST.OR.STARTN) THEN
DO 32 I=1,N
Z1(I)=0.D0
Z2(I)=0.D0
Z3(I)=0.D0
F1(I)=0.D0
F2(I)=0.D0
32 F3(I)=0.D0
ELSE
C3Q=H/HOLD
C1Q=C1*C3Q
C2Q=C2*C3Q
DO 35 I=1,N
AK1=CONT(I+N)
AK2=CONT(I+N2)
AK3=CONT(I+N3)
Z1I=C1Q*(AK1+(C1Q-C2M1)*(AK2+(C1Q-C1M1)*AK3))
Z2I=C2Q*(AK1+(C2Q-C2M1)*(AK2+(C2Q-C1M1)*AK3))
Z3I=C3Q*(AK1+(C3Q-C2M1)*(AK2+(C3Q-C1M1)*AK3))
Z1(I)=Z1I
Z2(I)=Z2I
Z3(I)=Z3I
F1(I)=TI11*Z1I+TI12*Z2I+TI13*Z3I
F2(I)=TI21*Z1I+TI22*Z2I+TI23*Z3I
35 F3(I)=TI31*Z1I+TI32*Z2I+TI33*Z3I
END IF
C *** *** *** *** *** *** ***
C LOOP FOR THE SIMPLIFIED NEWTON ITERATION
C *** *** *** *** *** *** ***
NEWT=0
FACCON=MAX(FACCON,UROUND)**0.8D0
THETA=ABS(THET)
40 CONTINUE
IF (NEWT.GE.NIT) GOTO 78
C --- COMPUTE THE RIGHT-HAND SIDE
DO 41 I=1,N
41 CONT(I)=Y(I)+Z1(I)
CALL FCN(N,X+C1*H,CONT,Z1)
DO 42 I=1,N
42 CONT(I)=Y(I)+Z2(I)
CALL FCN(N,X+C2*H,CONT,Z2)
DO 43 I=1,N
43 CONT(I)=Y(I)+Z3(I)
CALL FCN(N,XPH,CONT,Z3)
NFCN=NFCN+3
C --- SOLVE THE LINEAR SYSTEMS
IF (IMPLCT) THEN
IF (MLMAS.NE.N) THEN
DO 146 I=1,N
S1=0.0D0
S2=0.0D0
S3=0.0D0
J1B=MAX0(1,I-MLMAS)
J2B=MIN0(N,I+MUMAS)
DO 145 J=J1B,J2B
BB=FMAS(I-J+MBDIAG,J)
S1=S1-BB*F1(J)
S2=S2-BB*F2(J)
145 S3=S3-BB*F3(J)
A1=Z1(I)
A2=Z2(I)
A3=Z3(I)
Z1(I)=S1*FAC1 +TI11*A1+TI12*A2+TI13*A3
Z2(I)=S2*ALPHN-S3*BETAN+TI21*A1+TI22*A2+TI23*A3
146 Z3(I)=S3*ALPHN+S2*BETAN+TI31*A1+TI32*A2+TI33*A3
IF (BANDED) THEN
CALL SOLB(N,LDE1,E1,MLE,MUE,Z1,IP1)
CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
ELSE
CALL SOL(N,LDE1,E1,Z1,IP1)
CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
END IF
ELSE
DO 246 I=1,N
S1=0.0D0
S2=0.0D0
S3=0.0D0
DO 245 J=1,N
BB=FMAS(I,J)
S1=S1-BB*F1(J)
S2=S2-BB*F2(J)
245 S3=S3-BB*F3(J)
A1=Z1(I)
A2=Z2(I)
A3=Z3(I)
Z1(I)=S1*FAC1 +TI11*A1+TI12*A2+TI13*A3
Z2(I)=S2*ALPHN-S3*BETAN+TI21*A1+TI22*A2+TI23*A3
246 Z3(I)=S3*ALPHN+S2*BETAN+TI31*A1+TI32*A2+TI33*A3
CALL SOL(N,LDE1,E1,Z1,IP1)
CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
END IF
ELSE
DO 345 I=1,N
A1=Z1(I)
A2=Z2(I)
A3=Z3(I)
S2=-F2(I)
S3=-F3(I)
Z1(I)=-F1(I)*FAC1 +TI11*A1+TI12*A2+TI13*A3
Z2(I)=S2*ALPHN-S3*BETAN+TI21*A1+TI22*A2+TI23*A3
345 Z3(I)=S3*ALPHN+S2*BETAN+TI31*A1+TI32*A2+TI33*A3
IF (BANDED) THEN
CALL SOLB(N,LDE1,E1,MLE,MUE,Z1,IP1)
CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
ELSE
IF (NHESS.EQ.0) THEN
CALL SOL(N,LDE1,E1,Z1,IP1)
CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
ELSE
DO 140 MM=N-2,1,-1
MP=N-MM
MP1=MP-1
I=IPHES(MP)
IF (I.EQ.MP) GOTO 110
ZSAFE=Z1(MP)
Z1(MP)=Z1(I)
Z1(I)=ZSAFE
ZSAFE=Z2(MP)
Z2(MP)=Z2(I)
Z2(I)=ZSAFE
ZSAFE=Z3(MP)
Z3(MP)=Z3(I)
Z3(I)=ZSAFE
110 CONTINUE
DO 100 I=MP+1,N
E1IMP=FJAC(I,MP1)
Z1(I)=Z1(I)-E1IMP*Z1(MP)
Z2(I)=Z2(I)-E1IMP*Z2(MP)
100 Z3(I)=Z3(I)-E1IMP*Z3(MP)
140 CONTINUE
CALL SOLH(N,LDE1,E1,1,Z1,IP1)
CALL SOLHC(N,LDE1,E2R,E2I,1,Z2,Z3,IP2)
DO 240 MM=1,N-2
MP=N-MM
MP1=MP-1
DO 200 I=MP+1,N
E1IMP=FJAC(I,MP1)
Z1(I)=Z1(I)+E1IMP*Z1(MP)
Z2(I)=Z2(I)+E1IMP*Z2(MP)
200 Z3(I)=Z3(I)+E1IMP*Z3(MP)
I=IPHES(MP)
IF (I.EQ.MP) GOTO 240
ZSAFE=Z1(MP)
Z1(MP)=Z1(I)
Z1(I)=ZSAFE
ZSAFE=Z2(MP)
Z2(MP)=Z2(I)
Z2(I)=ZSAFE
ZSAFE=Z3(MP)
Z3(MP)=Z3(I)
Z3(I)=ZSAFE
240 CONTINUE
END IF
END IF
END IF
DO 51 I=1,N
F1(I)=F1(I)+Z1(I)
F2(I)=F2(I)+Z2(I)
51 F3(I)=F3(I)+Z3(I)
NSOL=NSOL+1
NEWT=NEWT+1
DYNO=0.D0
DO 57 I=1,N
DENOM=SCAL(I)
57 DYNO=DYNO+(Z1(I)/DENOM)**2+(Z2(I)/DENOM)**2
& +(Z3(I)/DENOM)**2
DYNO=DSQRT(DYNO/N3)
C --- BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
IF (NEWT.GE.2.AND.NEWT.LT.NIT) THEN
THETA=DYNO/DYNOLD
IF (THETA.LT.0.99D0) THEN
FACCON=THETA/(1.0D0-THETA)
DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)
QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH/FNEWT))
IF (QNEWT.GE.1.0D0) THEN
HHFAC=.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT))
H=HHFAC*H
REJECT=.TRUE.
IF (CALJAC) GOTO 20
GOTO 10
END IF
ELSE
GOTO 78
END IF
END IF
DYNOLD=DMAX1(DYNO,UROUND)
DO 58 I=1,N
F1I=F1(I)
F2I=F2(I)
F3I=F3(I)
Z1(I)=T11*F1I+T12*F2I+T13*F3I
Z2(I)=T21*F1I+T22*F2I+T23*F3I
58 Z3(I)=T31*F1I+ F2I
IF (FACCON*DYNO.GT.FNEWT) GOTO 40
C *** *** *** *** *** *** ***
C ERROR ESTIMATION
C *** *** *** *** *** *** ***
HEE1=DD1/H
HEE2=DD2/H
HEE3=DD3/H
IF (IMPLCT) THEN
DO 359 I=1,N
359 F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
IF (MLMAS.EQ.N) THEN
DO 361 I=1,N
SUM=0.D0
DO 360 J=1,N
360 SUM=SUM+FMAS(I,J)*F1(J)
361 F2(I)=SUM
ELSE
DO 363 I=1,N
SUM=0.D0
J1B=MAX0(1,I-MLMAS)
J2B=MIN0(N,I+MUMAS)
DO 362 J=J1B,J2B
362 SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J)
363 F2(I)=SUM
END IF
ELSE
DO 461 I=1,N
461 F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
END IF
DO 462 I=1,N
462 CONT(I)=F2(I)+Y0(I)
IF (BANDED) THEN
CALL SOLB(N,LDE1,E1,MLE,MUE,CONT,IP1)
ELSE
IF (NHESS.EQ.0) THEN
CALL SOL(N,LDE1,E1,CONT,IP1)
ELSE
DO 340 MM=N-2,1,-1
MP=N-MM
I=IPHES(MP)
IF (I.EQ.MP) GOTO 310
ZSAFE=CONT(MP)
CONT(MP)=CONT(I)
CONT(I)=ZSAFE
310 CONTINUE
DO 300 I=MP+1,N
300 CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
340 CONTINUE
CALL SOLH(N,LDE1,E1,1,CONT,IP1)
DO 440 MM=1,N-2
MP=N-MM
DO 400 I=MP+1,N
400 CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
I=IPHES(MP)
IF (I.EQ.MP) GOTO 440
ZSAFE=CONT(MP)
CONT(MP)=CONT(I)
CONT(I)=ZSAFE
440 CONTINUE
END IF
END IF
ERR=0.D0
DO 354 I=1,N
354 ERR=ERR+(CONT(I)/SCAL(I))**2
ERR=DMAX1(DSQRT(ERR/N),1.D-10)
IF ((FIRST.OR.REJECT).AND.ERR.GE.1.D0) THEN
DO 460 I=1,N
460 CONT(I)=Y(I)+CONT(I)
CALL FCN(N,X,CONT,F1)
NFCN=NFCN+1
DO 463 I=1,N
463 CONT(I)=F1(I)+F2(I)
IF (BANDED) THEN
CALL SOLB(N,LDE1,E1,MLE,MUE,CONT,IP1)
ELSE
IF (NHESS.EQ.0) THEN
CALL SOL(N,LDE1,E1,CONT,IP1)
ELSE
DO 540 MM=N-2,1,-1
MP=N-MM
I=IPHES(MP)
IF (I.EQ.MP) GOTO 510
ZSAFE=CONT(MP)
CONT(MP)=CONT(I)
CONT(I)=ZSAFE
510 CONTINUE
DO 500 I=MP+1,N
500 CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
540 CONTINUE
CALL SOLH(N,LDE1,E1,1,CONT,IP1)
DO 640 MM=1,N-2
MP=N-MM
DO 600 I=MP+1,N
600 CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
I=IPHES(MP)
IF (I.EQ.MP) GOTO 640
ZSAFE=CONT(MP)
CONT(MP)=CONT(I)
CONT(I)=ZSAFE
640 CONTINUE
END IF
END IF
ERR=0.D0
DO 364 I=1,N
364 ERR=ERR+(CONT(I)/SCAL(I))**2
ERR=DMAX1(DSQRT(ERR/N),1.D-10)
END IF
C --- COMPUTATION OF HNEW
C --- WE REQUIRE .2<=HNEW/H<=8.
FAC=DMIN1(SAFE,CFAC/(NEWT+2*NIT))
QUOT=DMAX1(.125D0,DMIN1(5.D0,ERR**.25D0/FAC))
HNEW=H/QUOT
C *** *** *** *** *** *** ***
C IS THE ERROR SMALL ENOUGH ?
C *** *** *** *** *** *** ***
IF (ERR.LT.1.D0) THEN
C --- STEP IS ACCEPTED
FIRST=.FALSE.
NACCPT=NACCPT+1
XOLD=X
HOLD=H
X=XPH
DO 75 I=1,N
Y(I)=Y(I)+Z3(I)
Z2I=Z2(I)
Z1I=Z1(I)
CONT(I+N)=(Z2I-Z3(I))/C2M1
AK=(Z1I-Z2I)/C1MC2
ACONT3=Z1I/C1
ACONT3=(AK-ACONT3)/C2
CONT(I+N2)=(AK-CONT(I+N))/C1M1
75 CONT(I+N3)=CONT(I+N2)-ACONT3
IF (ITOL.EQ.0) THEN
DO 81 I=1,N
81 SCAL(I)=ATOL(1)+RTOL(1)*DABS(Y(I))
ELSE
DO 82 I=1,N
82 SCAL(I)=ATOL(I)+RTOL(I)*DABS(Y(I))
END IF
IF (IOUT.NE.0) THEN
NRSOL=NACCPT+1
XSOL=X
XOSOL=XOLD
DO 77 I=1,N
77 CONT(I)=Y(I)
NSOLU=N
HSOL=HOLD
CALL SOLOUT(NRSOL,XOSOL,XSOL,CONT,NSOLU,IRTRN)
IF (IRTRN.LT.0) GOTO 79
END IF
CALJAC=.FALSE.
IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
H=HOPT
IDID=1
RETURN
END IF
CALL FCN(N,X,Y,Y0)
NFCN=NFCN+1
HNEW=POSNEG*DMIN1(DABS(HNEW),HMAXN)
HOPT=HNEW
IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
REJECT=.FALSE.
IF ((X+HNEW/QUOT1-XEND)*POSNEG.GT.0.D0) THEN
H=XEND-X
ELSE
QT=HNEW/H
HHFAC=H
IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) GOTO 30
H=HNEW
END IF
HHFAC=H
IF (THETA.LE.THET) GOTO 20
GOTO 10
ELSE
C --- STEP IS REJECTED
REJECT=.TRUE.
IF (FIRST) THEN
H=H*0.1D0
HHFAC=0.1D0
ELSE
HHFAC=HNEW/H
H=HNEW
END IF
IF (NACCPT.GE.1) NREJCT=NREJCT+1
IF (CALJAC) GOTO 20
GOTO 10
END IF
C --- UNEXPECTED STEP-REJECTION
78 CONTINUE
IF (IER.NE.0) THEN
WRITE (6,*) ' MATRIX IS SINGULAR, IER=',IER,' X=',X,' H=',H
NSING=NSING+1
IF (NSING.GE.6) GOTO 79
END IF
H=H*0.5D0
HHFAC=0.5D0
REJECT=.TRUE.
IF (CALJAC) GOTO 20
GOTO 10
C --- FAIL EXIT
79 WRITE (6,*) X,H,IER,NSTEP,IRTRN,NSING
979 FORMAT('X=',D14.7,' H=',D14.7,' IER=',I4,' NSTEP=')
IDID=-1
RETURN
END
C
C
REAL*8 FUNCTION CONTR5(I,X)
C ----------------------------------------------------------
C THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN
C APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X.
C IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
C THE LAST SUCCESSFULLY COMPUTED STEP (BY RADAU5).
C ----------------------------------------------------------
IMPLICIT REAL*8 (A-H,O-Z)
COMMON /CONT/NN,NN2,NN3,XSOL,HSOL,C2M1,C1M1,CONT(1)
S=(X-XSOL)/HSOL
CONTR5=CONT(I)+S*(CONT(I+NN)+(S-C2M1)*(CONT(I+NN2)
& +(S-C1M1)*CONT(I+NN3)))
RETURN
END
C