C ALGORITHM 689, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 17, NO. 2, PP. 167-177. JUNE, 1991. FILE_1: THIS INFO. FILE_2: SUBROUTINE COLDOC: DOCUMENTATION OF "COLVI2". FILE_3: SUBPROGRAMS OF "COLVI2"; "SOLVING" ROUTINES. FILE_4: SUBPROGRAMS OF "COLVI2"; "INITIALIZING" ROUTINES. FILE_5: SUBPROGRAMS OF "COLVI2"; "UTILITIES". FILE_6: SYSTEM DEPENDENT SUBPROGRAMS OF "COLVI2": SAVALL _ SAVE, IN CASE OF ERROR DURING COMPUTATIONAL PROCESS, ALL NEEDED VARIABLES ON A SEQUENTIAL UNFORMATTED FILE "COLSAV". RELOAD _ RELOAD VARIABLES WRITTEN BY "SAVALL" TO CONTINUE COMPUTATIONAL PROCESS AFTER ERRATIC EXIT. NCPJOB _ MONITOR THE CPU-TIME (CDC CYBER-750 VERSION). NCPJOB _ (VAX VERSION). INICMC _ INITIALIZE MACHINE CONSTANTS (CALLS "MACHAR" FROM W.J.CODY). FILE_7: ENVIRONMENT DEPENDENT SUBPROGRAMS OF "COLVI2" THAT INVOKE IMSL LIBRARY ROUTINES. DECLUF _ DECOMPOSE A MATRIX. SOLLUF _ SOLVE A LINEAR SYSTEM USING THE MATRIX DECOMPOSED BY DECLUF. ZERPOL _ COMPUTE REAL ZEROS OF A POLYNOMIAL. FILE_8: SAME AS ABOVE BUT FOR USE ON A MACHINE WITHOUT NUMERICAL LIBRARIES. ZERPOL CALLS THE (INCORPORATED) ACM-TOMS ROUTINE "RPOLY". FILE_9: MACHAR ROUTINE (REVISION DECEMBER 4, 1987) FROM W.J.CODY TO COMPUTE SOME MACHINE CONSTANTS. (MUST BE MODIFIED BEFORE COMPILING). FILE10: THREE UTILITY SUBROUTINES USED BY THE DRIVER PROGRAMS. FILE11: SUBROUTINES DEFINING THE VIE2 FOR PASS 1 - PASS 10. DRIVER PROGRAM PASS 1 - 6. DRIVER PROGRAM PASS 7 - 8. DRIVER PROGRAM PASS 9. DRIVER PROGRAM PASS 10 (RESUMING INTERRUPTED PASS 9). DRIVER PROGRAM PASS 11 + PROBLEM ROUTINES. DRIVER PROGRAM PASS 12 (RESUMING INTERRUPTED PASS 11). FILE12: OUTPUT PASS 1 - 6. OUTPUT PASS 7 - 8. OUTPUT PASS 9. OUTPUT PASS 10. OUTPUT PASS 11. OUTPUT PASS 12. N.B. RESULTS ARE OBTAINED ON A CYBER 750 (MACHINE PRECISION APPROXIMATELY 14 DIGITS). ------------------------------------------------------------------------ LIST OF SUBPROGRAM NAMES IN "COLVI2" + PARAMETER LIST (IN LEXICOGRAPHICAL ORDER). SUBROUTINE ADDABM (A, IA, IO, JO, N, S, B, IB) SUBROUTINE ADDABV (V, N, S, W) SUBROUTINE ADDV (V, N, S1, W1, S2, W2) SUBROUTINE ADJLSV (TN, HN, NEQN, KC, C,W, UN, LAGSAV) SUBROUTINE CHKFIL (CNTRL, IERROR) SUBROUTINE CHKINI + (NEQN,G, T0,TE, REQTOL, DEFOPT,IOPT,OPT,CNTRL, WKAREA,IW, + TN, HINIT, ZLEESM, IERROR) SUBROUTINE CHKOPT (CNTRL, IOPT, OPT, MCDEF, HMINFX, IERROR) SUBROUTINE CHKPTO (NEQN, TN, TE, REQTOL, TOLMIN, IERROR) SUBROUTINE CHKREC (IOPT, MCDEF, IERROR) SUBROUTINE CHKWKA (IW, TN, TE, HFX, IWCONS, IWSTEP, IERROR) SUBROUTINE COLCWL (COLPAR, M, C, W, LC, IERROR) SUBROUTINE COLDOC SUBROUTINE COLVI2 + (NEQN, G,KC,DKCDY,LINEAR, T0,TE, REQTOL, + DEFOPT, IOPT,OPT,CNTRL, WKAREA,IW, TNC, UE, GEE, IERROR) SUBROUTINE COMPLG (M, C, LCG, WKAREA) SUBROUTINE COMPLV (V, M, C, LCV) SUBROUTINE COMPUH (T, NEQN, T0, WKAREA, UH) SUBROUTINE COMPWL (M, S, L, C, W, LC) SUBROUTINE COPYV (V, N, W) SUBROUTINE DECLUF (A, N, IA, P, IERROR) SUBROUTINE DISWKS + (NEQN, NHC, RENTRY, NEWOPT, ZLEESM, WKAREA, IW, GECO, MAXNCO) SUBROUTINE ERRMSG (STRING) SUBROUTINE ESCRGS (NEQN, WKAREA,IW, T0,TE, TN, IERROR) SUBROUTINE GAUSSC (M, C, P0, P1, IERROR) SUBROUTINE GAUSS (M, C, WKAREA, IERROR) SUBROUTINE INICMC SUBROUTINE INILAG + (TN, NEQN,G,KC, C,W,LC,LOBAT, H, U, WKAREA, LAGSAV) SUBROUTINE INILGN (TN, NEQN, G, U, LC1, WKAREA, LAGN) SUBROUTINE INIVEC + (IU,IURG,IURL,IURN,ILEESM, NEQN, G, TN, TE, U,UR,URN,LEESUM) REAL FUNCTION INTEGL (J, U, M, C, WKAREA) SUBROUTINE ITRCOL + (TNP1, NEQN, G, KC, T0, H, C, W, U, WKAREA, LAGNP1, URNP1) REAL FUNCTION LAGPOL (J, V, M, C) REAL FUNCTION LEEWGT + (TN,HN, NEQN,KC, T0, C,W, CR,WR, UN, URN, WKAREA) SUBROUTINE LOBATC (M, C, P0, P1, IERROR) SUBROUTINE LOBATO (M, C, WKAREA, IERROR) INTEGER FUNCTION NCPJOB () SUBROUTINE RADAUC (M, C, P0, P1, IERROR) SUBROUTINE RADAU (M, C, WKAREA, IERROR) SUBROUTINE RELOAD + (NSAV, WKAREA,IW, DEFOPT,IOPT,OPT, TE, TN, IERROR) SUBROUTINE SAVALL (WKAREA,IW, DEFOPT, IOPT,OPT, TE, TN) SUBROUTINE SGEVI2 + (NEQN, G,KC,DKCDY,LINEAR, T0, CR,WR,LCR, H, UR, WKAREA, + TN, URN, IERROR) SUBROUTINE SLICE2 + (TN,HNM1,HN, NEQN,G,KC,DKCDY,LINEAR, C, CB, CR,WR,LCR,LOBATR, + LAGSAV, GLAGR,CORR,DSYS,MW, WKAREA, URN, IERROR) SUBROUTINE SLQCE2 + (TN, H, NEQN,G,KC,DKCDY,LINEAR, T0, C,W,LC, MM,SS,LL, LOBAT, + TOLCIT, LNP1FL, LAG, GLAG,CORR,DSYS,MW, WKAREA, U, LAGNP1, + IERROR) SUBROUTINE SOLLUF (A, N, IA, B, P) SUBROUTINE SOLNEW + (TN,HN, NEQN,KC,DKCDY,LINEAR, C,W,LC, MM,SS,LL, TOLCIT, + GLAG, CORR,DSYS,MW, WKAREA, UN, IERROR) SUBROUTINE SOLSYS + (TN,HN, NEQN,KC,DKCDY,LINEAR, C,W,LC, MM,SS,LL, TOLCIT, + GLAG, CORR,DSYS,MW, WKAREA, UN, IERROR) SUBROUTINE SOLVI2 + (NEQN, G,KC,DKCDY,LINEAR, T0,TE, C,W,LC, H, U, CR,WR,LCR, UR, + LAG, LAGSAV, LEESUM,ESTGEE,LEE, LC1, LC0,LCG,UN2,LAGN,LAGNP1 + URN, URNP1, WKAREA, TN, UN, GEE, IERROR) REAL FUNCTION UEEWGT + (TN,HN, NEQN,KC, T0,TE, C,W, CR,WR, UN, URN, + LEESUM, ESTGEE, WKAREA, LEE) SUBROUTINE UNITM (A, N) SUBROUTINE UTILIP (NEQN, UI, LC1, UIP1) REAL FUNCTION WMXNRM (ERR, SOL, NDIM) SUBROUTINE WRIRES (TN, HN, YNP1, UNP1, NEQN, Y) LOGICAL FUNCTION YPOLM + (TN,HN, NEQN,G,KC,DKCDY,LINEAR, C,W,LC, LCG,LC0, + LAGN,LAG,LAGNP1, UN2, NYPOLM, GLAG2,CORR,DSYS,MW, WKAREA, + UN, URN, IERROR) SUBROUTINE ZEROV (V, N) SUBROUTINE ZERPOL (C, N, S, IERROR) * ----------------------------------------------------------------------- SUBROUTINE COLDOC C C -------------------------------------------------------------------- I C I C HISTORY: 86/06/20: DATE WRITTEN. I C ======= 88/10/29: REVISION TOMS. I C I C -------------------------------------------------------------------- I C -------------------------------------------------------------------- I C I C DOCUMENTATION OF COLVI2 I C I C -------------------------------------------------------------------- I C I C SUBROUTINE COLVI2 I C + (NEQN, G,KC,DKCDY,LINEAR, T0,TE, REQTOL, I C + DEFOPT, IOPT,OPT,CNTRL, WKAREA,IW, TNC, UE, GEE, IERROR) I C I C PARAMETER SPECIFICATION: I C ----------------------- I C INTEGER NEQN, DEFOPT, IW, IERROR I C INTEGER IOPT(*), CNTRL(*) I C LOGICAL LINEAR I C REAL T0, TE, REQTOL, TNC I C REAL OPT(*), WKAREA(IW), UE(NEQN), GEE(NEQN) I C EXTERNAL G, KC, DKCDY I C I C -------------------------------------------------------------------- I C I C LANGUAGE: FORTRAN 77 I C ======== I C -------------------------------------------------------------------- I C I C PURPOSE: I C ======= I C SOLVE THE (SYSTEM OF) VOLTERRA INTEGRAL EQUATION(S) OF THE SECOND I C KIND (VIE2) I C T I C Y(T) = G(T) + INT K(T,S,Y(S)) DS, T IN [T0,TE]. I C T0 I C I C THE VIE2 IS APPROXIMATED BY A SYSTEM OF DISCRETIZED COLLOCATION I C EQUATIONS; THIS SYSTEM IS SOLVED BY AN ITERATION PROCESS THAT IS BY I C USER'S CHOICE FUNCTIONAL ITERATION OR NEWTON'S METHOD. I C I C IN THE ITERATIVE SOLUTION PROCESS EITHER A FIXED, USER-DEFINED I C OR A VARIABLE STEPSIZE IS USED. IN THE LATTER CASE AN ATTEMPT IS I C MADE TO KEEP THE GLOBAL ERROR IN "TE" LIMITED TO A PRESCRIBED I C TOLERANCE. I C FOR INFORMATION ON THE COLLOCATION SCHEMES, THE ERROR ESTIMATION AND I C STEPSIZE STRATEGY SEE: I C BLOM, J.G. AND BRUNNER, H., I C "THE NUMERICAL SOLUTION OF NONLINEAR VOLTERRA INTEGRAL EQUATIONS I C OF THE SECOND KIND BY COLLOCATION AND ITERATED COLLOCATION METHODS"I C SIAM J. SCI. STATIST. COMPUT. 8, 806-830 (1987). I C AND I C "DISCRETIZED COLLOCATION AND ITERATED COLLOCATION FOR NONLINEAR I C VOLTERRA INTEGRAL EQUATIONS OF THE SECOND KIND", REPORT NM-R8618, I C CWI, AMSTERDAM, 1986. (TO APPEAR IN ACM TOMS). I C I C -------------------------------------------------------------------- I C I C HOW TO USE: I C ========== I C INVOKE VIA A CALL OF THE SUBROUTINE "COLVI2". THE INPUT PARAMETERS I C ARE CHECKED ON LEGITIMACY AND CONSISTENCY; ERROR MESSAGES ARE I C WRITTEN TO A FILE (CF. PARAMETER "CNTRL(2)"). I C "COLVI2" OFFERS TWO DEFAULT SOLVERS FOR WHICH THE USER NEEDS TO I C SPECIFY ONLY THE PROBLEM DEPENDENT FUNCTIONS AND VARIABLES (PARAM- I C ETERS "NEQN" TO "TE"), THE REQUIRED TOLERANCE ("REQTOL"), THE PARAM- I C ETER "DEFOPT" TO SPECIFY THE SOLVER, WORKING STORAGE "WKAREA" OF I C SIZE "IW", AND THE OUTPUT PARAMETERS "TNC" TO "IERROR". ("IOPT", I C "OPT" AND "CNTRL" CAN BE DUMMY VARIABLES) I C I C PARAMETERS: I C ---------- I C NEQN DIMENSION OF THE SYSTEM OF VOLTERRA INTEGRAL EQUATIONS. I C G SUBROUTINE G(T,GV); REAL T, GV(NEQN). I C EVALUATES THE FORCING TERM "G" OF THE VIE2 IN "T". THE SECOND I C ARGUMENT OF THE SUBROUTINE IS AN ARRAY IN WHICH THE VALUE OF I C THE VECTOR "G(T)" SHOULD BE STORED. I C SHOULD BE DECLARED IN AN EXTERNAL STATEMENT IN THE CALLING I C PROGRAM. I C KC SUBROUTINE KC(T,S,Y,KV); REAL T,S, Y(NEQN),KV(NEQN). I C EVALUATES THE KERNEL "K" OF THE VIE2 IN "(T,S,Y)". THE LAST I C ARGUMENT OF THE SUBROUTINE IS AN ARRAY IN WHICH THE VALUE OF I C THE VECTOR "K(T,S,Y)" SHOULD BE STORED. I C SHOULD BE DECLARED EXTERNAL IN THE CALLING PROGRAM. I C DKCDY SUBROUTINE DKCDY(T,S,Y,DKV); REAL T,S, Y(NEQN),DKV(NEQN,NEQN).I C EVALUATES IN THE POINT "(T,S,Y)" THE JACOBIAN OF THE KERNEL I C W.R.T ITS THIRD ARGUMENT. THE LAST ARGUMENT OF THE SUBROUTINE I C IS A TWO-DIMENSIONAL ARRAY IN WHICH THE JACOBIAN SHOULD BE I C STORED; I.E., DKV(I,J) = DK_I/DY_J (T,S,Y) . I C SHOULD BE DECLARED EXTERNAL IN THE CALLING PROGRAM. I C ONLY NEEDED IF NEWTON'S METHOD IS USED TO SOLVE THE SYSTEM OF I C COLLOCATION EQUATIONS; IF FUNCTIONAL ITERATION IS USED A I C DUMMY ROUTINE SUFFICES. I C LINEAR LOGICAL VALUE. I C "TRUE" INDICATES THAT THE KERNEL IS A LINEAR FUNCTION W.R.T. I C ITS THIRD ARGUMENT. NOT USED IF FUNCTIONAL ITER. IS EMPLOYED. I C T0 REAL VALUE. I C LEFT ENDPOINT OF THE INTEGRATION INTERVAL. I C TE REAL VALUE. I C RIGHT ENDPOINT OF THE INTEGRATION INTERVAL. I C REQTOL REAL VALUE. I C REQUESTED TOLERANCE FOR GLOBAL ERROR (NOT USED IN CASE OF I C CONSTANT STEPSIZES). I C DEFOPT INTEGER VALUE. I C 0: USE NO DEFAULT SOLVERS, I C ?1: GAUSS 8 + ITERATED GAUSS COLLOCATION; I C ESCAPE IN CASE SOLUTION IS A POLYNOMIAL OF DEGREE < 8 TO:I C GAUSS 8 + [GAUSS 9 + (C10=1)] WITH I C LOCAL + UNIFORM ERROR CONTROL, I C ?2: LOBATTO 6 + LOBATTO 7 WITH I C GLOBAL AND UNIFORM ERROR CONTROL, I C WHERE "?" INDICATES THE METHOD OF CORRECTOR ITERATION: I C 0: NEWTON'S METHOD; UPDATE JACOBIAN EACH NEWTON ITERATION, I C 1: NEWTON'S METHOD; EVALUATE JACOBIAN ONCE PER STEP, I C 2: FUNCTIONAL ITERATION (NO "DKCDY" NEEDED!). I C IF DEFOPT > 0, THE DEFAULTS FOR "IOPT", "OPT" AND "CNTRL" I C VALUES ARE USED (IN SOFAR NOT CONTRADICTING THE ABOVE). I C IN THIS CASE THESE VECTORS ARE NOT USED SO THERE ARE NO I C RESTRAINTS WITH RESPECT TO THE VALUE OR THE LENGTH OF THESE I C ARRAYS. I C IOPT INTEGER VALUED OPTION VECTOR OF LENGTH AT MOST 9 (TRUE LENGTH I C DEPENDENT ON "DEFOPT", RESP. IOPT(2), IOPT(3)); I C (0): DEFAULT VALUES. I C 1. KIND AND NUMBER OF COLLOCATION POINTS. I C 0: 8 POINT GAUSS COLLOCATION (ORDER = 8), I C OTHER: VALUE WITH DECIMAL EXPANSION "MC", I C WHERE "C" SPECIFIES THE KIND AND "M" THE # OF I C COLLOCATION POINTS. I C C COLL. POINTS M GLOBAL ORDER I C 1 GAUSS 2<=M M I C 2 (M-1) GAUSS + (CM=1) 3<=M 2M-2 I C 3 LOBATTO 2<=M 2M-2 I C 4 RADAU 2<=M 2M-1 I C 2. STEPSIZE CHOICE. I C 0: VARIABLE STEPSIZE, I C 1: FIXED STEPSIZE. I C 3. GLOBAL ERROR IN ENDPOINT REQUIRED? I C 0: GLOBAL ERROR ESTIMATION IN "TE", I C 1: NO GLOBAL ERROR ESTIMATION. I C 4. DEFINES ERROR WEIGHTS. I C 0: MIXED ERROR (1 / MAX(1.0,!SOL!)), I C 1: ABSOLUTE ERROR (1.0), I C 2: RELATIVE ERROR (1 / !SOL!). I C 5. INDICATES METHOD OF CORRECTOR ITERATION IN THE PROCESS OF I C SOLVING THE SYSTEM OF COLLOCATION EQUATIONS. I C 0: NEWTON'S METHOD; UPDATE JACOBIAN EACH NEWTON ITERATION, I C 1: NEWTON'S METHOD; EVALUATE JACOBIAN ONCE PER STEP, I C 2: FUNCTIONAL ITERATION (NO "DKCDY" NEEDED!). I C 6. MAXIMUM # KERNEL EVALUATIONS ALLOWED. I C 0: NO MAXIMUM. I C 7. MAXIMUM # CPU-SECONDS ALLOWED. I C 0: NO MAXIMUM. I C 8. KIND AND NUMBER OF COLLOCATION POINTS FOR REFERENCE I C SOLUTION. (NEEDS ONLY TO BE SPECIFIED IF IOPT(2)=0 OR I C IOPT(3)=0). I C 0: 8 POINT ITERATED GAUSS COLLOCATION (ORDER = 16), I C OTHER: VALUE WITH DECIMAL EXPANSION "MC" WHERE "C" AND "M" I C CAN HAVE THE VALUES AS SPECIFIED UNDER IOPT(1) WITH THE I C EXCEPTION THAT THE VALUE "M1", WITH 2<=M, INDICATES I C M POINT ITERATED GAUSS COLLOCATION (ORDER = 2M). I C NOTE: ITERATED COLLOCATION IS ONLY ALLOWED IN COMBINATION I C ---- WITH GAUSS COLLOCATION. I C 9. STEPSIZE STRATEGY CONTROLLER (NEEDED ONLY IF IOPT(2)=0). I C IF REF.SOL. IS COMPUTED BY ITERATED COLLOC.: A VALUE I C WITH DECIMAL EXPANSION "PT" WHERE "P" CAN HAVE THE I C VALUE 0,1 OR 2 AND "T" 0 OR 1. I C P=0: CHECK IF SOLUTION IS A POLYNOMIAL OF DEGREE < M, I C AND IF SO, ESCAPE TO LOCAL + UNIFORM ERROR CONTROL I C WITH COMPUTATION OF THE REFERENCE SOLUTION BY I C (M+1) GAUSS + (C[M+2]=1). I C IF IOPT(3)=0, THE GLOBAL ERROR IN "TE" WILL BE I C ESTIMATED BY THE SUM OF THE LOCAL ERRORS IN "TE". I C P=1: CHECK IF SOLUTION IS A POLYNOMIAL OF DEGREE < M, I C IF SO, RETURN TO CALLING PROGRAM. I C P=2: NO CHECK. I C T=0: INCREASE TOLERANCE BY A FACTOR "TOLREL" (SEE BELOW I C UNDER "CONSTANTS USED") IF THE ERROR RESULTING I C FROM A STEP WITH VALUE "HMIN" IS GREATER THAN THE I C TOLERANCE. I C T=1: RETURN TO CALLING PROGRAM IF TOLERANCE CAN NOT BE I C SATISFIED. I C OTHERWISE: A VALUE WITH DECIMAL EXPANSION "PGUT", WHERE I C "P" AND "T" ARE DESCRIBED ABOVE; "P" IS USED ONLY IF I C THE APPROX. METHOD IS GAUSS COLLOCATION, IF G=0 AND IF I C ORDER REF.SOL. <= ORDER GAUSS QUADR. FORMULA. I C IN CASE P=0 AND IF THE SOLUTION BEHAVES AS A POLYNOMIAL I C AN ESCAPE IS MADE TO LOCAL+UNIFORM ERROR CONTROL WITH I C COMPUTATION OF THE REF.SOL. BY THE SAME METHOD AS I C BEFORE BUT WITH AN ADEQUATE # OF COLLOC. POINTS. IF I C IOPT(3)=0, THE GLOBAL ERROR IN "TE" WILL BE ESTIMATED I C BY THE SUM OF THE LOCAL ERRORS IN "TE". I C EACH OF THE DIGITS "G" AND "U" HAS THE VALUE 0 OR 1. I C G=0: STEPSIZE CONTROL USING GLOBAL ERROR. I C G=1: STEPSIZE CONTROL USING LOCAL ERROR. I C U=0: MODIFIED STEPSIZE CONTROL; UNIFORM ERROR CONTROL I C OVER REMAINING INTERVAL. I C U=1: NO MODIFICATION OF STEPSIZE CONTROL. I C NOTE: TO USE LOCAL OR UNIFORM ERROR CONTROL IN I C ---- COMBINATION WITH AN M-POINT GAUSS QUADRATURE I C FORMULA REQUIRES AN ORDER OF THE REFERENCE SOLUTION OF I C AT LEAST 2M+1. I C NOTE: BOTH THE VALUES IOPT(6) AND IOPT(7) CAN BE EXCEEDED BY THE I C ---- NUMBER OF KERNEL EVALUATIONS, RESP. CPU-SECONDS NEEDED TO I C SOLVE THE SYSTEM OF COLLOCATION EQUATIONS IN AN INTERVAL. I C OPT REAL VALUED OPTION VECTOR OF LENGTH AT MOST 4 (TRUE LENGTH I C DEPENDENT ON "DEFOPT", RESP. IOPT(2), IOPT(9)); I C (0.0): DEFAULT VALUES. I C 1. INITIAL VALUE FOR STEPSIZE. IF IOPT(2)=1 THIS LENGTH I C WILL REMAIN FIXED THROUGHOUT THE COMPUTATION AND SHOULD I C NOT BE ZERO. I C 0.0: DEFAULT VALUE. IF FIRST CALL: HINIT = MIN(TE-T0,1.0), I C OTHERWISE HINIT IS SET TO THE GUESS OF THE LENGTH OF I C THE NEXT SUBINTERVAL MADE IN THE PREV. CALL OF COLVI2.I C 2. MINIMUM STEPSIZE. (NEEDED ONLY IF IOPT(2)=0). I C 0.0: DEFAULT VALUE HMIN = MAX(SUNFLO,HINIT*HMINFC), I C WITH HMINFC = 1E-5 (SEE BELOW UNDER "CONSTANTS USED").I C MIN. LENGTH OF N-TH SUBINTERVAL IS: I C HMINN = MAX(HMIN,SRELPR*!TN!) I C WITH SUNFLO = SMALLEST F.P. NUMBER AND I C SRELPR = F.P. MACHINE PRECISION, I C (SEE UNDER "MACHINE CONSTANTS"). I C 3. MAXIMUM STEPSIZE. (NEEDED ONLY IF IOPT(2)=0). I C 0.0: DEFAULT VALUE HMAX = 1.0. I C IF "COLVI2" IS REQUIRED TO CHECK IF THE SOLUTION BEHAVES I C AS A POLYNOMIAL HMAX SHOULD BE <= 1.0. I C 4. INTERVAL LENGTH HC. (NEEDED ONLY IF UNIFORM ERROR CONTROL I C IS REQUIRED); AT STEP "TN" LOCAL ERR. CONTR. IS PERFORMED I C IN TE, (-HC), TN+HN. I C 0.0: DEFAULT VALUE HC = HMAX. I C CNTRL INTEGER VALUED CONTROL VECTOR OF LENGTH AT MOST 4. I C (IF "DEFOPT">0 "IOPT", "OPT" AND "CNTRL" ARE NOT USED). I C (0): DEFAULT VALUES. I C 1. RE-ENTRY INDICATOR (CF. "GENERAL COMMENTS" SUB "ENTRY"). I C 0: FIRST ENTRY, I C 1: RE-ENTRY, NEW OPTIONS, I C 2: RE-ENTRY AFTER SAVE, NEW OPTIONS, I C 3: RE-ENTRY, OLD OPTIONS, I C 4: RE-ENTRY AFTER SAVE, OLD OPTIONS. I C 2. LOGICAL UNIT NUMBER OF FILE FOR ERROR MESSAGES. I C (CF. "GENERAL COMMENTS" SUB "ERROR MESSAGES"). I C 0: ERROR MESSAGES ARE WRITTEN TO THE STANDARD OUTPUT FILE I C (ADDRESSED BY THE "PRINT" STATEMENT). I C 1 <= CNTRL(2) <= IMXLUN (SEE "MACHINE CONSTANTS"). I C 3. CONTROL ON WRITING OF RESULTS IN ALL STEP POINTS. I C (CF. "GENERAL COMMENTS" SUB "WRITE ALL"). I C 0: NO INTERMEDIATE WRITING. I C 1 <= CNTRL(2) <= IMXLUN : LOG. UNIT # OF FILE. I C 4. INDICATOR TO SAVE VARIABLES FOR RE-ENTRY AFTER ERROR. I C (CF. "GENERAL COMMENTS" SUB "SAVE"). I C 0: NO SAVE. I C 1 <= CNTRL(2) <= IMXLUN : LOG. UNIT # OF FILE TO I C SAVE COMMON BLOCKS, WKAREA, IOPT AND OPT FOR RE-ENTRY. I C WKAREA REAL WKAREA(IW). I C TEMPORARY WORKING STORAGE. THE EXACT AMOUNT NEEDED IS I C SPECIFIED BELOW UNDER "DISTRIBUTION WKAREA" AND "STORAGE I C OCCUPIED". HERE, A FEW SPECIFIC CASES ARE TREATED. I C FIRST AN UPPERBOUND WILL BE GIVEN FOR THE CASE "DEFOPT" > 0. I C LET MAXNC = MAX. # (SUCCESSFUL) STEPS POSSIBLE (THIS NUMBER I C WILL BE CALCULATED BY "COLVI2" FROM THE USER I C SUPPLIED INPUT VALUES), I C THEN AN UPPERBOUND FOR THE DIMENSION OF "WKAREA" IS GIVEN BY: I C IF "DEFOPT" = ?1 I C 1456+101.NEQN.NEQN+(61+2.INT(TE-T0)).NEQN+(1+8.NEQN).MAXNC,I C IF "DEFOPT" = ?2 I C 429+37.NEQN.NEQN+(23+2.INT(TE-T0)).NEQN+(1+11.NEQN).MAXNC. I C NOTE THAT FOR "DEFOPT" = 0, IT IS POSSIBLE TO JUST TAKE SOME I C VALUE FOR "IW", SAY 2000, SET "CNTRL(4) > 0" AND LET THE I C PROGRAM RUN. IF THE CODE LACKS WORKING SPACE IT DUMPS ALL I C NEEDED VARIABLES ON FILE AND AN ERROR EXIT IS TAKEN. THE I C COMPUTATION CAN THEN BE RESTARTED FROM THE LAST REACHED POINT I C BY A NEW CALLING PROGRAM WHICH OFFERS MORE WORKING SPACE TO I C "COLVI2". I C NOW SOME EXAMPLES IN WHICH "IW" IS CALCULATED EXACTLY: I C LET M = # COLLOCATION PARAMETERS (SEE IOPT(1)), I C S = # QUADRATURE POINTS (IF IOPT(1)=?2: M-1, OTHERWISE M),I C L = 2 IF IOPT(1)=?3, OTHERWISE 1, I C AND DEFINE ANALOGOUSLY MR, SR, LR FROM IOPT(8). I C LET X = MAX(M,MR), I C NHC = INT((TE-T0)/HC), I C NIRVEC = 2.NEQN IF FUNCTIONAL ITERATION IS USED, I C OTHERWISE: NEQN.(1+NEQN), I C NWKSYS(N) = DIM.WKAREA SOLVER_FOR_LIN.SYS.OF_DIM.(N.NEQN), I C F(N) = IF FUNCT.IT.: MAX(NIRVEC,2.N-N.NEQN), I C IF NEWTON'S METHOD IS USED AND THE JACOBIAN IS I C UPDATED EACH ITERATION: I C (N.NEQN)**2 + MAX(NIRVEC,NWKSYS(N)), I C IF NEWTON'S METH. IS USED WITHOUT UPDATING THE JAC.I C (N.NEQN)**2 + NIRVEC+NWKSYS(N), I C THEN THE DIMENSION OF WKAREA SHOULD BE I C IF DEFOPT = ?1: I C 625+30.NEQN + (1+8.NEQN).MAXNC + I C (NO ESCAPE:) F(8) I C (ESCAPE:) 831+(20+2.NHC).NEQN+F(10). I C IF DEFOPT = ?2: I C 429+(16+2.NHC).NEQN+F(6)+ (1+11.NEQN).MAXNC. I C IN CASE OF CONSTANT STEPSIZES, S=M,L=1, NO ERROR EST. IN "TE":I C 1+3.M+M.M.M+(1+2.M).NEQN+F(M) + (1+M.NEQN).MAXNC. I C IN CASE OF ITERATED COLLOCATION WITHOUT ESCAPE: I C 1+6.M+M.M+M.M.M+(6+3.M).NEQN+F(M) + (1+M.NEQN).MAXNC. I C IN CASE OF LOC.+UNIF.ERR.CONTR.(NO GAUSS, NO ERR.EST.IN "TE"):I C 1+4.X+2.X.X.X+(4+5.X+2.NHC).NEQN+F(X) + I C (1+(M-L+1).NEQN).MAXNC. I C IN CASE OF GLOBAL + UNIFORM ERROR CONTROL (NO GAUSS): I C 1+4.X+2.X.X.X+(4+2.X+2.NHC).NEQN+F(X) + I C (1+(M-L+1+MR-LR+1).NEQN).MAXNC. I C IW INTEGER VALUE. I C DIMENSION OF "WKAREA" AS DECLARED IN MAIN PROGRAM. I C TNC REAL VALUE. I C EXIT: RIGHT ENDPOINT OF LAST INTERVAL ON WHICH SOLUTION HAS I C BEEN COMPUTED (SHOULD BE EQUAL TO "TE"). I C UE REAL UE(NEQN). I C EXIT: COMPUTED SOLUTION OF THE VIE2 AT "TE". I C GEE REAL GEE(NEQN). I C EXIT: CONTAINS THE GLOBAL ERROR EST. AT "TE" IF EITHER I C IOPT(3)=0 OR IF GLOBAL OR UNIFORM ERROR CONTROL I C HAS BEEN SPECIFIED. I C IERROR INTEGER VALUE. I C ERROR COMPLETION CODE. I C IF AN ERROR HAS BEEN DETECTED, INFO IS WRITTEN TO THE I C ERROR_MESSAGE_FILE (LOG. UNIT # : CNTRL(2)). I C 0: NO ERRORS. I C 1: "DEFOPT" INCORRECT OR I C INITIAL FILE STATUS WRONG OF ERROR_MESSAGE_FILE, I C INFO ON STANDARD OUTPUT FILE. I C 2: INCORRECT INPUT. I C INFO ON ERROR_MESSAGE_FILE. I C 3: FAILED TO COMPUTE COLLOCATION PARAMETERS. I C 11: FAILURE TO MEET TOLERANCE "REQTOL" WITH STEPSIZE "HMINN".I C 12: WORKING STORAGE NEEDED EXCEEDS "IW" (VAR. STEPSIZE). I C 13: TOTAL # KERNEL EVALUATIONS USED > IOPT(6). I C 14: TOTAL # CPU-SECONDS USED > IOPT(7). I C 15: POLYNOMIAL SOLUTION (GAUSS). I C 16: TOLERANCE WOULD BE RELAXED TO A VALUE > 1.0. I C 21: CORRECTOR ITER. PROCESS DID NOT CONVERGE WITHIN "MAXFIT" I C (FUNCTIONAL ITERATION) OR "MAXNIT" (NEWTON'S METHOD) I C ITERATIONS (CF. "CONSTANTS USED"). (FIXED STEPSIZE). I C 31: CORRECTOR ITER. PROCESS DID NOT CONVERGE WITHIN "MAXFIT" I C (FUNCTIONAL ITERATION) OR "MAXNIT" (NEWTON'S METHOD) I C ITERATIONS (CF. "CONSTANTS USED"). (WITHIN "YPOLM"). I C 113: TOTAL # KERNEL EV. USED > IOPT(6) -+ WHILE COMP. I C 114: TOTAL # CPU-SEC. USED > IOPT(7) ---+ GLOB. ERR. EST. TE.I C OTHER: ERROR COMPLETION CODE FROM ONE OF THE ROUTINES I C "DECLUF", "SOLLUF" OR "ZERPOL". I C I C -------------------------------------------------------------------- I C I C GENERAL COMMENTS: I C ================ I C I C EXAMPLE PROGRAMS: WITH THE PACKAGE GOES A SET OF DRIVER PROGRAMS TO I C ---------------- DEMONSTRATE THE USE OF "COLVI2" AND THREE UTILITY I C ROUTINES USED BY THE DRIVERS. FOR AN ELABORATE DESCRIPTION OF THE I C DRIVERS AND THE VIE2'S THEY SOLVE WE REFER TO THE DOCUMENTATION IN I C THE PROGRAMS. A SHORT DESCRIPTION OF THE DRIVERS AND THE UTILITY I C ROUTINES FOLLOWS: I C I C FIRST DRIVER PROGRAM: I C PASS 1 - PASS 6: DEMONSTRATE SIMPLE USE OF "COLVI2". I C USE DEFOPT = 21,22, 1,2 OR SMALL CHANGES ON THE DEFAULT VALUES. I C SECOND DRIVER PROGRAM: I C PASS 7 - PASS 8: DEMONSTRATE RE-ENTRY FACILITY OF "COLVI2". I C WRITE INTERMEDIATE RESULTS TO FILE TO SHOW TRANSITION. I C 7: SAME OPTIONS AS PASS 5; DIVIDE INTEGRATION INTERVAL IN TWO I C PARTS. I C 8: AS PASS 6 BUT WITHOUT AUTOMATIC ESCAPE; NEEDS MORE WORKING I C SPACE BECAUSE GLOB.ERR.EST. IN "TE" USES SEPARATELY COMPUTED I C REF.SOL. I C THIRD DRIVER PROGRAM: I C PASS 9: DEMONSTRATE SAVE FACILITY OF "COLVI2". I C SAME OPTIONS AS PASS 5; WORKING STORAGE DIMINISHED TO FORCE I C AN EXIT BECAUSE OF A LACK OF WORKING STORAGE. I C FOURTH DRIVER PROGRAM: I C PASS 10: DEMONSTRATE RE-ENTRY_AFTER_SAVE FACILITY OF "COLVI2". I C SAME OPTIONS AS PREV. CALL; WORKING STORAGE ENLARGED. I C FIFTH DRIVER PROGRAM: I C PASS 11: DEMONSTRATE SAVE FACILITY OF "COLVI2" AFTER ARITHMETIC I C ERROR OTHER PROBLEM, LOOOONG INTEGR. PATH.; USE FUNCTIONAL ITER.I C NO DKCDY NEEDED. EMPLOY RADAU INTEGRATION; HMAX = 10.0 AND I C UNIFORM ERROR CONTROL WITH HC=10.0 (=HMAX). I C FOR THIS DEMONSTRATION WE NEED A SYSTEM ROUTINE THAT ALLOWS I C THE USER TO REGAIN CONTROL AFTER ARITHMETIC MODE ERRORS. I C SIXTH DRIVER PROGRAM: I C PASS 12: DEMONSTRATE SAVE FACILITY OF "COLVI2" AFTER ARITHMETIC I C ERROR IN PASS 11; ENLARGE HMAX TO 50., NO UNIFORM ERROR CONTROL.I C I C THE THREE UTILITY ROUTINES: I C REAL FUNCTION AERE (Y, YA) I C REAL Y, YA I C COMPARISON OF THE APPROXIMATED SOLUTION "YA" WITH THE EXACT I C SOLUTION "Y". IF Y > 1.0 THE NUMBER OF CORRECT SIGNIFICANT I C DIGITS IS RETURNED, OTHERWISE THE NUMBER OF CORRECT DIGITS. I C SUBROUTINE SUMARY (NOUT, NEQN, WKAREA, YE, T, UE, GEE, IERROR) I C INTEGER NOUT, NEQN, IERROR I C REAL T I C REAL WKAREA(*), YE(NEQN), UE(NEQN), GEE(NEQN) I C EXTRACT STATISTICS FROM "COLVI2" COMMON BLOCKS AND WRITE SUMMARYI C OF RESULTS TO FILE WITH LUN "NOUT". I C SUBROUTINE ACVSUM (IERROR, WKAREA, TN, T0) I C INTEGER IERROR I C REAL TN, T0 I C REAL WKAREA(*) I C ENTRY SCVSUM I C ACCUMULATE COUNTING VALUES IN "COLVI2" COMMON BLOCK /COLCMI/ I C THAT ARE ZEROED WHEN "COLVI2" IS CALLED MORE THAN ONCE. I C ENTRY SCVSUM: STORE ACCUMULATED VALUES IN /COLCMI/. I C I C SOLUTION IN ARBITRARY POINT: THE PACKAGE CONTAINS A SUBROUTINE TO I C --------------------------- COMPUTE THE APPROXIMATION "UH" IN AN I C ARBITRARY POINT "T" BETWEEN "T0" AND "TNC" BY LAGRANGE INTERPOL. I C USING THE ARRAY OF APPROXIMATIONS "U" STORED IN "WKAREA". I C THE HEADING OF THIS ROUTINE IS: I C SUBROUTINE COMPUH (T, NEQN, T0, WKAREA, UH) I C INTEGER NEQN I C REAL T, T0 I C REAL WKAREA(*), UH(NEQN) I C "WKAREA" SHOULD CONTAIN THE ENTIRE, UNALTERED, "WKAREA" ARRAY AS I C RETURNED BY "COLVI2". I C I C ERROR MESSAGES: ERROR MESSAGES GENERATED BY "COLVI2" ARE WRITTEN I C -------------- TO A SEQUENTIAL FORMATTED FILE. IF DEFOPT > 0 OR I C CNTRL(2)=0 THE STANDARD OUTPUT FILE IS USED. OTHERWISE CNTRL(2) I C DEFINES THE LOGICAL UNIT NUMBER OF THE FILE AND THIS FILE SHOULD I C BE OPENED IN THE MAIN PROGRAM. I C I C THE NEXT THREE PARAGRAPHS ARE IRRELEVANT WHEN "COLVI2" WILL BE CALLEDI C WITH DEFOPT > 0. I C I C SAVE: THE PACKAGE HAS THE OPTION (CNTRL(4)>0) TO SAVE ALL NECESSARY I C ---- VARIABLES ON A FILE IN CASE OF AN ERROR DURING THE I C COMPUTATIONAL PROCESS, I.E. AFTER THE CONTROL AND INITIALIZATION I C PHASE. I C THE COMMON BLOCK VARIABLES, WORKING STORAGE AND OPTION VECTORS I C ARE WRITTEN TO A SEQUENTIAL UNFORMATTED FILE NAMED "COLSAV". THE I C LOGICAL UNIT NUMBER USED IS "CNTRL(4)". IF THE FILE "COLSAV" I C ALREADY EXISTS IT WILL BE OVERWRITTEN. I C THIS FILE SHOULD BE AVAILABLE IF "COLVI2" IS CALLED WITH I C CNTRL(1)=2 OR 4. I C (SEE ALSO UNDER "CONSTANTS USED" AND "OTHER MACHINE DEPENDENCIES") I C I C ENTRY: "COLVI2" ACKNOWLEDGES A NUMBER OF DIFFERENT ENTRY OPTIONS. I C ----- THE FIRST TIME "COLVI2" IS CALLED TO SOLVE A SPECIFIC VIE2 I C CNTRL(1) SHOULD BE 0 AND ALL OPTION VECTORS SHOULD HAVE LEGITIMATE I C VALUES. I C IT IS POSSIBLE TO CALL "COLVI2" A SECOND TIME IN THE SAME MAIN I C PROGRAM TO CONTINUE THE PROCESS OF SOLVING THE VIE2 AFTER A NORMAL I C EXIT OR AFTER "COLVI2" RETURNED WITH "IERROR" >= 10. IN THE LATTER I C CASE THE USER SHOULD REACT APPROPRIATELY ON THE GIVEN ERROR. I C THE ARRAY "WKAREA" SHOULD BE UNCHANGED OR COPIED INTO A NEW I C WORKING STORAGE; "TNC" SHOULD BE UNALTERED. I C CNTRL(1)=1 INDICATES THAT NEW OPTION VECTORS HAVE BEEN DEFINED, I C CNTRL(1)=3 THAT THE OLD OPTIONS SHOULD BE USED. I C A RE-ENTRY TO CONTINUE IN THE SAME OR A NEW JOB IS ALSO POSSIBLE I C AFTER AN ERRATIC EXIT (WITH SAVE CONTROL ON), BOTH WITH NEW OPTION I C VECTORS (CNTRL(1)=2) AND WITH THE OLD OPTIONS (CNTRL(1)=4). I C THE FILE "COLSAV" (SEE ABOVE) SHOULD BE AVAILABLE, AND "IW" SHOULD I C BE >= DIMENSION ARRAY "WKAREA" IN THE PREVIOUS CALL. I C NOTE: IN CASE OF RE-ENTRY THE NUMBER OF COLLOCATION PARAMETERS AND I C ---- THE COLLOCATION METHOD TO APPROXIMATE THE SOLUTION (IOPT(1)) I C SHOULD BE THE SAME AS IN THE PREVIOUS CALL OF "COLVI2". IF GLOBAL I C ERROR CONTROL IS USED THIS HOLDS ALSO FOR IOPT(8). I C IT IS NOT ALLOWED TO ASK FOR GLOBAL ERROR CONTROL IF THIS WAS NOT I C USED IN THE PREVIOUS CALL OF COLVI2 . I C IT IS NOT POSSIBLE TO RE-ENTER THE PROCESS OF COMPUTING THE GLOBAL I C ERROR ESTIMATE IN "TE" (HOWEVER, ONE CAN COPY THE CALL OF "SGEVI2" I C FROM "COLVI2"). I C IF UNIFORM ERROR CONTROL IS REQUIRED AND IN THE PREVIOUS CALL THIS I C EITHER WAS NOT THE CASE, OR THE ENDPOINT "TE" OR THE VALUE OF I C OPT(4) WERE DIFFERENT, THEN AN ESTIMATION OF THE SUM OF THE LOCAL I C ERRORS IS MADE BASED ON THE LOCAL ERRORS IN THE FIRST POINT I C COMPUTED. I C I C WRITE ALL: IF CNTRL(3) IS NON-ZERO THE RESULTS IN EACH STEP POINT I C --------- ARE WRITTEN TO A SEQUENTIAL FORMATTED FILE WITH LOGICAL I C UNIT NUMBER CNTRL(3). THIS FILE SHOULD BE OPENED IN THE CALLING I C PROGRAM. I C A SUBROUTINE YEXACT(T,YV) SHOULD EXIST AND DELIVER THE EXACT I C SOLUTIONS AT "T" IN YV(1:NEQN). IF NO SOLUTION IS AVAILABLE I C "YEXACT" STILL HAS TO BE PROVIDED AND SHOULD RETURN LEGITIMATE F.P.I C NUMBERS, SAY 0.0, IN YV. I C I C -------------------------------------------------------------------- I C I C SUBPROGRAMS: I C =========== I C I C SOLVING ROUTINES: I C ---------------- I C COLVI2 _ ENVELOPING ROUTINE. I C CHECK INPUT, INITIALIZE COMMON AND COLL. VARIABLES, I C DIGEST RESULTS OF "SOLVI2". I C SOLVI2 _ SUPERVISOR. I C SOLVE VIE2 WITH CHOSEN COLLOCATION METHOD. MONITOR STEPS, I C PERFORM STEP ACCEPTANCE / REJECTANCE. I C SGEVI2 _ COMPUTE REF.SOL. IN "TE" WITH CHOSEN COLL. METHOD AND STEP- I C SIZES AS USED IN THE COMPUTATION OF THE APPROXIMATION. I C SLQCE2 _ SET UP SYSTEM OF COLLOC. EQUATIONS IN A SUBINTERVAL, I C COMPUTE LAG TERMS WITH QUADRATURE. I C SLICE2 _ SET UP SYSTEM OF COLLOC. EQUATIONS IN A SUBINTERVAL, I C COMPUTE LAG TERMS WITH INTERPOLATION. I C SOLSYS _ SOLVE THE SYSTEM OF COLLOCATION EQUATIONS IN A SUBINTERVAL. I C SOLNEW _ SOLVE THE SYSTEM OF COLLOCATION EQUATIONS IN A SUBINTERVAL, I C CALLED BY "SOLSYS" IF NEWTON ITERATION IS TO BE USED. I C I C YPOLM _ CHECK IF SOLUTION BEHAVES AS A POLYNOMIAL OF DEGREE < M. I C I C UTILIP _ COMPUTE COLL. APPR. IN A STEPPOINT WITH LAGR. INTERPOL. I C COMPUH _ COMPUTE COLL. APPR. IN AN ARBITRARY POINT WITH LAGR. INTERP.I C ITRCOL _ COMPUTE THE ITERATED COLL. APPROX. IN A STEPPOINT. I C I C INITIALIZING ROUTINES: I C --------------------- I C CHKINI _ CHECK INPUT PARAMETERS AND INITIALIZE COMMON BLOCKS. I C DISWKS _ DISTRIBUTE WORKING STORAGE. I C INIVEC _ INITIALIZE VECTORS WITH SOLUTION AND SUM OF LOCAL ERRORS. I C ESCRGS _ ALTER COMMON BLOCK VALUES AND RE-DISTRIBUTE WORK SPACE IN I C CASE OF ESCAPE WHEN SOLUTION BEHAVED AS A POLYNOMIAL. I C INILAG _ INITIALIZE LAG TERM VECTORS IN CASE LAG TERMS REF. SOL. ARE I C COMPUTED BY INTERPOLATION. I C INILGN _ INITIALIZE LAG TERM VECTOR IN STARTING POINT IN CASE IT HAS I C TO BE CHECKED IF SOLUTION IS POLYNOMIAL. I C I C COLCWL _ INITIALIZE THE SET OF COLLOCATION PARAMETERS, I C THE ASSOCIATED WEIGHT FACTORS FOR THE QUADRATURE FORMULA I C AND THE LAGRANGIAN INTERPOLATION COEFFICIENTS NEEDED. I C GAUSS _ COMPUTE GAUSS-LEGENDRE COLL. PARAMETERS (2 <= M). I C GAUSSC _ COMPUTE GAUSS-LEGENDRE COLL. PARAMETERS (6 <= M). I C LOBATO _ COMPUTE LOBATTO COLL. PARAMETERS (2 <= M). I C LOBATC _ COMPUTE LOBATTO COLL. PARAMETERS (8 <= M). I C RADAU _ COMPUTE RADAU COLL. PARAMETERS (2 <= M). I C RADAUC _ COMPUTE RADAU COLL. PARAMETERS (4 <= M). I C COMPWL _ COMPUTE WEIGHT FACTORS FOR QUADRATURE FORMULA AND I C LAGRANGIAN INTERPOLATION COEFFICIENTS IN (CI.CJ). I C COMPLV _ COMPUTE LAGR. INTERPOL. COEFF. IN A GIVEN POINT. I C COMPLG _ GIVEN THE POINTS 0.0, C(1:M), 1.0 COMPUTE THE LAGR. INTERP. I C COEFF. IN CJ/2 ("C" ARE THE GAUSS COLL. PAR.). I C INTEGL _ U I C COMPUTE INT L_J(V) DV . I C 0 I C LAGPOL _ M I C COMPUTE L_J(V) = PROD (V-CI)/(CJ-CI) . I C I=1,I/=J I C I C UTILITIES: I C --------- I C ADJLSV _ ADJUST LAG TERM VECTORS IN CASE REF.SOL. LAGTERMS ARE COMP. I C BY INTERPOLATION. I C I C LEEWGT _ COMPUTE (TN+HN-T0)/HN * LOCAL ERROR IN (TN+HN). I C UEEWGT _ COMPUTE MAXIMUM OF GLOBAL ERRORS IN [TN,TE]. THE GLOBAL I C ERROR IS APPROXIMATED BY A SUM OF LOCAL ERRORS ON [T0,TN] I C AND AN ESTIMATION OF THE LOCAL ERRORS ON [TN+HN,TE] BASED I C ON THE LOCAL ERROR IN (TN+HN). I C WMXNRM _ COMPUTE THE MAXIMUM NORM OF A GIVEN ERROR VECTOR, USING I C WEIGHTS AS PRESCRIBED BY IOPT(4). I C I C THE NEXT FOUR SUBROUTINES INSPECT THE STARTING CONDITIONS FOR I C "COLVI2". IF THE INITIAL FILE STATUS OF THE ERROR_MESSAGE_FILE IS I C WRONG, A MESSAGE IS WRITTEN TO THE STANDARD OUTPUT FILE AND "COLVI2" I C RETURNS. OTHERWISE EACH SUBROUTINE CHECKS AS MUCH AS POSSIBLE IN ITS I C FIELD. ALL THE ERRORS FOUND, AS WELL AS THE ERRORS FOUND DURING AN I C EVENTUAL RELOAD ARE WRITTEN TO THE ERROR_MESSAGE_FILE. I C CHKFIL _ CHECK STATUS OF FILES THAT THE USER SHOULD HAVE OPENED. I C CHKPTO _ CHECK DIMENSION, INTEGRATION BOUNDS OF VIE2 AND TOLERANCE I C PARAMETERS, AND THE ORDER OF THE REF.SOL.METHOD. I C CHKOPT _ CHECK VALIDITY OF OPTION- AND CONTROL-VECTORS. I C CHKREC _ CHECK CONSISTENCY PARAMETERS WITH PREVIOUS CALL OF COLVI2. I C CHKWKA _ CHECK (AS FAR AS POSSIBLE) SIZE OF WORKING STORAGE AREA. I C I C ERRMSG _ WRITE ERROR MESSAGE TO A FILE WITH LOGICAL UNIT # CNTRL(2). I C WRIRES _ WRITE INTERMEDIATE RESULTS TO FILE WITH LOGICAL I C UNIT NUMBER CNTRL(3). I C I C ADDABM _ A = A + B, A AND B MATRICES. I C ADDABV _ V = V+W, V AND W VECTORS. I C ADDV _ V = W1+W2, V, W1 AND W2 VECTORS. I C COPYV _ COPY VECTOR. I C UNITM _ INITIALIZE MATRIX ON UNIT MATRIX. I C ZEROV _ ZERO VECTOR. I C I C THE NEXT SEVEN ROUTINES ARE SYSTEM OR ENVIRONMENT DEPENDENT. FOR A I C DESCRIPTION OF THESE SEE BELOW UNDER "OTHER MACHINE DEPENDENCIES". I C SAVALL _ SAVE, IN CASE OF ERROR DURING COMPUTATIONAL PROCESS, I C ALL NEEDED VARIABLES ON A SEQ. UNFORM. FILE "COLSAV". I C (COMMON, WKAREA, DEFOPT, IOPT, OPT, TE, TN). I C RELOAD _ RELOAD VARIABLES WRITTEN BY "SAVALL" TO CONTINUE I C COMPUTATIONAL PROCESS AFTER ERRATIC EXIT. I C NCPJOB _ RETURN NUMBER OF CPU SECONDS USED IN THIS JOB. I C INICMC _ INITIALIZE THE COMMON BLOCKS WITH MACHINE CONSTANTS. I C DECLUF _ DECOMPOSE FULL MATRIX INTO "LU". I C SOLLUF _ SOLVE LU.X = B . I C ZERPOL _ COMPUTE REAL ZEROS OF POLYNOMIAL AND SORT THESE. I C I C ---------------------------------------------------------------------I C I C DESCRIPTION OF VARIABLES AND CONSTANTS: I C ====================================== I C I C COMMON BLOCKS: I C ------------- I C THROUGHOUT THE PACKAGE TWO NAMED COMMON BLOCKS ARE USED THAT HOLD I C MACHINE CONSTANTS, I C I C ONE CONTAINING INTEGER VARIABLES: I C COMMON /COLMCI/ IBETA, IOVFLO, NSDEC, IMXLUN I C I C ONE CONTAINING FLOATING POINT VARIABLES: I C COMMON /COLMCR/ SRELPR, SOVFLO, SUNFLO I C I C FOUR OTHER NAMED COMMON BLOCKS ARE USED THAT HOLD METHOD PARAMETERS I C AND STATISTICS, I C I C ONE CONTAINING INTEGER VARIABLES: I C COMMON /COLCMI/ METH, M, S, L, ORDER, ORDERQ, I C + METHR, MR, SR, LR, ORDERR, I C + ERRWGT, NHFAIL, I C + NERR, NWIR, NSAV, I C + MAXNC, MAXKEV, MAXCPS, I C + N, NCIT, NKEV, NCPS I C I C ONE CONTAINING LOGICAL VARIABLES: I C COMMON /COLCML/ VS, GSSCKM, ESCGSS, GEC, ULEC, RLXTOL, GEETE, I C + FUNCIT, NEWTON I C I C A THIRD CONTAINING FLOATING POINT VARIABLES: I C COMMON /COLCMR/ TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC I C I C THE FOURTH CONTAINING INDEX VARIABLES FOR THE WORK SPACE "WKAREA" I C COMMON /COLIXW/ IC1,IC2,IC3,IC4,IC5,IC6,IC7,IC8,ICE, I C + IV1,IV2,IV3,IVE, IL1,IL2,IL3,IL4,IL5,ILAG,IL6,ILE I C I C COMMON VARIABLES: I C ---------------- I C FOR A DESCRIPTION OF THE VARIABLES IN THE COMMON BLOCKS /COLMCI/ AND I C /COLMCR/ SEE BELOW UNDER "MACHINE CONSTANTS". I C INTEGER VALUES: I C METH COLLOCATION METHOD FOR APPROXIMATION. I C M # COLLOCATION PARAMETERS FOR APPROXIMATION. I C S # QUADRATURE POINTS FOR APPROXIMATION. I C L LOWER BOUND LOOP; L=2 IF METH -> LOBATTO, L=1 OTHERWISE. I C ORDER GLOBAL ORDER OF COLLOCATION METHOD FOR APPROXIMATION. I C ORDERQ ORDER OF QUADRATURE TO COMPUTE LAG TERMS FOR APPROXIMATION. I C METHR + I C MR I I C SR I AS METH, M, S, L, ORDER BUT FOR COLLOCATION METHOD TO I C LR I COMPUTE THE REFERENCE SOLUTION. I C ORDERR + I C ERRWGT ERROR WEIGHT INDICATOR (= IOPT(4)+1). I C NHFAIL TOTAL # STEPS FAILED. I C NERR LOGICAL UNIT # OF ERROR_MESSAGE_FILE (=CNTRL(2)). I C NWIR LOGICAL UNIT # OF FILE FOR INTERMEDIATE RESULTS (=CNTRL(3)). I C NSAV LOG.UN.# OF FILE FOR SAVING INFO IN CASE OF ERROR (=CNTRL(4)).I C MAXNC MAX. # SUBINTERVALS ALLOWED BY SIZE OF WORKING STORAGE. I C MAXKEV TOTAL # KERNEL EVALUATIONS ALLOWED. I C MAXCPS TOTAL # CPU-SECONDS ALLOWED. I C N # CURRENT SUBINTERVAL (TN,TN+HN]. I C OK_EXIT: # LAST INTERVAL + 1 . I C NCIT TOTAL # CORRECTOR ITERATIONS USED. I C NKEV TOTAL # KERNEL EVALUATIONS USED. I C NCPS # CPU-SECONDS USED WHEN "COLVI2" WAS CALLED. I C ON EXIT: # CPU-SECONDS USED BY "COLVI2". I C LOGICAL VALUES: I C VS IF TRUE, CONTROL ERROR, OTHERWISE USE FIXED STEPSIZE. I C GSSCKM IF TRUE, CHECK IF SOLUTION IS POLYNOMIAL OF DEGREE < M. I C ESCGSS IF TRUE, ESCAPE TO OTHER REF.SOL.APPROX. IN CASE SOL. IS POL. I C GEC IF TRUE, USE GLOBAL ERROR CONTROL. I C ULEC IF TRUE, PERFORM UNIFORM ERROR CONTROL. I C RLXTOL IF TRUE, RELAX TOLERANCE IN CASE OF FAILURE WITH MIN.STEPSIZE.I C GEETE IF TRUE, PROVIDE GLOBAL ERROR ESTIMATION IN ENDPOINT "TE". I C FUNCIT IF TRUE, USE FUNCTIONAL ITERATION AS CORRECTOR ITERATION I C WHILE SOLVING THE COLLOCATION SYSTEM. I C NEWTON IF TRUE, USE NEWTON'S METHOD; UPDATE JACOBIAN EACH ITERATION. I C IF FALSE AND IF FUNCIT IS FALSE: USE NEWTON'S METHOD; UPDATE I C JACOBIAN ONLY ONCE PER STEP. I C F.P. NUMBERS I C TOLLE REQUESTED TOLERANCE FOR GLOBAL ERROR (=REQTOL). I C TOLCIA TOLERANCE FOR ERROR IN CORRECTOR ITERATION PROCESS WHILE I C SOLVING COLLOCATION EQUATIONS TO COMPUTE THE APPROXIMATION. I C TOLCIR TOLERANCE FOR ERROR IN CORRECTOR ITERATION PROCESS WHILE I C SOLVING COLLOCATION EQUATIONS TO COMPUTE THE REFERENCE SOL. I C HMIN MIN. STEPSIZE ALLOWED. I C HMAX MAX. STEPSIZE ALLOWED. I C HC INTERVAL LENGTH FOR UNIFORM ERROR CONTROL (TE,(-HC),TN+HN). I C I C FOR A DESCRIPTION OF THE VARIABLES OF COMMON BLOCK "COLIXW" I C SEE BELOW UNDER "DISTRIBUTION WKAREA". I C I C I C RECURRING LOCAL VARIABLES: I C ------------------------- I C VAROUT INTERNAL FILE. I C USED TO CONVERT ARITHMETIC VALUES TO CHARACTER FORMAT. I C GAUSS LOGICAL. I C TRUE IF GAUSS COLLOCATION IS USED TO APPROXIMATE SOLUTION. I C LEC LOGICAL. I C TRUE IF LOCAL ERROR CONTROL IS USED. I C LOBAT LOGICAL. I C TRUE IF LOBATTO COLLOCATION IS USED TO APPROXIMATE SOLUTION. I C LOBATR LOGICAL. I C TRUE IF LOBATTO COLLOCATION IS USED TO COMPUTE REF.SOL. I C RSITCL LOGICAL. I C TRUE IF ITERATED COLLOCATION IS USED TO COMPUTE REF.SOL. I C C REAL C(M). I C CONTAINS THE COLLOCATION PARAMETERS. I C W REAL W(S). I C CONTAINS QUADRATURE WEIGHTS. I C LC REAL LC(M,L:M,L:S). I C LC(I,J,K) = L_I(CJ.CK). I C IDEM CR, WR AND LCR FOR COLLOCATION METHOD TO COMPUTE REF.SOL. I C LC1 REAL LC1(M). I C LC1(I) = L_I(1.0). I C LC0 REAL LC0(M). I C LC0(I) = L_I(0.0). I C LCG REAL LCG(0:M+1,M). I C LCG(I,J) = L_I(CJ/2) (BASED ON 0,CJ,1). I C LAGN REAL LAGN(NEQN). I C CONTAINS FC_N(TN). I C LAG REAL LAG(0:NEQN*(M-L+1)-1). I C CONTAINS FC_N(TN+CJ.HN). I C LAGNP1 REAL LAGNP1(NEQN). I C CONTAINS FC_N(TNP1). I C LAGSAV REAL LAGSAV(0:NEQN*(2*M-L+1)-1). I C CONTAINS FC_N(TNM1+CJ.HNM1) AND FC_N(TN+CJ.HN). I C GLAG REAL GLAG(0:NEQN*(M-L+1)-1). I C CONTAINS G(TN+CJ.HN) + FC_N(TN+CJ.HN). I C LEE REAL LEE(NEQN*((TE-T0)/HC+1)). I C CONTAINS LOCAL ERROR ESTIMATES OVER CURRENT INTERVAL IN I C TI = TE (-HC) T0 . I C LEESUM REAL LEESUM(NEQN*((TE-T0)/HC+1)). I C CONTAINS SUM OF LOCAL ERROR ESTIMATES OVER ALL INTERVALS UPTO I C THE CURRENT ONE IN TI = TE (-HC) T0 . I C NYPOLM REAL NYPOLM(NEQN). I C CONTAINS # CONSEC. TIMES A COMPONENT OF THE SOL. IS POLYN. I C H REAL H(0:MAXNC). I C CONTAINS SUBINTERVAL LENGTHS. I C H. REAL. I C "H." = H(.) . I C T. REAL. .-1 I C "T." = T0 + SUM H(I) . I C I=0 I C U REAL U(-NEQN:NEQN*(M-L+1)*MAXNC-1). I C CONTAINS APPROXIMATED SOLUTION IN COLLOCATION POINTS. I C UR REAL UR(-NEQN:NEQN*(MR-LR+1)*MAXNC-1). I C CONTAINS REFERENCE SOLUTION IN COLLOCATION POINTS. I C UN2 REAL UN2(0:NEQN*M-1). I C U(TN+CJ.HN/2). I C URN REAL URN(NEQN). I C UR(TN). I C URNP1 REAL URNP1(NEQN). I C UR(TN+HN). I C I C CONSTANTS USED: (OTHER THAN MACHINE CONSTANTS) I C -------------- I C GSSFAC = 2.0, FACTOR USED TO RELAX THE ORDER DEMAND IN THE CHECK ON I C POLYNOMIAL BEHAVIOR. LET UN2 BE THE SOLUTION COMPUTED I C OVER HALF THE N-TH INTERVAL THEN I C GEE(UN2)/GEE(UN) <= GSSFAC/(2**M) IMPLIES THAT THE I C SOLUTION IS NOT A POLYNOMIAL OF DEGREE < M. I C ("YPOLM"). I C HFAC = 0.9, REDUCTION FACTOR TO GET CONSERVATIVE GUESS OF THE I C STEPSIZE ("SOLVI2"). I C HFLFAC = 0.25, PENALTY REDUCTION FACTOR OF STEPSIZE. I C ("SOLVI2"). I C HMINFC = 1E-5, LIMIT FACTOR STEPSIZE: H(N)>=HINIT.HMINFC. I C ("CHKINI"). I C HRLFAC = 2.0, FACTOR TO LIMIT HNEW: 1/HRLFAC <= HNEW/HOLD <= HRLFAC.I C ("SOLVI2"). I C LSBITS = 128, 7 BITS FOR COMPUTATIONAL LOSS; USED IN TOLERANCE I C MATTERS ("CHKINI","ESCRGS"). I C MAXFIT = 15, MAX. NUMBER OF FUNCTIONAL ITERATIONS. I C ("SOLSYS"). I C MAXNIT = 10, MAX. NUMBER OF NEWTON ITERATIONS. I C ("SOLNEW"). I C NPGESC = 2, NUMBER OF CONSECUTIVE TIMES IT IS ALLOWED TO FIND I C POLYNOMIAL SOLUTION; IN CASE OF AUTOMATIC ESCAPE THE I C LAST NPGESC+1 STEPS ARE DISCARDED ("ESCRGS", "YPOLM").I C SAVFIL = 'COLSAV', FILE NAME ASSOCIATED WITH LOGICAL UNIT # CNTRL(4).I C ("SAVALL", "RELOAD"). I C TOLFRS = 0.1, FACTOR BY WHICH THE TOLERANCE FOR THE CORRECTOR ITER. I C PROCESS TO SOLVE THE COLLOC.EQ. FOR THE APPROX. IS I C MULTIPLIED TO GET THE TOLERANCE FOR THE CORR.IT.PROC. I C TO COMPUTE THE REF.SOL. ("CHKINI","SOLVI2","ESCRGS"). I C TOLMIN = LSBITS*SRELPR, MINIMUM TOLERANCE POSSIBLE (FOR "SRELPR" SEE I C BELOW) ("CHKINI","ESCRGS"). I C TOLMAX = 1.0, MAXIMUM VALUE TO WHICH TOLERANCE MAY BE RELAXED. I C ("SOLVI2"). I C TOLREL = 4.0, FACTOR TO RELAX THE TOLERANCE IN CASE IT CANNOT BE I C SATISFIED WITH MIN. STEPSIZE ("SOLVI2"). I C HDR = ' ERROR COLVI2...', ("CHKINI", "CHKFIL"). I C MCDEF = 81, ("CHKINI"). I C NCMI = 23, # VARIABLES IN COMMON /COLCMI/. ("SAVALL", "RELOAD").I C NCML = 9, # VARIABLES IN COMMON /COLCML/. ("SAVALL", "RELOAD").I C NCMR = 6, # VARIABLES IN COMMON /COLCMR/. ("SAVALL", "RELOAD").I C NCMIX = 21, # VARIABLES IN COMMON /COLIXW/. ("SAVALL", "RELOAD").I C I C I C DISTRIBUTION WKAREA: I C ------------------- I C NOTE: (...!...!:...!...!...) STANDS FOR: I C ---- IF ... THEN ... ELSE IF ... THEN ... ELSE ... ENDIF I C GAUSS = METH .EQ. 1 I C GAUSS1 = METH .EQ. 2 I C LEC = VS .A. .N.GEC I C LOBAT = METH .EQ. 3 I C LOBATR = METHR .EQ. 3 I C IG = METHR .EQ. 1 I C FUNCIT = TRUE IF FUNCTIONAL IT. IS USED TO SOLVE COLLOC. SYSTEM I C NEWTON = TRUE IF NEWTON'S METH. WITH JACOBIAN UPDATING IS USED I C I C L = ( LOBAT! 2! 1 ) ---- I C S = ( GAUSS1! M-1! M ) I I C ML = M-L+1 I IDEM LR,SR,MRL,SRL I C SL = S-L+1 ---- I C MW = MAX(MRL,ML).NEQN I C MB = M+ML I C NHC = (TE-T0)/HC I C NWKSYS = MW I C NIRVEC = ( FUNCIT! 2.NEQN! NEQN.(1+NEQN) ) I C I C 1 C(M) I C IC1= 1+M W(S) I C IC2=IC1+S LC(M) (C(L:M).C(L:S)) I C IC3=IC2+M.ML.SL (GAUSS! LC(M) (1.0)) I C ((VS.O.GEETE).A..N.IG I C IC4=IC3+(GAUSS!M) ! CR(MR) I C IC5=IC4+(VS.O.GEETE.A..N.IG!MR) WR(SR) I C IC6=IC5+(VS.O.GEETE.A..N.IG!SR) LCR(MR) (CR(LR:MR).CR(LR:SR)) I C ) I C (GSSCKM I C IC7=IC6+(VS.O.GEETE.A..N.IG ! LC(M) (0.0) I C !MR.MRL.SRL) I C IC8=IC7+(GSSCKM!M) LCG(0:M+1) (C(M)/2) I C ) I C I C ICE=IC8+(GSSCKM!(M+2).M) (ULEC!LEESUM(NEQN*(NHC+1)) I C IV1=ICE+(ULEC!NEQN.(NHC+1)) H(0:MAXNC) I C IV2=IV1+MAXNC+1 U(-NEQN:NEQN*ML*MAXNC-1) I C (.N.IG.A.(GEC.O.GEETE) I C IV3=IV2+NEQN.(1+ML.MAXNC) !UR(-NEQN:NEQN*MRL*MAXNC-1) I C !:LEC!URN(-NEQN:NEQN*MRL-1) I C ) I C (GSSCKM I C IVE=IV3+ !URN(NEQN) I C (.N.IG.A.(GEC.O.GEETE) I C !NEQN.(1+MRL.MAXNC) I C !:LEC!NEQN.(MRL+1)) I C IL1=IVE+(GSSCKM!NEQN) LAGN(NEQN) I C IL2=IL1+(GSSCKM!NEQN) LAGNP1(NEQN) I C ) I C IL3=IL2+(GSSCKM!NEQN) (IG!URNP1(NEQN)) I C NOTE: (.NOT.GSSCKM!URNP1->LAGNP1) I C ---- I C IL4=IL3+(IG!NEQN) (ULEC!LEE(NEQN*(NHC+1)) I C IL5=IL4+(ULEC!NEQN.(NHC+1)) (LEC!LAGSAV(0:NEQN*MB-1)) I C ILAG=IL5+(LEC LAG(0:NEQN*ML-1) I C !NEQN.M) NOTE: (LEC!LAG -> SECOND PART I C ---- OF LAGSAV) I C IL6=IL5+(LEC (GSSCKM! UN2(0:NEQN*M-1); I C !NEQN.MB) NOTE: (.N.LEC!UN2 -> LAG) I C ---- I C ILE=IL6+(GSSCKM!M.NEQN) (GSSCKM!NYPOLM(NEQN)) I C ITE=ILE+(GSSCKM!NEQN) GLAG(NEQN*ML),(G)LAGR(NEQN*MRL) I C AND (GSSCKM! GLAG(1/2)(NEQN*M) ) I C IW1=ITE+MW CORR(MW) AND (LEC! CB(MB) ) I C IW2=IW1+MW (.N.FUNCIT! DSYS(MW,MW)) I C (NEWTON I C IW3=IW2+(.N.FUNCIT!MW.MW) !WKAREA(MAX(NWKSYS,NIRVEC)) I C !:FUNCIT! WKAREA(NIRVEC) I C !WKAREA(NWKSYS+NIRVEC)) I C (FOR LIN.SYS.SOL. AND TEMP. STOR.)I C NOTE: (LEC! IW-IW1 SHOULD BE >= MB) I C ---- I C I C U(TI+CJ.HI)(K) -> U(NEQN*ML*I + NEQN*(J-L) + K-1) I C I=-1,J=M,K=1:NEQN; I=0:MAXNC-1,J=1:M,K=1:NEQN I C I C STORAGE OCCUPIED: I C ---------------- I C M+S+M.ML.SL+1+(1+NEQN.ML).MAXNC+NEQN+MW+MW + I C (GAUSS! M) + I C (VS.O.GEETE.A..N.IG! MR+SR+MR.MRL.SRL) + I C (GSSCKM! M+(M+2).M+NEQN.(4+M)) + I C (ULEC! 2.NEQN.(NHC+1)) + I C (.N.IG.A.(GEC.O.GEETE)! NEQN.(1+MRL.MAXNC)) + I C (LEC.A..N.GEETE! NEQN.(1+MRL)) + I C (IG! NEQN) + I C (LEC! NEQN.MB) + I C (FUNCIT.A.LEC! MAX(NIRVEC,MB-MW)) + I C (FUNCIT.A..N.LEC! NIRVEC) + I C (.N.FUNCIT! MW.MW) + I C (NEWTON! MAX(NWKSYS,NIRVEC) + I C (.N.FUNCIT.A..N.NEWTON! NWKSYS+NIRVEC)) I C I C -------------------------------------------------------------------- I C I C !!!!!!!! IMPORTANT !!!!!!!!! IMPORTANT !!!!!!!!! IMPORTANT !!!!!!!! !! C ========= ========= ========= !! C MACHINE CONSTANTS: !! C ----------------- !! C THE FOLLOWING MACHINE DEPENDENT CONSTANTS ARE USED IN THE PACKAGE: !! C !! C IBETA RADIX OF THE FLOATING-POINT REPRESENTATION. !! C IOVFLO LARGEST INTEGER VALUE "I" SUCH THAT ALL INTEGERS IN [-I,+I] !! C ARE REPRESENTABLE INTEGER NUMBERS. !! C NSDEC NUMBER OF SIGNIFICANT DECIMAL DIGITS. !! C IMXLUN LARGEST LOGICAL UNIT NUMBER ALLOWED BY THIS COMPILER. !! C DEFINED AS 999; ANSI STANDARD: IOVFLO. !! C SRELPR SMALLEST REAL VALUE "X" FOR WHICH 1.0-X < 1.0 < 1.0+X . !! C (STORED VALUES). !! C SOVFLO LARGEST REAL VALUE "X" SUCH THAT -X AND +X ARE !! C REPRESENTABLE F.P. NUMBERS. !! C SUNFLO SMALLEST REAL VALUE "X" SUCH THAT -X AND +X ARE !! C REPRESENTABLE F.P. NUMBERS. !! C !! C THEY ARE STORED IN THE COMMON BLOCKS /COLMCI/ AND /COLMCR/. THESE !! C COMMON BLOCKS ARE INITIALIZED BY THE ROUTINE "INICMC" (SEE BELOW). !! C !! C OTHER MACHINE DEPENDENCIES: !! C -------------------------- !! C !! C _ SINCE THE PACKAGE HAS BEEN DEVELOPED ON A MACHINE WITH A RATHER !! C LARGE WORD LENGTH, NO USE HAS BEEN MADE OF TYPE DOUBLE PRECISION !! C TO REPRESENT FLOATING POINT NUMBERS. IF DOUBLE PRECISION IS !! C REQUIRED CHANGE ALL TYPE REAL DECLARATIONS AND SPECIFICATIONS TO !! C TYPE DOUBLE PRECISION. NO UNDECLARED OR UNSPECIFIED VARIABLES ARE !! C USED. FOR F.P. INTRINSIC FUNCTIONS THE GENERIC NAMES ARE CHOSEN. !! C THE ROUTINES "DECLUF", "SOLLUF" AND "ZERPOL" SHOULD BE ADJUSTED !! C AND THE DOUBLE PRECISION VERSION OF THE "MACHAR" ROUTINE SHOULD !! C BE USED. !! C !! C _ TO SAVE THE NEEDED VARIABLES IN CASE OF A DETECTED ERROR AND IF !! C CNTRL(4)>0 AN EXTERNAL FILE WITH THE NAME "COLSAV" IS USED. !! C IF THIS IS NOT A LEGITIMATE FILE NAME CHANGE THE RELEVANT !! C PARAMETER STATEMENTS IN SUBROUTINES "SAVALL" AND "RELOAD". !! C !! C _ IF THE RECORD LENGTH (IN WORDS) OF A FILE OPENED FOR UNFORMATTED !! C I/O AND SEQUENTIAL ACCESS IS SMALLER THAN THE SIZE OF THE WORKING !! C STORAGE "WKAREA" + 13 ADDITIONAL WORDS, THEN THE ROUTINES "SAVALL"!! C AND "RELOAD" HAVE TO BE ADJUSTED. !! C !! C _ TO MONITOR THE CPU-TIME USED THE SUBROUTINES "COLVI2", "SOLVI2", !! C "SGEVI2" AND "CHKINI" INVOKE AN !! C INTEGER FUNCTION NCPJOB () !! C THAT SHOULD RETURN THE NUMBER OF CPU SECONDS USED SINCE A SPECIFIC!! C TIME, E.G. THE START OF THE JOB. !! C IN THE PACKAGE TWO VERSIONS OF "NCPJOB" ARE INCORPORATED. ONE FOR !! C A CDC CYBER-750 CALLING A REAL FUNCTION "SECOND". THE OTHER FOR A !! C VAX CALLING THE SYSTEM ROUTINE "ETIME". !! C !! C _ TO INITIALIZE THE MACHINE CONSTANTS "CHKINI" CALLS A: !! C SUBROUTINE INICMC !! C THAT INVOKES THE ROUTINE "MACHAR" OF W.J. CODY TO AUTOMATICALLY !! C DETERMINE SOME OF THE CONSTANTS. BOTH THE SINGLE AND THE DOUBLE !! C PRECISION VERSION OF "MACHAR" ARE INCORPORATED. !! C "IOVFLO" IS SET TO 2**31-1 (WHICH IS CORRECT FOR 32-BIT INTEGERS; !! C AND ACCEPTABLE FOR MACHINES WITH LARGER INTEGERS, SINCE IT IS ONLY!! C USED AS "A LARGE VALUE"). "IMXLUN" IS SET TO 999 WHICH IS ALSO !! C ACCEPTABLE FOR MOST OTHER MACHINES. !! C !! C _ TO DECOMPOSE THE JACOBIAN AND TO SOLVE THE SYSTEM OF LINEAR !! C EQUATIONS IN THE NEWTON PROCESS THE SUBROUTINES "DECLUF" AND !! C "SOLLUF" ARE CALLED BY "SOLNEW". THE HEADERS ARE: !! C SUBROUTINE DECLUF (A, N, IA, WKAREA, IERROR) !! C INTEGER N, IA, IERROR !! C REAL A(IA,*), WKAREA(*) !! C AND !! C SUBROUTINE SOLLUF (A, N, IA, B, WKAREA, IERROR) !! C INTEGER N, IA, IERROR !! C REAL A(IA,*), B(*), WKAREA(*) !! C THERE ARE TWO VERSIONS OF THESE ROUTINES. ONE THAT ASSUMES THE !! C AVAILABILITY OF AN IMSL LIBRARY ("DECLUF" CALLS "LUDATF" AND !! C "SOLLUF" CALLS "LUELMF"). THE OTHER MAKES USE OF A SIMPLE !! C (INCORPORATED) AX=B SOLVER (GAUSS ELIMINATION). !! C IF ANOTHER SYSTEM SOLVER NEEDS MORE WORK SPACE SOME STATEMENTS !! C HAVE TO BE ALTERED IN ROUTINES "CHKINI" AND "ESCRGS". !! C !! C _ TO DETERMINE THE HIGHER ORDER COLLOCATION PARAMETERS THE ROUTINE !! C "ZERPOL" IS CALLED: !! C SUBROUTINE ZERPOL (C, N, S, IERROR) !! C INTEGER N, IERROR !! C REAL C(0:N), S(N) !! C THAT SHOULD RETURN IN "S" THE SORTED REAL ZEROS OF THE POLYNOMIAL !! C C(0).Z**N + C(1).Z**(N-1) +...+ C(N-1).Z + C(N) = 0. !! C ONE VERSION OF "ZERPOL" CALLS THE SUBROUTINE "ZPOLR" OF THE IMSL !! C LIBRARY. THE SECOND VERSION INVOKES THE ACM-TOMS ROUTINE "RPOLY" !! C (ALSO INCORPORATED). !! C NB: DUE TO INTERNALLY DECLARED ARRAYS THE MAXIMUM DEGREE OF THE !! C POLYNOMIAL IS 100. !! C !! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! C I C -------------------------------------------------------------------- I C RETURN END SUBROUTINE COLVI2 + (NEQN, G,KC,DKCDY,LINEAR, T0,TE, REQTOL, + DEFOPT, IOPT,OPT,CNTRL, WKAREA,IW, TNC, UE, GEE, IERROR) C C ---------------------------------------------------------------------I C PURPOSE: SOLVE SYSTEM OF VIE2'S. FOR A DESCRIPTION OF "COLVI2" SEE I C ------- SUBROUTINE "COLDOC". I C I C COMMON VARIABLES: I C ---------------- I INTEGER METH, M, S, L, ORDER, ORDERQ, METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, N, NCIT, NKEV, NCPS COMMON /COLCMI/ METH, M, S, L, ORDER, ORDERQ, + METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, + NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, + N, NCIT, NKEV, NCPS * LOGICAL VS, GSSCKM, ESCGSS, GEC, ULEC, RLXTOL, GEETE, + FUNCIT, NEWTON COMMON /COLCML/ VS, GSSCKM, ESCGSS, GEC, ULEC, RLXTOL, GEETE, + FUNCIT, NEWTON * REAL TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC COMMON /COLCMR/ TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC * INTEGER IC1,IC2,IC3,IC4,IC5,IC6,IC7,IC8,ICE, IV1,IV2,IV3,IVE, + IL1,IL2,IL3,IL4,IL5,ILAG,IL6,ILE COMMON /COLIXW/ IC1,IC2,IC3,IC4,IC5,IC6,IC7,IC8,ICE, + IV1,IV2,IV3,IVE, IL1,IL2,IL3,IL4,IL5,ILAG,IL6,ILE * SAVE /COLCMI/, /COLCML/, /COLCMR/, /COLIXW/ C I C PARAMETER SPECIFICATION: I C ----------------------- I INTEGER NEQN, DEFOPT, IW, IERROR INTEGER IOPT(*), CNTRL(*) LOGICAL LINEAR REAL T0, TE, REQTOL, TNC REAL OPT(*), WKAREA(IW), UE(NEQN), GEE(NEQN) EXTERNAL G, KC, DKCDY C I C INVOKED BY: USER PROGRAM I C ---------- I C I C CHANGES IN COMMON VARIABLES: I C --------------------------- I C NCPS # CP SECONDS USED TO SOLVE VIE2 I C I C LOCAL VARIABLES: I C --------------- I CHARACTER*15 VAROUT LOGICAL FENTRY, GAUSS, LEC, RSITCL, ZLEESM REAL HINIT C FENTRY TRUE IF THIS IS THE FIRST CALL OF "COLVI2" I C ZLEESM TRUE IF "LEESUM" PART OF "WKAREA" HAS TO BE ZEROED; I.E. I C IF UNIFORM ERROR CONTROL IS USED AND NO PREVIOUSLY COMPUTED I C RELEVANT SUM OF LOCAL ERRORS IS AVAILABLE. I C HINIT INITIAL GUESS FOR LENGTH OF FIRST SUBINTERVAL. I C I C ---------------------------------------------------------------------I C INTEGER NCPJOB EXTERNAL NCPJOB * IERROR = 0 C C C CHECK PARAMETERS AND INITIALIZE COMMON; REACT ON RE-ENTRY C (RE)DISTRIBUTE WORKING STORAGE CALL CHKINI (NEQN, G, T0,TE, REQTOL, DEFOPT,IOPT,OPT,CNTRL, + WKAREA,IW, TNC, HINIT, ZLEESM, IERROR) C RETURN IF SOMETHING WAS WRONG IF (IERROR .NE. 0) RETURN C C FENTRY = N .EQ. 0 GAUSS = METH .EQ. 1 LEC = VS .AND. .NOT. GEC RSITCL = METHR .EQ. 1 C C C IF FIRST ENTRY, OR IF NEW OPTION VECTORS ARE GIVEN: C INITIALIZE COLLOCATION PARAMETERS, WEIGHT FACTORS AND C LAGRANGIAN INTERPOLATION COEFFICIENTS. IF (DEFOPT .EQ. 0) THEN C RE-ENTRY, OLD OPTION VECTIONS IF (CNTRL(1) .GT. 2) GOTO 40 C FIRST CALL OF COLVI2 IF (FENTRY) GOTO 10 C RE-ENTRY, NEW OPTION VECTORS IF (.NOT. GEC) GOTO 20 GOTO 30 ENDIF * C INITIALIZE COLL.PARS. FOR APPROX. OF SOL. 10 CALL COLCWL (METH, M, WKAREA(1),WKAREA(IC1),WKAREA(IC2), IERROR) IF (IERROR .NE. 0) RETURN C CALCULATE LAGR.INT.POL. COEFF. U(TN+1) IF (GAUSS) CALL COMPLV (1.0, M, WKAREA(1), WKAREA(IC3)) * 20 IF ((VS .OR. GEETE) .AND. .NOT. RSITCL) THEN C ERROR ESTIMATION REQUIRED, NO ITERATED COLLOCATION C INITIALIZE COLL.PARS. REF. SOL. CALL COLCWL (METHR, MR, WKAREA(IC4), WKAREA(IC5), WKAREA(IC6), + IERROR) IF (IERROR .NE. 0) RETURN ENDIF * 30 IF (GSSCKM) THEN C CALCULATE LAGR.INT. COEFF. U(TN) CALL COMPLV (0.0, M, WKAREA(1), WKAREA(IC7)) C CALCULATE LAGR.INT. COEFF. UR(TN) CALL COMPLG (M, WKAREA(1), WKAREA(IC8), WKAREA(ILE)) ENDIF C C C SOLVE VIE2 WITH CHOSEN COLLOCATION METHOD. C C STORE INITIAL STEPSIZE IN H(N) 40 WKAREA(IV1+N) = HINIT C C INITIALIZE, IF NEEDED, SOLUTION AND REF. SOL. VECTORS, AND C THE VECTOR CONTAINING THE SUM OF THE LOCAL ERRORS CALL INIVEC (FENTRY, + FENTRY .AND. (.NOT.RSITCL .AND. GEC), LEC, + GSSCKM, ZLEESM, NEQN, G, TNC, TE, WKAREA(IV2), + WKAREA(IV3), WKAREA(IVE), WKAREA(ICE)) C C SOLVE VIE2; IF EXIT OK: TNC=TE, UE=U(TE), C GEE=GLOBAL ERROR IN TE, IF (GEC) CALL SOLVI2 (NEQN, G,KC,DKCDY,LINEAR, T0,TE, + WKAREA(1),WKAREA(IC1),WKAREA(IC2), + WKAREA(IV1),WKAREA(IV2), + WKAREA(IC4),WKAREA(IC5),WKAREA(IC6),WKAREA(IV3), + WKAREA(ILAG), WKAREA(IL5), + WKAREA(ICE), ZLEESM, WKAREA(IL4), + WKAREA(IC3), WKAREA(IC7), WKAREA(IC8), + WKAREA(IL6), WKAREA(IL1), WKAREA(IL2), + WKAREA(IVE), WKAREA(IL3), WKAREA(ILE), + TNC, UE, GEE, IERROR) IF ((IERROR .EQ. 15) .AND. ESCGSS) THEN C SOLUTION IS POLYNOMIAL; ESCAPE TO HIGHER ORDER METHOD FOR REF.SOL. C AND LOCAL + UNIFORM ERROR CONTROL IERROR = 0 CALL ESCRGS (NEQN, WKAREA,IW, T0, TE, TNC, IERROR) IF (IERROR .NE. 0) GOTO 910 C C SOLVE VIE2; IF EXIT OK: TNC=TE, UE=U(TE) ZLEESM = .TRUE. CALL SOLVI2 (NEQN, G,KC,DKCDY,LINEAR, T0,TE, + WKAREA(1),WKAREA(IC1),WKAREA(IC2), + WKAREA(IV1),WKAREA(IV2), + WKAREA(IC4),WKAREA(IC5),WKAREA(IC6), WKAREA(IV3), + WKAREA(ILAG), WKAREA(IL5), + WKAREA(ICE), ZLEESM, WKAREA(IL4), + WKAREA(IC3), WKAREA(IC7), WKAREA(IC8), + WKAREA(IL6), WKAREA(IL1), WKAREA(IL2), + WKAREA(IVE), WKAREA(IL3), WKAREA(ILE), + TNC, UE, GEE, IERROR) ENDIF C C C CHECK IF ALL WENT WELL C IF (IERROR .NE. 0) GOTO 910 C C C IF GLOBAL ERROR CONTROL HAS BEEN USED "GEE" CONTAINS GLOBAL ERROR C IN "TE"; OTHERWISE C IF GLOBAL ERROR EST. IN TE IS REQUIRED PERFORM (ITERATED) COLL. TO C COMPUTE REFERENCE SOLUTION IN "TE"; IF NOT, STORE IN "GEE" C SUM OF LOCAL ERRORS IN "TE', IF AVAILABLE. IF (.NOT.GEC) THEN IF (GEETE) THEN C COMPUTE REFERENCE SOLUTION IN TE; ON EXIT: GEE=UR(TE) IF (RSITCL) THEN C ITERATED COLLOCATION, RESTORE OLD N VALUE N = N-1 CALL ITRCOL (TE, NEQN, G,KC, T0, WKAREA(IV1), WKAREA(1), + WKAREA(IC1),WKAREA(IV2),WKAREA(ILE),GEE,GEE) N = N+1 ELSE C SOLVE VIE2 WITH COL.PARS. FOR REF.SOL. USING THE SAME C STEPSIZES AS BEFORE C C STORE Y(T0) IN UR CALL COPYV (WKAREA(IV2), NEQN, WKAREA(IV3)) * CALL SGEVI2 (NEQN, G,KC,DKCDY,LINEAR, T0, + WKAREA(IC4), WKAREA(IC5), WKAREA(IC6), + WKAREA(IV1), WKAREA(IV3), + WKAREA(ILE), TNC, GEE, IERROR) IF (IERROR .NE. 0) GOTO 900 ENDIF C COMPUTE GLOBAL ERROR IN TE CALL ADDABV (GEE, NEQN, -1.0, UE) ELSE IF (ULEC) THEN C STORE SUM OF LOCAL ERRORS IN "TE", LEESUM(0:NEQN-1), IN GEE CALL COPYV (WKAREA(ICE), NEQN, GEE) ENDIF ENDIF C C C COMPUTE NUMBER OF CP-SECONDS USED NCPS = NCPJOB() - NCPS RETURN C C ALAS! 900 CONTINUE CALL ERRMSG ('PROBLEMS WITH COMPUTATION OF GLOBAL ERROR IN "TE"') 910 CONTINUE WRITE(VAROUT,'(E15.5)') TNC CALL ERRMSG ('ENDPOINT NOT REACHED, LAST T-VALUE :'//VAROUT) C C IF REQUIRED SAVE VARIABLES ON FILE FOR RE-ENTRY IF (NSAV .NE. 0) CALL SAVALL (WKAREA,IW, DEFOPT,IOPT,OPT, TE, TNC) * C C COMPUTE NUMBER OF CP-SECONDS USED NCPS = NCPJOB() - NCPS RETURN END SUBROUTINE SOLVI2 + (NEQN, G,KC,DKCDY,LINEAR, T0,TE, C,W,LC, H, U, CR,WR,LCR, UR, + LAG, LAGSAV, LEESUM,ESTGEE,LEE, LC1, LC0,LCG,UN2,LAGN,LAGNP1, + URN, URNP1, WKAREA, TN, UN, GEE, IERROR) C C ---------------------------------------------------------------------I C PURPOSE: SUPERVISE PROCESS OF SOLVING VIE2; PERFORM ERROR CONTROL, I C ------- MONITOR STEPS, CHECK ON POLYNOMIAL SOLUTION, IF REQUIRED. I C I C COMMON VARIABLES: I C ---------------- I INTEGER METH, M, S, L, ORDER, ORDERQ, METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, N, NCIT, NKEV, NCPS COMMON /COLCMI/ METH, M, S, L, ORDER, ORDERQ, + METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, + NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, + N, NCIT, NKEV, NCPS * LOGICAL VS, GSSCKM, ESCGSS, GEC, ULEC, RLXTOL, GEETE, + FUNCIT, NEWTON COMMON /COLCML/ VS, GSSCKM, ESCGSS, GEC, ULEC, RLXTOL, GEETE, + FUNCIT, NEWTON * REAL TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC COMMON /COLCMR/ TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC * REAL SRELPR, SOVFLO, SUNFLO COMMON /COLMCR/ SRELPR, SOVFLO, SUNFLO * SAVE /COLCMI/, /COLCML/, /COLCMR/, /COLMCR/ C I C PARAMETER SPECIFICATION: I C ----------------------- I INTEGER NEQN, IERROR LOGICAL LINEAR, ESTGEE REAL T0, TE, TN REAL C(M), W(S), LC(M,L:M,L:S), + H(0:MAXNC), U(-NEQN:NEQN*(M-L+1)*MAXNC-1), + CR(*),WR(*),LCR(MR,LR:MR,LR:*), UR(-NEQN:*), + LAG(0:NEQN*(M-L+1)-1), LAGSAV(0:*), LEESUM(*), LEE(*), + LC1(*),LC0(*), LCG(0:M+1,*), UN2(0:*), LAGN(*), LAGNP1(*), + URN(*), URNP1(*), + WKAREA(*), UN(NEQN), GEE(NEQN) EXTERNAL G, KC, DKCDY C C - COLL. VARS NEEDED I C W I-> TO APPROXIMATE I C LC - SOLUTION. I C H ENTRY: H(0:N-1) SUBINT. LENGTHS OF ALL PREVIOUS STEPS TAKEN I C H(N) INITIAL GUESS FOR STEPSIZE IN N-TH INTERVAL I C EXIT: H(0:NC-1) SUBINT. LENGTHS OF ALL STEPS TAKEN (TNC = TN)I C H(NC) INITIAL GUESS FOR LENGTH NEXT SUBINTERVAL I C U ENTRY: SHOULD CONTAIN U(TI+CJ.HI) I=0:N-1,J=1:M I C NOTE: U(-NEQN:-1) SHOULD CONTAIN Y(T0) I C EXIT: CONTAINS U(TI+CJ.HI) I=0:NC-1,J=1:M I C U(TIJ) -> U(NEQN*((M-L+1)*I+(J-L))+(0:NEQN-1)) I C CR - COLL. VARS FOR REF. SOL (NOT NEEDED IF CONSTANT STEPSIZES I C WR I-> WILL BE USED OR IF THE REF. SOL. IS COMPUTED WITH I C LCR - ITERATED COLLOCATION. I C UR IF GLOBAL ERROR CONTROL IS REQUIRED AND THE REF.SOL. METHOD I C IS NOT ITERATED COLLOCATION: I C ENTRY: SHOULD CONTAIN UR(TI+CRJ.HI) I=0:N-1,J=1:M I C NOTE: UR(-NEQN:-1) SHOULD CONTAIN Y(T0) I C EXIT: CONTAINS UR(TI+CRJ.HI) I=0:NC-1,J=1:MR I C UR(TIJ) -> UR(NEQN*((MR-LR+1)*I+(J-LR))+(0:NEQN-1)) I C IF LOCAL ERROR CONTROL HAS BEEN SELECTED: I C ENTRY: UR(-NEQN:-1) SHOULD CONTAIN UR(TN) I C EXIT: UR(-NEQN:-1) CONTAINS UR(TNC) I C UR(0:NEQN*(MR-LR+1)-1) CONT. UR(T(NC-1)+CRJ.H(NC-1))I C J=LR:MR I C LAG WORKING STORAGE FOR LAGTERM FCN(TN+CJ.HN) J=L:M I C FCN(TNJ) -> LAG(NEQN*(J-L)+(0:NEQN-1)) I C LAGSAV IN CASE OF LOCAL ERROR CONTROL: I C WORKING STORAGE FOR LAGTERM FCN(TN-H(N-1)+CJ.H(N-1)) AND I C FCN(TN+CJ.HN) J=1:M I C FCN(TNM1+CJ.HNM1) -> LAGSAV(NEQN*(J-1)+(0:NEQN-1)) J=1:M I C FCN(TNJ) -> LAGSAV(NEQN*(M+J-L)+(0:NEQN-1)) J=L:M I C NOTE: IN THIS CASE LAG SHARES THE STORAGE LOCATIONS WITH I C ---- LAGSAV(NEQN*M:NEQN*(2*M-L+1)-1) I C LEESUM IF UNIFORM ERROR CONTROL HAS BEEN SELECTED: I C ENTRY: LEESUM(I*NEQN+(1:NEQN)) SHOULD CONTAIN EITHER 0.0 OR I C (K=0,N-1) SUM LEE_K(I*NEQN+(1:NEQN)) I=0,... I C (LOCAL ERRORS OVER (TK,TK+1) IN TI) T=TI=TE,(-HC),TN I C (CF. CHKINI) I C EXIT: APPROX.OF (K=0,NC-1) SUM LEE_K(I*NEQN+(1:NEQN)) I=0,...I C ESTGEE IF UNIFORM ERROR CONTROL HAS BEEN SELECTED: I C ENTRY: TRUE IF NO ERROR ESTIMATES OVER [T0,TN] ARE AVAILABLE I C LEE IF UNIFORM ERROR CONTROL: I C WORKING STORAGE FOR LOCAL ERROR EST. IN T=TI FOR I=0,... I C LC1 IF APPROX. METH. IS GAUSS: I C COLL. VAR. FOR COMPUT. OF SOL. IN STEPPOINTS I C -I IF CHECK ON POL. SOL. IS REQUIRED: I C LC0 I LAGR.COEFF.TO COMPUTE SOL. IN TN I C LCG I LAGR.COEFF.TO COMPUTE FCN(TN+CJ.HN/2) I C UN2 I WORKING STORAGE FOR U(TN+CJ.HN/2) (0:NEQN*M-1) I C LAGN I WORKING STORAGE FOR FCN(TN) I C LAGNP1 I STORAGE FOR FCN(TNP1) I C URN -I ENTRY: UR(TN) I C URNP1 IF REF.SOL. IS COMPUTED WITH ITERATED COLLOCATION: I C WORKING STORAGE FOR UR(TN+HN) I C WKAREA WORKING STORAGE FOR INTERMEDIATE VECTOR RESULTS AND I C FOR SOLVING THE COLL. SYSTEM (CF. SOLSYS) I C TN ENTRY: LEFT ENDPOINT OF N-TH SUBINTERVAL I C EXIT: RIGHT ENDPOINT OF LAST SUBINTERVAL I C UN EXIT: APPROX.SOL. IN "TN" I C GEE EXIT: IN CASE OF GLOBAL ERROR CONTROL: I C GLOBAL ERROR ESTIMATE IN "TN" I C IERROR ENTRY: 0 I C EXIT: 0: OK I C 11: FAILURE WITH MINIMUM STEPSIZE I C 12: SIZE WORKING STORAGE AREA TOO SMALL I C 13: # KERNEL EVAL. > MAX. # KERN. EV. ALLOWED I C 14: # CPU-SECONDS > MAX. # CPU-SEC. ALLOWED I C 15: POLYNOMIAL SOLUTION I C 16: TOLERANCE WOULD BE RELAXED TO A VALUE > 1.0 I C OTHER: ERRORS FROM "SOLSYS" I C I C INVOKED BY: COLVI2 I C ---------- I C I C CHANGES IN COMMON VARIABLES: I C --------------------------- I C NHFAIL ADDED 1 AFTER EACH FAILED STEP I C N ADDED 1 AFTER EACH SUCCESSFUL STEP I C TOLLE + INCREASED BY A FACTOR TOLREL, IF STEP I C TOLCIA I FAILED WITH MINIMUM STEPSIZE I C TOLCIR + I C I C CONSTANTS: I C --------- I REAL HFAC, HFLFAC, HRLFAC, TOLFRS, TOLMAX, TOLREL PARAMETER (HFAC = 0.9) PARAMETER (HFLFAC = 0.25) PARAMETER (HRLFAC = 2.0) PARAMETER (TOLFRS = 0.1) PARAMETER (TOLMAX = 1.0) PARAMETER (TOLREL = 4.0) C I C LOCAL VARIABLES: I C --------------- I CHARACTER*10 VAROUT INTEGER IHFAIL, INDEXN, INDXRN, INDRN1, ITE, IW1,IW2,IW3, MW, NC LOGICAL GAUSS, LAST, LEC, LOBAT, LOBATR, RSITCL, WRINT REAL GELIML,GELIMU, HMINN, HNM1,HN, INVP,INVQ, + LELIML,LELIMU, WGEE, WLEE, WULEE C IHFAIL # OF STEPS FAILED IN THE CURRENT SUBINTERVAL I C INDEXN POINTER TO SOLUTION IN 1-ST COLLOC.POINT IN N-TH INTERVAL I C INDXRN POINTER TO REF.SOL. IN 1-ST COLLOC.POINT IN N-TH INTERVAL I C INDRN1 POINTER TO REF.SOL. IN TN+HN I C ITE, ..., IW3 POINTERS TO WKAREA LOCATIONS (CF. "COLDOC" SUB I C "DISTRIBUTION WKAREA" I C MW MAX. DIMENSION COLLOCATION SYSTEM TO BE SOLVED I C NC # STEPS TO BE TAKEN IN CASE OF CONSTANT STEPSIZES I C LAST TRUE, IF TN+HN = TE I C WRINT TRUE, IF INTERMEDIATE RESULTS ARE REQUIRED I C GELIML - LIMITS TO ERROR TERM IN GLOBAL ERROR CONTROLLED STEPSIZE I C GELIMU / STRATEGY, SO THAT RESULTING FACTOR LIES BETWEEN 0.5 AND 2.0 I C HMINN MIN. STEPSIZE FOR CURRENT SUBINTERVAL I C HNM1 LENGTH PREVIOUS SUBINTERVAL I C HN GUESS FOR LENGTH CURRENT SUBINTERVAL, STORED IN H(N) I C INVP 1/ORDER - POWERS USED IN I C INVQ 1/ORDERQ / STEPSIZE STRATEGY I C LELIML - LIMITS TO ERROR TERM IN LOCAL ERROR CONTROLLED STEPSIZE I C LELIMU / STRATEGY, SO THAT RESULTING FACTOR LIES BETWEEN 0.5 AND 2.0 I C WGEE WEIGHTED NORM OF GLOBAL ERROR IN TN+HN I C WLEE WEIGHTED NORM OF LOCAL ERROR IN TN+HN I C WULEE MAX. OF WEIGHTED NORMS OF ERRORS OVER [TN+HN,TE] I C I C ---------------------------------------------------------------------I C INTEGER J, ML, MRL, MWR, MWS, NCPJOB, NL2, NLR2 LOGICAL ACCEPT, RSFAIL, POLY, YPOLM REAL FACH, LEEWGT, R, UEEWGT, WMXNRM EXTERNAL LEEWGT, NCPJOB, UEEWGT, WMXNRM, YPOLM * ML = M-L+1 MRL = MR-LR+1 C C C DISTRIBUTE WORKING STORAGE "WKAREA" FOR "YPOLM", C LAG TERMS AND LINEAR SYSTEM SOLVER MWS = ML*NEQN MWR = MRL*NEQN MW = MAX(MWS, MWR) IF (GSSCKM) THEN ITE = 1+NEQN ELSE ITE = 1 ENDIF IW1 = ITE + MW IW2 = IW1 + MW IF (FUNCIT) THEN IW3 = IW2 ELSE IW3 = IW2 + MW*MW ENDIF C C C INITIALIZE LOOP CONSTANTS AND VARIABLES GAUSS = METH .EQ. 1 LAST = .FALSE. LEC = VS .AND. .NOT.GEC LOBAT = L .EQ. 2 LOBATR = LR .EQ. 2 RSITCL = METHR .EQ. 1 WRINT = (NWIR .NE. 0) * IHFAIL = 0 IF (.NOT. VS) THEN R = (TE-TN)/H(N) NC = INT(R) IF (R-NC .GT. 0.0) NC = NC+1 NC = NC + N ENDIF INDXRN = NEQN*(1-LR) NL2 = NEQN*(L -2) NLR2 = NEQN*(LR-2) * GELIML = TOLLE / (HRLFAC**ORDER) GELIMU = TOLLE * (HRLFAC**ORDER) HN = H(N) IF (LEC) THEN IF (N .EQ. 0) THEN HNM1 = HN ELSE HNM1 = H(N-1) ENDIF ENDIF INVP = 1.0 / ORDER INVQ = 1.0 / ORDERQ LELIML = TOLLE / (HRLFAC**ORDERQ) LELIMU = TOLLE * (HRLFAC**ORDERQ) * IF (GSSCKM) THEN C ZERO NYPOLM (# CONSEC. TIMES SOLUTION BEHAVES AS A POLYNOMIAL) CALL ZEROV (WKAREA(1), NEQN) C STORE FCN(TN) IN LAGN CALL INILGN (TN, NEQN,G, U, LC1, WKAREA(IW3), LAGN) ENDIF * IF (LEC) C IN CASE OF LOC. ERR. CONTR. COMPUTE FCN(T(N-1)+CJ.HN-1) J=1,...,M C AND STORE IN LAGSAV + CALL INILAG (TN, NEQN,G,KC, C,W,LC,LOBAT, H, U, + WKAREA(IW3), LAGSAV) GOTO 100 C C C LOOP ENTRY IN CASE A FAILURE WITH MINIMUM STEPSIZE OCCURRED. 90 IF (.NOT.RLXTOL .OR. N .EQ. 0) GOTO 920 * C RELAX TOLERANCE; USE H(N-1) AS INITIAL GUESS FOR STEPSIZE TOLLE = TOLLE * TOLREL IF (TOLLE .GT. TOLMAX) C TOLERANCE RELAXED UP TO 1; NO USE TO GO ON + GOTO 930 * WRITE(VAROUT,'(E10.3)') TN CALL ERRMSG ('CANNOT MEET REQTOL, FAILURE WITH MIN. H IN T='// + VAROUT) WRITE(VAROUT,'(E10.3)') TOLLE CALL ERRMSG (' REQTOL RELAXED TO '//VAROUT) IF (TOLCIR/TOLCIA .LE. TOLFRS) TOLCIR = TOLCIR * TOLREL TOLCIA = TOLCIA * TOLREL GELIML = TOLLE / (HRLFAC**ORDER) GELIMU = TOLLE * (HRLFAC**ORDER) LELIML = TOLLE / (HRLFAC**ORDERQ) LELIMU = TOLLE * (HRLFAC**ORDERQ) HN = H(N-1) C C C LOOP UNTIL "TE" IS REACHED 100 CONTINUE C C LIMIT HN C COMPUTE MIN. STEPSIZE FOR THIS INTERVAL HMINN = MAX(HMIN,SRELPR*ABS(TN)) HN = MAX(HN,HMINN) HN = MIN(HN,HMAX) IF (TE-(TN+HN) .LE. ABS(TE)*SRELPR) THEN HN = TE-TN LAST = .TRUE. ELSE IF (TE-(TN+1.5*HN) .LE. ABS(TE)*SRELPR) THEN C AVOID A VERY SMALL LAST STEPSIZE (ERROR ESTIMATION PROBLEMS) HN = (TE-TN)/1.5 ENDIF C C C ENTRY LOOP IN CASE OF CONSTANT STEPSIZES 110 CONTINUE IERROR = 0 INDEXN = NEQN*(1-L+ML*N) C C STORE INITIAL GUESS OF HN H(N) = HN C C CHECK # KEV IF (NKEV .GT. MAXKEV) GOTO 900 C C CHECK CPU-TIME IF (NCPJOB()-NCPS .GT. MAXCPS) GOTO 910 C C CCCCCCCCCC SOLVE COLLOCATION EQUATIONS ON [TN,TN+HN] C C INITIALIZE U(TN+CJ.HN) ON U(TN) DO 200 J = L, M CALL COPYV (U(INDEXN+NL2), NEQN, U(INDEXN+NEQN*(J-1))) 200 CONTINUE * CALL SLQCE2 (TN,H, NEQN, G,KC,DKCDY,LINEAR, T0, + C,W,LC, M,S,L, LOBAT, TOLCIA, .FALSE., + LAG, WKAREA(ITE),WKAREA(IW1),WKAREA(IW2),MWS, + WKAREA(IW3), U, LAGNP1, IERROR) C C CHECK IF ALL WENT WELL IF (IERROR .NE. 0) THEN C CONTINUE ONLY IF CONVERGENCE PROBLEMS IN VARIABLE STEP SIZE CASE IF (.NOT.(VS .AND. IERROR .EQ. 21)) GOTO 980 C TRY AGAIN IHFAIL = IHFAIL + 1 IF (HN .LE. HMINN) THEN C IF FAILED WITH MINIMUM STEPSIZE RELAX TOLERANCE (IF ALLOWED) GOTO 90 ELSE C IF NOT, DECREASE STEPSIZE HN = HN*HFLFAC ENDIF GOTO 100 ENDIF C C C OK; STORE SOLUTION IN TN+HN IN UN IF (GAUSS) THEN C COMPUTE COLLOCATION SOLUTION IN TN+HN WITH LAGR. INTERP. CALL UTILIP (NEQN, U(INDEXN), LC1, UN) ELSE CALL COPYV (U(INDEXN+NEQN*(M-1)), NEQN, UN) ENDIF * * IF (.NOT. VS) THEN C C CCCCCCCCCC CONSTANT STEPSIZES; ADJUST LOOPVARIABLES C IF (WRINT) CALL WRIRES (TN, HN, UN, UN, NEQN, WKAREA(ITE)) C TAKE NEXT STEP TN = TN+HN N = N+1 IF (N .LT. NC) GOTO 110 ELSE C C CCCCCCCCCC VARIABLE STEPSIZES, ERROR CONTROL C CCCCCCCCCC COMPUTE REFERENCE SOLUTION C RSFAIL = .FALSE. * IF (LEC) THEN C LOCAL ERROR ESTIMATION C COMPUTE LAGTERM BY INTERPOLATION C INITIALIZE UR(TN+CRJ.HN) ON UR(TN) DO 210 J = LR, MR CALL COPYV (UR(-NEQN), NEQN, UR(NEQN*(J-LR))) 210 CONTINUE * 220 CALL SLICE2 (TN, HNM1,HN, NEQN, G,KC,DKCDY,LINEAR, + C, WKAREA(IW1), CR, WR, LCR, LOBATR, + LAGSAV, WKAREA(ITE), WKAREA(IW1), + WKAREA(IW2),MWR,WKAREA(IW3), UR(INDXRN),IERROR) IF (IERROR .EQ. 21 .AND. .NOT. RSFAIL) THEN C CONVERGENCE PROBLEMS, TRY AGAIN WITH BETTER ESTIMATE TO START WITH RSFAIL = .TRUE. IERROR = 0 GOTO 220 ENDIF ELSE IF (RSITCL) THEN C GLOBAL ERROR ESTIMATION C COMPUTE REF. SOL. IN TN+HN WITH ITER. COLL. CALL ITRCOL (TN+HN, NEQN, G,KC, T0, H, C,W, U, WKAREA(IW3), + LAGNP1, URNP1) ELSE C COMPUTE REF. SOL. WITH HIGHER ORDER COLLOCATION METHOD C COMPUTE LAG TERM WITH QUADRATURE INDXRN = NEQN*(1-LR+MRL*N) C INITIALIZE UR(TN+CRJ.HN) ON UR(TN) DO 230 J = LR, MR CALL COPYV (UR(INDXRN+NLR2),NEQN, UR(INDXRN+NEQN*(J-1))) 230 CONTINUE * 240 CALL SLQCE2 (TN,H, NEQN, G,KC,DKCDY,LINEAR, T0, CR,WR,LCR, + MR,SR,LR,LOBATR, TOLCIR, GSSCKM, WKAREA(ITE), + WKAREA(ITE),WKAREA(IW1),WKAREA(IW2),MWR, + WKAREA(IW3), UR, LAGNP1, IERROR) IF (IERROR .EQ. 21 .AND. .NOT. RSFAIL) THEN C CONVERGENCE PROBLEMS, TRY AGAIN WITH BETTER ESTIMATE TO START WITH RSFAIL = .TRUE. IERROR = 0 GOTO 240 ENDIF ENDIF C CHECK IF ALL WENT WELL IF (IERROR .NE. 0) GOTO 940 * INDRN1 = INDXRN+NEQN*(MR-1) C C IF REQUIRED, GIVE INFO IF (WRINT) THEN IF (RSITCL) THEN CALL WRIRES (TN, HN, URNP1, UN, NEQN, WKAREA(ITE)) ELSE CALL WRIRES (TN, HN, UR(INDRN1), UN, NEQN, WKAREA(ITE)) ENDIF ENDIF C C CCCCCCCCCC CONTROL ERROR C IF (GEC) THEN C CONTROL GLOBAL ERROR; GEE = UR(TN+HN)-U(TN+HN) IF (RSITCL) THEN CALL ADDV (GEE, NEQN, 1.0, URNP1, -1.0, UN) WGEE = WMXNRM (GEE, URNP1, NEQN) ELSE CALL ADDV (GEE, NEQN, 1.0, UR(INDRN1), -1.0, UN) WGEE = WMXNRM (GEE, UR(INDRN1), NEQN) ENDIF C ACCEPT STEP IF GLOBAL ERROR <= TOLLE AND C IN CASE OF UNIFORM ERROR CONTROL IF ALSO UNIFORM ERROR <= TOLLE ACCEPT = WGEE .LE. TOLLE IF (ULEC) THEN WULEE = UEEWGT (TN,HN, NEQN,KC, T0,TE, C,W, + CR,WR, U(INDEXN), UR(INDXRN), + LEESUM, ESTGEE, WKAREA(IW3), LEE) ACCEPT = ACCEPT .AND. WULEE .LE. TOLLE ENDIF ELSE C CONTROL LOCAL ERROR; LEE = (TN,TN+HN) INT K(TN+HN,.,.) - SUM ... WLEE = LEEWGT (TN,HN, NEQN,KC, T0, C,W, CR,WR, + U(INDEXN), UR(INDXRN), WKAREA(IW3)) IF (ULEC) THEN C CONTROL ERROR OVER WHOLE INTERVAL WULEE = UEEWGT (TN,HN, NEQN,KC, T0,TE, C,W, + CR,WR, U(INDEXN), UR(INDXRN), + LEESUM, ESTGEE, WKAREA(IW3), LEE) WLEE = MAX(WLEE, WULEE) ENDIF ACCEPT = WLEE .LE. TOLLE ENDIF * IF (ACCEPT) THEN C IF (GSSCKM) THEN C CHECK IF SOLUTION BEHAVES AS A POLYNOMIAL OF DEGREE < M RSFAIL = .FALSE. 250 POLY = YPOLM (TN,HN, NEQN,G,KC,DKCDY,LINEAR, C,W,LC, + LCG,LC0, LAGN,LAG,LAGNP1, UN2, + WKAREA(1), WKAREA(ITE), WKAREA(IW1), + WKAREA(IW2),MWS, WKAREA(IW3), + U(INDEXN), URN, IERROR) IF (IERROR .EQ. 21 .AND. .NOT. RSFAIL) THEN C CONVERGENCE PROBLEMS, TRY AGAIN WITH BETTER ESTIMATE TO START WITH RSFAIL = .TRUE. IERROR = 0 GOTO 250 ENDIF IF (IERROR .NE. 0) GOTO 950 * IF (POLY) C SOL. VIE2 FOUND TO BE POLYNOMIAL IN TN "NPGESC" C CONSECUTIVE TIMES; IT IS ASSUMED THAT THE SOL. IS A POL. OF C DEGREE < M; ESCAPE FROM GAUSS+REFSOL METHOD + GOTO 970 ENDIF C CCCCC STEP ACCEPTED, ADJUST (LOOP) VARIABLES C CHECK RESERVE WORKING STORAGE IF (N+1 .GE. MAXNC) GOTO 960 C NHFAIL = NHFAIL + IHFAIL IHFAIL = 0 IF (GSSCKM) THEN C COPY UR(TN+HN) INTO URN IF (RSITCL) THEN CALL COPYV (URNP1, NEQN, URN) ELSE CALL COPYV (UR(INDRN1), NEQN, URN) ENDIF C STORE FCN(TN+HN) IN LAGN CALL G(TN+HN, WKAREA(IW3)) CALL ADDV (LAGN, NEQN, 1.0,URN, -1.0,WKAREA(IW3)) ENDIF IF (LEC) THEN C ADJUST LAGSAV ARRAY FOR NEXT STEP CALL ADJLSV (TN,HN, NEQN,KC, C,W, U(INDEXN), + LAGSAV) HNM1 = HN C COPY UR(TN+HN) TO UR(-NEQN:) CALL COPYV (UR(INDRN1), NEQN, UR(-NEQN)) ENDIF IF (ULEC) THEN C ADD LEE(I) TO SUM OF LOCAL ERRORS IN TI, TI=TE,-HC,TN+HN CALL ADDABV (LEESUM, NEQN*INT((TE-TN-HN)/HC+1),1.0, LEE) ESTGEE = .FALSE. ENDIF * TN = TN + HN N = N+1 ELSE C CCCCC STEP REJECTED IHFAIL = IHFAIL + 1 C CHECK IF FAILED WITH MIN. STEPSIZE IF (HN .LE. HMINN) THEN C IF SO, RELAX TOLERANCE IF ALLOWED; GOTO 90 ELSE IF (MOD(IHFAIL,2) .EQ. 0) THEN C IF FAILED REPEATEDLY, DECREASE THE STEPSIZE WITH EXTRA FACTOR HN = HN*HFLFAC GOTO 100 ENDIF LAST = .FALSE. ENDIF C CCCCCCCCCC COMPUTE NEW STEPSIZE C IF (GEC) THEN R = MAX(GELIML,WGEE) R = MIN(GELIMU,R) FACH = (TOLLE/R) ** INVP IF (ULEC) THEN R = MAX(LELIML,WULEE) R = MIN(LELIMU,R) FACH = MIN(FACH, (TOLLE/R) ** INVQ) ENDIF HN = HN*HFAC*FACH ELSE R = MAX(LELIML,WLEE) R = MIN(LELIMU,R) FACH = (TOLLE/R) ** INVQ HN = HN*HFAC*FACH ENDIF C C LOOP IF "TE" NOT REACHED IF (.NOT. LAST) GOTO 100 C C FINISHED C STORE GUESS FOR LENGTH NEXT INTERVAL H(N) = HN ENDIF RETURN C C C ERROR RETURNS C C NUMBER OF KERNEL EVALUATIONS TOO LARGE 900 CALL ERRMSG ('NUMBER OF KERNEL EVALUATIONS EXCEEDS IOPT(6)') IERROR = 13 RETURN C C TOO MUCH CPU-TIME USED 910 CALL ERRMSG ('CPU-TIME USED EXCEEDS IOPT(7)') IERROR = 14 RETURN C C FAILED TO MEET TOLERANCE WITH MINIMUM STEPSIZE 920 NHFAIL = NHFAIL + IHFAIL CALL ERRMSG ('CANNOT MEET REQTOL, FAILURE WITH MIN. STEPSIZE') IERROR = 11 RETURN C C TOLERANCE INCREASED UNACCEPTABLY 930 NHFAIL = NHFAIL + IHFAIL CALL ERRMSG ('RELAXATION WOULD RESULT IN A TOLERANCE > 1.0') IERROR = 16 RETURN C C COMPUTATION REFERENCE SOLUTION FAILED 940 NHFAIL = NHFAIL + IHFAIL CALL ERRMSG (' ERROR OCCURRED WHILE COMPUTING REF.SOL.'// + 'TO APPROXIMATE ERROR') RETURN C C CHECK ON POLYNOMIAL SOLUTION FAILED 950 NHFAIL = NHFAIL + IHFAIL CALL ERRMSG (' ERROR OCCURRED WHILE CHECKING WHETHER'// + ' SOLUTION IS POLYNOMIAL') RETURN C C SIZE WORKING STORAGE TOO SMALL 960 WRITE(VAROUT,'(I10)') N+2 CALL ERRMSG ('SIZE WORKING STORAGE TOO SMALL FOR'//VAROUT// + ' SUBINTERVALS') IERROR = 12 RETURN C C SOLUTION POLYNOMIAL; PROBLEMS IN CASE OF GAUSS COLL. PARS. 970 IERROR = 15 980 RETURN * END SUBROUTINE SGEVI2 + (NEQN, G,KC,DKCDY,LINEAR, T0, CR,WR,LCR, H, UR, WKAREA, + TN, URN, IERROR) C C ---------------------------------------------------------------------I C PURPOSE: COMPUTE REFERENCE SOLUTION BY HIGHER ORDER COLLOCATION I C ------- METHOD. INVOKED TO ESTIMATE THE GLOBAL ERROR IN TN+H(N). I C I C COMMON VARIABLES: I C ---------------- I INTEGER METH, M, S, L, ORDER, ORDERQ, METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, N, NCIT, NKEV, NCPS COMMON /COLCMI/ METH, M, S, L, ORDER, ORDERQ, + METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, + NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, + N, NCIT, NKEV, NCPS * LOGICAL VS, GSSCKM, ESCGSS, GEC, ULEC, RLXTOL, GEETE, + FUNCIT, NEWTON COMMON /COLCML/ VS, GSSCKM, ESCGSS, GEC, ULEC, RLXTOL, GEETE, + FUNCIT, NEWTON * REAL TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC COMMON /COLCMR/ TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC * SAVE /COLCMI/, /COLCML/, /COLCMR/ C I C PARAMETER SPECIFICATION: I C ----------------------- I INTEGER NEQN, IERROR LOGICAL LINEAR REAL T0, TN REAL CR(MR), WR(SR), LCR(MR,LR:MR,LR:SR), H(0:MAXNC), + UR(-NEQN:NEQN*(MR-LR+1)*MAXNC-1), WKAREA(*), URN(NEQN) EXTERNAL G, KC, DKCDY C H ENTRY: H(0:N-1) SUBINT. LENGTHS OF STEPS TAKEN IN THE I C COMPUTATION OF THE APPROX. SOL. BY "SOLVI2" I C UR ENTRY: UR(-NEQN:-1) SHOULD CONTAIN G(T0) I C EXIT: CONTAINS UR(TI+CRJ.HI) I=0:N-1, J=1:MR I C UR(TIJ) -> UR(NEQN*((MR-LR+1)*I+(J-LR))+(0:NEQN-1)) I C WKAREA WORKING STORAGE FOR SOLVING THE COLLOC. SYSTEM (CF.SOLSYS) I C AND FOR STORAGE OF INTERMEDIATE VECTOR RESULTS IN "SLQCE2" I C TN ENTRY: RIGHT ENDPOINT OF LAST SUBINTERVAL I C EXIT: AFTER NORMAL EXIT UNCHANGED, OTHERWISE LEFT ENDPOINT I C OF SUBINTERVAL AT WHICH ERROR OCCURRED I C URN EXIT: REFERENCE SOLUTION IN TN I C IERROR ENTRY: 0 I C EXIT: 0: OK I C 113: # KERNEL EVAL. > MAX. # KERN. EV. ALLOWED I C 114: # CPU-SECONDS > MAX. # CPU-SEC. ALLOWED I C OTHER: ERRORS FROM "SOLSYS" I C I C INVOKED BY: COLVI2 I C ---------- I C I C CHANGES IN COMMON VARIABLES: I C --------------------------- I C N USED TO KEEP TRACK OF THE # SUBINTERVALS ON WHICH THE I C INTEGRATION ALREADY HAS BEEN PERFORMED. ON EXIT "N" HAS THE I C SAME VALUE AS ON ENTRY. I C I C LOCAL VARIABLES: I C --------------- I INTEGER INDEXN, IW1,IW2,IW3, MW, NC LOGICAL LOBATR REAL HN C INDEXN POINTER TO REF.SOL. IN 1-ST (IF LOBATTO: 2-ND) COLLOC. POINT I C IN N-TH INTERVAL I C IW1,...,IW3 POINTERS TO WKAREA LOCATIONS (CF. "COLDOC" SUB I C "DISTRIBUTION WKAREA") I C MW DIMENSION OF SYSTEM OF COLLOCATION EQUATIONS I C NC # INTERVALS IN WHICH INTEGRATION INTERVAL HAS BEEN DIVIDED I C I C ---------------------------------------------------------------------I C INTEGER J, MRL, NCPJOB LOGICAL RSFAIL EXTERNAL NCPJOB * MRL = MR-LR+1 C C C DISTRIBUTE WORKING STORAGE "WKAREA" FOR LAGTERMS AND LIN.SYS.SOLVER MW = MRL*NEQN IW1 = 1 + MW IW2 = IW1 + MW IF (FUNCIT) THEN IW3 = IW2 ELSE IW3 = IW2 + MW*MW ENDIF C C C INITIALIZE LOOP CONSTANTS AND VARIABLES LOBATR = LR .EQ. 2 * NC = N N = 0 TN = T0 C C C LOOP UNTIL "TE" IS REACHED 100 CONTINUE RSFAIL = .FALSE. HN = H(N) C C CHECK # KEV IF (NKEV .GT. MAXKEV) GOTO 900 C C CHECK CPU-TIME IF (NCPJOB()-NCPS .GT. MAXCPS) GOTO 910 C INDEXN = NEQN*MRL*N C C C COMPUTE REF. SOL. WITH HIGHER ORDER COLLOCATION METHD DO 200 J = LR, MR CALL COPYV (UR(INDEXN-NEQN), NEQN, UR(INDEXN+NEQN*(J-LR))) 200 CONTINUE * 210 CALL SLQCE2 (TN,H, NEQN, G,KC,DKCDY,LINEAR, T0, + CR,WR,LCR, MR,SR,LR, LOBATR, TOLCIR, .FALSE., + WKAREA(1),WKAREA(1),WKAREA(IW1),WKAREA(IW2), MW, + WKAREA(IW3), UR, WKAREA(IW3), IERROR) IF (IERROR .EQ. 21 .AND. .NOT. RSFAIL) THEN C CONVERGENCE PROBLEMS, TRY AGAIN WITH BETTER ESTIMATE TO START WITH RSFAIL = .TRUE. IERROR = 0 GOTO 210 ELSE IF (IERROR .NE. 0) THEN GOTO 920 ENDIF C C LOOP IF "TE" NOT REACHED TN = TN+HN N = N+1 IF (N .LT. NC) GOTO 100 C C FINISHED CALL COPYV (UR(NEQN*MRL*N-NEQN), NEQN, URN) RETURN C C C ERROR RETURNS C C NUMBER OF KERNEL EVALUATIONS TOO LARGE 900 CALL ERRMSG ('NUMBER OF KERNEL EVALUATIONS EXCEEDS IOPT(6)') IERROR = 113 GOTO 920 C C TOO MUCH CPU-TIME USED 910 CALL ERRMSG ('CPU-TIME USED EXCEEDS IOPT(7)') IERROR = 114 C C ERROR WHILE SOLVING COLLOCATION EQUATION 920 CALL ERRMSG ('ERROR WHILE COMPUTING GLOBAL ERROR IN "TE"') N = NC * RETURN * END SUBROUTINE SLQCE2 + (TN, H, NEQN,G,KC,DKCDY,LINEAR, T0, C,W,LC, MM,SS,LL, LOBAT, + TOLCIT, LNP1FL, LAG, GLAG,CORR,DSYS,MW, WKAREA, U, LAGNP1, + IERROR) C C ---------------------------------------------------------------------I C PURPOSE: SOLVE SYSTEM OF COLLOC. EQ. FOR (REF) SOL. IN SUBINTERVAL I C ------- [TN,TN+HN]. APPROXIMATE LAG TERM WITH QUADRATURE. I C NOTE: IT IS POSSIBLE THAT "LAG" AND "GLAG" , AS WELL AS "WKAREA" AND I C ---- "LAGNP1" SHARE THE SAME MEMORY LOCATIONS. I C I C COMMON VARIABLES: I C ---------------- I INTEGER METH, M, S, L, ORDER, ORDERQ, METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, N, NCIT, NKEV, NCPS COMMON /COLCMI/ METH, M, S, L, ORDER, ORDERQ, + METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, + NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, + N, NCIT, NKEV, NCPS * SAVE /COLCMI/ C I C PARAMETER SPECIFICATION: I C ----------------------- I INTEGER NEQN, MM, SS, LL, MW, IERROR LOGICAL LINEAR, LOBAT, LNP1FL REAL TN, T0, TOLCIT REAL H(0:MAXNC), C(MM), W(SS), LC(MM,LL:MM,LL:SS), + LAG(0:NEQN*(MM-LL+1)-1), GLAG(MW), CORR(MW), DSYS(*), + WKAREA(*), U(-NEQN:NEQN*(MM-LL+1)*MAXNC-1), LAGNP1(*) EXTERNAL G, KC, DKCDY C I C TN LEFT ENDPOINT OF CURRENT SUBINTERVAL I C H LENGTH OF SUBINTERVALS [TI,TI+HI] I=0,...,N I C C -I COLLOC. VARS; I C W I EITHER C,W,LC OR CR,WR,LCR I C LC -I I C MM -I DIMENSION PARAMETERS FOR COLL. VARS; I C SS I EITHER M,S,L OR MR,SR,LR I C LL -I I C LOBAT TRUE IF COLL. PARAMETERS ARE LOBATTO POINTS I C TOLCIT TOLERANCE FOR CORRECTOR ITERATION PROCESS TO SOLVE COLL. EQ. I C EITHER TOLCIA OR TOLCIR. I C LNP1FL TRUE IF FCN(TN+HN) HAS TO BE STORED IN LAGNP1 I C TO BE USED FOR CHECK ON POL.SOL. IF APPROX.METH. IS GAUSS I C LAG WORKING STORAGE FOR THE APPROX. OF THE LAG TERM IN TN+CJ.HN, I C FOR J=LL,...,MM I C FCN(TNJ) -> LAG(NEQN*(J-LL)+(0:NEQN-1)) I C GLAG WORKING STORAGE FOR G(TNJ)+FCN(TNJ), J=LL,...,MM I C G(TNJ)+FCN(TNJ) -> GLAG(NEQN*(J-LL)+(1:NEQN)) I C CORR WORKING STORAGE FOR THE RIGHT HAND SIDE, OR THE SOLUTION, OF I C THE LINEAR SYSTEM IN THE CORRECTOR ITERATION PROCESS. I C DSYS WORKING STORAGE FOR THE DERIVATIVE OF THE SYSTEM OF LIN. EQ. I C IN THE NEWTON PROCESS. I C MW DIMENSION OF LIN. SYSTEM (NEQN.(MM-LL+1)) I C WKAREA WORKING STORAGE FOR TEMPORARY VECTORS (2*NEQN) I C ALSO USED AS WORKING STORAGE FOR SOLSYS. I C U ENTRY: CONTAINS THE APPROXIMATED SOLUTION IN TI+CJ.HI I C FOR I = 0,...,N-1; J=1,...,MM I C SHOULD CONTAIN AN INITIAL APPROXIMATION OF U(TN+CJ.HN) I C EXIT: APPROX. SOL. OF IN TI+CJ.HI, I=0,...,N, J=1,...,MM I C U(TI+CJ.HI) -> U(NEQN*((MM-LL+1)*I+(J-LL))+(0:NEQN-1)) I C LAGNP1 EXIT: CONTAINS APPROX. OF FCN(TN+HN), IF LNP1FL=.TRUE. I C OTHERWISE NOT USED I C IERROR ERROR COMPLETION CODE I C ENTRY: SHOULD CONTAIN 0 I C EXIT: 0: NO ERRORS I C OTHER: ERROR COMPLETION CODE OF SOLSYS I C I C INVOKED BY: SOLVI2, SGEVI2 I C ---------- I C I C CHANGES IN COMMON VARIABLES: I C --------------------------- I C NKEV ADDED: # KERNEL EVALUATIONS NEEDED TO COMPUTE LAG TERMS IN I C CURRENT SUBINTERVAL I C I C LOCAL VARIABLES: I C --------------- I INTEGER INDEXI REAL HI, HN, TI, TNJ C INDEXI POINTER TO APPROX. IN 1-ST (IF LOBATTO: 2-ND) COLLOC. POINT I C OF I-TH SUBINTERVAL I C TNJ TN + C(J).HN I C I C ---------------------------------------------------------------------I C INTEGER I, INDEXJ, I1, J, K, MMLL REAL CJHN * I1 = NEQN+1 MMLL = MM-LL+1 C C TN C APPROXIMATE LAG TERM FCN(TNJ) = INT KC(TN+CJ.HN),S,Y(S)) DS BY C T0 C N-1 S C LAG_J = SUM HI SUM WK.KC(TN+CJ.HN,TI+CK.HI,U_IK) C I=0 K=1 C COMPUTE GLAG_J = G(TN+CJ.HN) + LAG_J J=LL,...,MM C STORE KERNEL VECTORS IN WKAREA C NOTE: IT IS POSSIBLE THAT SOME G-FUNCTIONS (E.G. G(TN+HN)) ALREADY C ---- HAVE BEEN EVALUATED; PITY. HN = H(N) IF (.NOT. LOBAT) THEN DO 10 J = 1, MM INDEXJ = NEQN*(J-1) CJHN = C(J)*HN TNJ = TN + CJHN * CALL ZEROV (LAG(INDEXJ), NEQN) TI = T0 DO 20 I = 0, N-1 INDEXI = NEQN*MMLL*I HI = H(I) DO 30 K = 1, SS CALL KC(TNJ,TI+C(K)*HI,U(INDEXI+NEQN*(K-1)),WKAREA) CALL ADDABV (LAG(INDEXJ), NEQN, HI*W(K), WKAREA) 30 CONTINUE TI = TI + HI 20 CONTINUE IF (LNP1FL .AND. J.EQ.MM) + CALL COPYV (LAG(INDEXJ), NEQN, LAGNP1) CALL G(TNJ,WKAREA) CALL ADDV (GLAG(INDEXJ+1), NEQN, 1.0,WKAREA,1.0,LAG(INDEXJ)) 10 CONTINUE NKEV = NKEV + MM*N*SS * ELSE IF (N .GT. 0) THEN C LOBATTO, N > 0 DO 50 J = 2, MM INDEXJ = NEQN*(J-2) CJHN = C(J)*HN TNJ = TN + CJHN * CALL KC(TNJ,T0,U(-NEQN),WKAREA) CALL ADDV (LAG(INDEXJ), NEQN, H(0)*W(1),WKAREA, 0.0,WKAREA) TI = T0 DO 60 I = 0, N-2 INDEXI = NEQN*MMLL*I HI = H(I) DO 70 K = 2, MM-1 CALL KC(TNJ,TI+C(K)*HI,U(INDEXI+NEQN*(K-2)),WKAREA) CALL ADDABV (LAG(INDEXJ), NEQN, HI*W(K), WKAREA) 70 CONTINUE CALL KC(TNJ,TI+HI,U(INDEXI+NEQN*(MM-2)),WKAREA) CALL ADDABV (LAG(INDEXJ), NEQN, (HI+H(I+1))*W(MM),WKAREA) TI = TI + HI 60 CONTINUE INDEXI = NEQN*MMLL*(N-1) HI = H(N-1) DO 80 K = 2, MM CALL KC(TNJ,TI+C(K)*HI,U(INDEXI+NEQN*(K-2)),WKAREA) CALL ADDABV (LAG(INDEXJ), NEQN, HI*W(K), WKAREA) 80 CONTINUE C C WKAREA(1:NEQN) CONTAINS K(TNJ,TN,U_N1) CALL G(TNJ,WKAREA(I1)) CALL ADDV (GLAG(INDEXJ+1),NEQN, 1.0,WKAREA(I1), + 1.0,LAG(INDEXJ)) C C ADD FIRST TERM OF SUM THAT APPROXIMATES (TN,TN+HN) INT K(TNJ,..) C I.E. CJ.HN.W1.K(TNJ,TN,U(TN)) CALL ADDABV (GLAG(INDEXJ+1), NEQN, CJHN*W(1), WKAREA) 50 CONTINUE NKEV = NKEV + (MM-1)*(1+N*(MM-1)) ELSE C LOBATTO, FIRST STEP DO 90 J = 2, MM INDEXJ = NEQN*(J-2) CJHN = C(J)*HN TNJ = TN + CJHN * CALL ZEROV (LAG(INDEXJ), NEQN) CALL KC(TNJ,TN,U(-NEQN),WKAREA) CALL G(TNJ,GLAG(INDEXJ+1)) C C ADD FIRST TERM OF SUM THAT APPROXIMATES (TN,TN+HN) INT K(TNJ,..) C I.E. CJ.HN.W1.K(TNJ,TN,U(TN)) CALL ADDABV (GLAG(INDEXJ+1), NEQN, CJHN*W(1), WKAREA) 90 CONTINUE NKEV = NKEV + MM-1 ENDIF C C C SOLVE, BY FUNCTIONAL OR NEWTON ITER., THE SYSTEM OF COLLOC. EQUATIONS C U_NJ - GLAG_J - C S C - HN.SUM CJ.WK.KC(TN+CJ.HN,TN+CJ.CK.HN,U(TN+CJ.CK.HN)) = 0 C K=1 J=LL,...,MM C CALL SOLSYS (TN,HN, NEQN,KC,DKCDY,LINEAR, C,W,LC, MM,SS,LL, + TOLCIT, GLAG,CORR,DSYS,MW, WKAREA, + U(NEQN*(MMLL*N+(1-LL))), IERROR) * RETURN END SUBROUTINE SLICE2 + (TN,HNM1,HN, NEQN,G,KC,DKCDY,LINEAR, C, CB, CR,WR,LCR,LOBATR, + LAGSAV, GLAGR,CORR,DSYS,MW, WKAREA, URN, IERROR) C C ---------------------------------------------------------------------I C PURPOSE: SOLVE SYSTEM OF COLLOC. EQ. FOR REF.SOL. IN [TN,TN+HN]. I C ------- COMPUTE LAG TERM IN (TN+CRJ.HN) WITH LAGRANGE INTERPOLATION I C OVER TWO INTERVALS. I C NOTE: CB AND CORR (AND POSS. DSYS AND WKAREA) SHARE MEMORY LOCATIONS I C ---- I C I C COMMON VARIABLES: I C ---------------- I INTEGER METH, M, S, L, ORDER, ORDERQ, METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, N, NCIT, NKEV, NCPS COMMON /COLCMI/ METH, M, S, L, ORDER, ORDERQ, + METHR, MR, SR, LR, ORDERR, + ERRWGT, NHFAIL, + NERR, NWIR, NSAV, + MAXNC, MAXKEV, MAXCPS, + N, NCIT, NKEV, NCPS * REAL TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC COMMON /COLCMR/ TOLLE, TOLCIA, TOLCIR, HMIN, HMAX, HC * SAVE /COLCMI/, /COLCMR/ C I C PARAMETER SPECIFICATION: I C ----------------------- I INTEGER NEQN, MW, IERROR LOGICAL LINEAR, LOBATR REAL TN, HNM1, HN REAL C(M), CB(2*M-L+1), CR(MR), WR(SR), LCR(MR,LR:MR,LR:SR), + LAGSAV(0:NEQN*(2*M-L+1)-1), GLAGR(MW), CORR(MW), + DSYS(*), WKAREA(*), URN(0:NEQN*MR-1) EXTERNAL G, KC, DKCDY C I C TN LEFT ENDPOINT OF CURRENT SUBINTERVAL I C HNM1 LENGTH OF SUBINTERVAL [T(N-1),TN] I C HN LENGTH OF CURRENT SUBINTERVAL I C CB WORKING STORAGE FOR COLL.POINTS NEEDED FOR LAG TERM INTERPOL. I C LAGSAV LAG TERM APPROX. IN THE PREVIOUS AND CURRENT INTERVAL. I C ENTRY: I C FCN(TNM1+CJ.HNM1) -> LAGSAV(NEQN*(J-1):NEQN*J-1) J=1..M I C FCN(TN+CJ.HN) -> LAGSAV(NEQN*(M+J-L):NEQN*(M+J-L+1)-1) J=L..M I C GLAGR WORKING STORAGE FOR G(TN+CRJ.HN) + FCN(TN+CRJ.HN) J=LR..MR I C CORR WORKING STORAGE FOR THE RIGHT HAND SIDE, OR THE SOLUTION, OF I C THE LINEAR SYSTEM IN THE CORRECTOR ITERATION PROCESS. I C DSYS WORKING STORAGE FOR THE DERIVATIVE OF THE SYSTEM OF LIN. EQ. I C I