C ALGORITHM 688, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 17, NO. 2, PP. 153-166. JUNE, 1991. C The programs on the tape are as follows : C C 1. README This file. C C 2. ex1drvsp.f A single precision driving routine for example 1 C of Sincovec and Madsen. Includes the subroutines : C C MAIN C SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) C SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) C SUBROUTINE UINIT(X,U,NPDE) C SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) C C 3. ex2drvsp.f A single precision driving routine for example 2 C of Sincovec and Madsen. Includes the subroutines. C C MAIN C SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) C SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) C SUBROUTINE UINIT(X,U,NPDE) C SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) C REAL FUNCTION EXACTV(X,T) C C 4. pdecolsp.f The single precision version of the modified pdecol. C C The subroutines are in the following order : (only the eight routines C above the line of asterisks differ from the original code, with C CRDCMP and CRSLVE replacing DECB and SOLB. C C SUBROUTINE PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF, C SUBROUTINE ADDA(PW,A,BC,NPDE) C SUBROUTINE DIFFUN (N, T, Y, YDOT, IFLAG, PW, IPIV, WORK, IWORK) C SUBROUTINE INITAL(K,A,RHS,X,XT,XC,PW,IPIV,ILEFT) C SUBROUTINE PSETIB (C, PW, N0, CON, MITER, IFLAG, A, ILEFT, XC, C SUBROUTINE STIFIB (N0,Y,YMAX,ERROR,SAVE1,SAVE2,SAVE3, C SUBROUTINE CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK, C SUBROUTINE CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK, C C ************************************************************** C C SUBROUTINE BSPLVD ( XT, K, X, ILEFT, VNIKX, NDERIV ) C SUBROUTINE BSPLVN ( XT, JHIGH, INDEX, X, ILEFT, VNIKX ) C SUBROUTINE INTERV ( XT, LXT, X, ILEFT, MFLAG ) C SUBROUTINE COLPNT(X, XC, XT) C SUBROUTINE COSET (METH, NQ, EL, TQ) C SUBROUTINE DIFFF(T,X,IPT,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE,CMAX, C SUBROUTINE EVAL(ICPT,NPDE,C,UVAL,A,ILEFT) C SUBROUTINE GFUN (T,C,UDOT,NPDE,NCPTS,A,BC,DBDU,DBDUX,DZDT, C SUBROUTINE INTERP (TOUT, Y, N0, Y0) C SUBROUTINE RES(T,H,C,V,R,NPDE,NCPTS,A,ILEFT,BC,DBDU,DBDUX,DZDT, C SUBROUTINE VALUES(X,USOL,SCTCH,NDIM1,NDIM2,NPTS,NDERV,WORK) C C 5. ex1drvdp.f A double precision driving routine for example 1 C of Sincovec and Madsen. Includes the subroutines. C C MAIN C SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) C SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) C SUBROUTINE UINIT(X,U,NPDE) C SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) C C 6. ex2drvdp.f A double precision driving routine for example 2 C of Sincovec and Madsen. Includes the subroutines. C C MAIN C SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) C SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) C SUBROUTINE UINIT(X,U,NPDE) C SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) C DOUBLE PRECISION FUNCTION EXACTV(X,T) C C C 7. pdecoldp.f The double precision version of the modified pdecol. C C The subroutines are in the following order : (only the eight routines C above the line of asterisks differ from the original code, with C CRDCMP and CRSLVE replacing DECB and SOLB. C C SUBROUTINE PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF, C SUBROUTINE ADDA(PW,A,BC,NPDE) C SUBROUTINE DIFFUN (N, T, Y, YDOT, IFLAG, PW, IPIV, WORK, IWORK) C SUBROUTINE INITAL(K,A,RHS,X,XT,XC,PW,IPIV,ILEFT) C SUBROUTINE PSETIB (C, PW, N0, CON, MITER, IFLAG, A, ILEFT, XC, C SUBROUTINE STIFIB (N0,Y,YMAX,ERROR,SAVE1,SAVE2,SAVE3, C SUBROUTINE CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK, C SUBROUTINE CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK, C C ************************************************************** C C SUBROUTINE BSPLVD ( XT, K, X, ILEFT, VNIKX, NDERIV ) C SUBROUTINE BSPLVN ( XT, JHIGH, INDEX, X, ILEFT, VNIKX ) C SUBROUTINE INTERV ( XT, LXT, X, ILEFT, MFLAG ) C SUBROUTINE COLPNT(X, XC, XT) C SUBROUTINE COSET (METH, NQ, EL, TQ) C SUBROUTINE DIFFF(T,X,IPT,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE,CMAX, C SUBROUTINE EVAL(ICPT,NPDE,C,UVAL,A,ILEFT) C SUBROUTINE GFUN (T,C,UDOT,NPDE,NCPTS,A,BC,DBDU,DBDUX,DZDT, C SUBROUTINE INTERP (TOUT, Y, N0, Y0) C SUBROUTINE RES(T,H,C,V,R,NPDE,NCPTS,A,ILEFT,BC,DBDU,DBDUX,DZDT, C SUBROUTINE VALUES(X,USOL,SCTCH,NDIM1,NDIM2,NPTS,NDERV,WORK) C C THIS SINGLE PRECISION DRIVER MAY BE USED TO RUN THE FIRST EXAMPLE C IN SECTION 7 OF SINCOVEC AND MADSEN, TOMS, 5, 1979. C THE INPUT REQUIRED CONSISTS OF : C NINT : THE NUMBER OF SUBINTERVALS TO BE USED IN THE SPACE C VARIABLE. C KORD : THE ORDER OF THE B-SPLINES TO BE USED ( = DEGREE-1). C SEVERAL SETS OF INPUT VALUES MAY BE GIVEN IN FORMAT 2I4. C REAL XLEFT,DTUSED,XBKPT(121),U(2,121),SCTCH(121), * WORK(8000),T0,TOUT,DT,EPS,DX INTEGER NQ,NSTEPS,NF,NJ,IWORK(1000),NPDE,NINT,NPTS,KORD,NCC, * MF,I,K,INDEX COMMON /ENDPT/ XLEFT COMMON /GEAR0/ DTUSED,NQ,NSTEPS,NF,NJ NPDE = 2 1 CONTINUE READ(5,9100,END=9000) NINT, KORD NPTS = NINT + 1 NCC = 2 T0 = 0.0 TOUT = 1.E-3 DT = 1.E-7 EPS = 1.0E-4 MF = 21 WRITE(6,9200)NINT,KORD,EPS INDEX = 1 IWORK(1) = 8000 IWORK(2) = 1000 DX = 1.0/DBLE(NPTS-1) DO 10 I = 1,NPTS XBKPT(I) = DBLE(I-1) * DX 10 CONTINUE XLEFT = XBKPT(1) 20 CONTINUE CALL PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,INDEX, * WORK,IWORK) IF ( INDEX .NE. 0 ) GO TO 70 WRITE(6,9150) TOUT,DTUSED,NSTEPS 9150 FORMAT(' T = ',E10.3,' DT = ',E10.3,' TOTAL STEPS = ',I6) CALL VALUES(XBKPT,U,SCTCH,NPDE,NPTS,NPTS,1,WORK) DO 60 K = 1,NPDE WRITE(6,9300)K WRITE(6,9400)(U(K,I),I=1,NPTS) 60 CONTINUE TOUT = TOUT * 10.0 IF ( TOUT .LT. 11.0 ) GO TO 20 70 WRITE(6,9500) INDEX GO TO 1 9000 STOP 9100 FORMAT(2I4) 9200 FORMAT(' NO. OF SUBINTERVALS = ',I3,' KORD = ',I2, *' EPS = ',D10.2) 9300 FORMAT(10X,'PDE COMPONENT = ',I3) 9400 FORMAT(10X,5E12.4) 9500 FORMAT(' INDEX = ',I3) END SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) C C THIS IS THE USER SUPPLIED SUBROUTINE TO SPECIFY THE DIFFERENTIAL C EQUATIONS. C REAL U(NPDE),UX(NPDE),UXX(NPDE),FVAL(NPDE),T,X INTEGER NPDE FVAL(1) = U(2)*U(2)*UXX(1) - U(1)*U(2) - U(1)*U(1) + 10.0 * + 2.0*U(2)*UX(2)*UX(1) FVAL(2) = U(1)*U(1)*UXX(2) + U(1)*U(2) - U(2)*U(2) * + UXX(1) + 2.0*U(1)*UX(1)*UX(2) RETURN END SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) C C THIS ROUTINE SPECIFIES THE BOUNDARY CONDITION EQUATIONS. C INTEGER NPDE REAL T,X,U(NPDE),UX(NPDE),DZDT(NPDE), * DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE), * XLEFT COMMON /ENDPT/ XLEFT IF ( X .NE. XLEFT ) GO TO 10 DBDU(1,1) = 1.0 DBDU(1,2) = 0.0 DBDU(2,1) = 0.0 DBDU(2,2) = 1.0 DBDUX(1,1) = 0.0 DBDUX(1,2) = 0.0 DBDUX(2,1) = 0.0 DBDUX(2,2) = 0.0 DZDT(1) = 0.0 DZDT(2) = 0.0 RETURN 10 CONTINUE DBDU(1,1) = U(2)*COS(U(1)*U(2)) DBDU(1,2) = U(1)*COS(U(1)*U(2)) DBDU(2,1) = U(2)*SIN(U(1)*U(2)) DBDU(2,2) = U(1)*SIN(U(1)*U(2)) DBDUX(1,1) = 1.0 DBDUX(1,2) = 0.0 DBDUX(2,1) = 0.0 DBDUX(2,2) = 1.0 DZDT(1) = 0.0 DZDT(2) = 0.0 RETURN END SUBROUTINE UINIT(X,U,NPDE) C C UINIT GIVES THE INITIAL CONDITIONS AT T=T0. C INTEGER NPDE REAL X,U(NPDE) U(1) = 0.5*(X + 1.0) U(2) = 4.0*ATAN(1.0) RETURN END SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) C C THIS IS THE OPTIONAL ROUTINE PROVIDED IF THE USER WISHES TO C SUPPLY AN ANALYTIC JACOBIAN. C INTEGER NPDE REAL T,X,U(NPDE),UX(NPDE),UXX(NPDE), * DFDU(NPDE,NPDE),DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE) DFDU(1,1) = -U(2) -2.0*U(1) DFDU(1,2) = 2.0*U(2)*UXX(1) - U(1) + 2.0*UX(2)*UX(1) DFDU(2,1) = 2.0*U(1)*UXX(2) + U(2) + 2.0*UX(1)*UX(2) DFDU(2,2) = U(1) - 2.0*U(2) DFDUX(1,1) = 2.0*U(2)*UX(2) DFDUX(1,2) = 2.0*U(2)*UX(1) DFDUX(2,1) = 2.0*U(1)*UX(2) DFDUX(2,2) = 2.0*U(1)*UX(1) DFDUXX(1,1) = U(2)*U(2) DFDUXX(1,2) = 0.0 DFDUXX(2,1) = 1.0 DFDUXX(2,2) = U(1)*U(1) RETURN END C C THIS SINGLE PRECISION DRIVER MAY BE USED TO RUN THE SECOND EXAMPLE C IN SECTION 7 OF SINCOVEC AND MADSEN, TOMS, 5, 1979. C THE INPUT REQUIRED CONSISTS OF : C NINT : THE NUMBER OF SUBINTERVALS TO BE USED IN THE C VARIABLE. C KORD : THE ORDER OF THE B-SPLINES TO BE USED ( = DEG C SEVERAL SETS OF INPUT VALUES MAY BE GIVEN IN FORMAT 2I4. C REAL XLEFT,DTUSED,XBKPT(121),U(1,121),SCTCH(121), * WORK(8000),T0,TOUT,DT,EPS,DX,PI,EXACTV,ERROR INTEGER NQ,NSTEPS,NF,NJ,IWORK(1000),NPDE,NINT,NPTS,KORD,NCC, * MF,I,K,INDEX COMMON /ENDPT/ XLEFT COMMON /GEAR0/ DTUSED,NQ,NSTEPS,NF,NJ COMMON /PIVAL/ PI PI = ACOS(-1.0E+0) NPDE = 1 NCC = 2 1 CONTINUE READ(5,9100,END=9000)NINT,KORD NPTS = NINT + 1 T0 = 0.0 TOUT = 1.0E+0 DT = 0.1E-8 EPS = 1.E-4 MF = 21 INDEX = 1 IWORK(1) = 8000 IWORK(2) = 1000 DX = 1.0/REAL(NPTS-1) DO 10 I = 1,NPTS XBKPT(I) = REAL(I-1) * DX 10 CONTINUE XLEFT = XBKPT(1) WRITE(6,9200) NINT, KORD, EPS 20 CONTINUE CALL PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,INDEX, * WORK,IWORK) IF ( INDEX .NE. 0 ) GO TO 70 WRITE(6,9150) TOUT,DTUSED,NSTEPS CALL VALUES(XBKPT,U,SCTCH,NPDE,NPTS,NPTS,1,WORK) ERROR = 0.E0 DO 60 K = 1,NPDE DO 41 I = 1,NPTS ERROR = ERROR + (U(K,I)-EXACTV(XBKPT(I),TOUT))**2 41 CONTINUE 60 CONTINUE WRITE(6,9300)SQRT(ERROR) TOUT = TOUT * 10.0 IF ( TOUT .LT. 1.0) GO TO 20 70 CONTINUE WRITE(6,9400) INDEX GO TO 1 9000 STOP 9100 FORMAT(2I4) 9150 FORMAT(' T = ',E10.3,' DT = ',E10.3,' TOTAL STEPS = ',I6) 9200 FORMAT(' NINT = ',I3, ' KORD = ', I3, ' EPS = ', E10.3) 9300 FORMAT(' L2 ERROR = ', D24.15) 9400 FORMAT(' INDEX = ',I3) END SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) C C THIS IS THE USER SUPPLIED SUBROUTINE TO SPECIFY THE DIFFERENTIAL C EQUATIONS. C REAL U(NPDE),UX(NPDE),UXX(NPDE),FVAL(NPDE),T,X INTEGER NPDE REAL PI COMMON/PIVAL/PI FVAL(1) = UXX(1) + PI * PI * SIN( PI * X ) RETURN END SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) C C THIS ROUTINE SPECIFIES THE BOUNDARY CONDITION EQUATIONS. C INTEGER NPDE REAL T,X,U(NPDE),UX(NPDE),DZDT(NPDE), * DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE), * XLEFT COMMON /ENDPT/ XLEFT IF ( X .NE. XLEFT ) GO TO 10 DBDU(1,1) = 1.0 DBDUX(1,1) = 0.0 DZDT(1) = 0.0 RETURN 10 CONTINUE DBDU(1,1) = 1.0 DBDUX(1,1) = 0.0 DZDT(1) = 0.0 RETURN END SUBROUTINE UINIT(X,U,NPDE) C C UINIT GIVES THE INITIAL CONDITIONS AT T=T0. C INTEGER NPDE REAL X,U(NPDE) U(1) = 1.0 RETURN END SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) C C THIS IS THE OPTIONAL ROUTINE PROVIDED IF THE USER WISHES TO C SUPPLY AN ANALYTIC JACOBIAN. C INTEGER NPDE REAL T,X,U(NPDE),UX(NPDE),UXX(NPDE), * DFDU(NPDE,NPDE),DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE) DFDU(1,1) = 0.0 DFDUX(1,1) = 0.0 DFDUXX(1,1) = 1.0 RETURN END REAL FUNCTION EXACTV(X,T) C C THIS FUNCTION PROVIDES THE EXACT ANSWER C TO THE DIFFERENTIAL EQUATION. C REAL X,T REAL PI COMMON/PIVAL/PI EXACTV = 1.0 + SIN(PI*X)*(1.0 - EXP(-PI*PI*T)) RETURN END SUBROUTINE PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF, * INDEX,WORK,IWORK) C C----------------------------------------------------------------------- C----------------------------------------------------------------------- C C THIS IS THE AUG 1987 VERSION OF PDECOL, WHICH MAKES USE OF THE ALMOST C BLOCK-DIAGONAL LINEAR EQUATION SOLVER COLROW (COMPRISING CRDCMP, FOR C THE FACTORIZATION, AND CRSLVE FOR THE SOLUTION PHASE). C C THIS PACKAGE WAS CONSTRUCTED SO AS TO CONFORM TO AS MANY ANSI-FORTRAN C RULES AS WAS CONVIENTLY POSSIBLE. THE FORTRAN USED VIOLATES ANSI C STANDARDS IN THE TWO WAYS LISTED BELOW.... C C 1. SUBSCRIPTS OF THE GENERAL FORM C*V1 + V2 + V3 ARE USED C (POSSIBLY IN A PERMUTED ORDER), WHERE C IS AN INTEGER CONSTANT C AND V1, V2, AND V3 ARE INTEGER VARIABLES. C C 2. ARRAY NAMES APPEAR SINGLY IN DATA STATEMENTS IN THE ROUTINES C BSPLVN AND COSET. C C----------------------------------------------------------------------- C----------------------------------------------------------------------- C C----------------------------------------------------------------------- C PDECOL IS THE DRIVER ROUTINE FOR A SOPHISTICATED PACKAGE OF C SUBROUTINES WHICH IS DESIGNED TO SOLVE THE GENERAL SYSTEM OF C NPDE NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS OF AT MOST SECOND C ORDER ON THE INTERVAL (XLEFT,XRIGHT) FOR T .GT. T0 WHICH IS OF THE C FORM.... C C DU/DT = F( T, X, U, UX, UXX ) C C WHERE C C U = ( U(1), U(2), ... , U(NPDE) ) C UX = ( UX(1), UX(2), ... , UX(NPDE) ) C UXX = (UXX(1),UXX(2), ... ,UXX(NPDE) ) . C C EACH U(K) IS A FUNCTION OF THE SCALAR QUANTITIES T AND X. C UX(K) REPRESENTS THE FIRST PARTIAL DERIVATIVE OF U(K) WITH RESPECT C TO THE VARIABLE X, UXX(K) REPRESENTS THE SECOND PARTIAL DERIVATIVE C OF U(K) WITH RESPECT TO THE VARIABLE X, AND DU/DT IS THE VECTOR OF C PARTIAL DERIVATIVES OF U WITH RESPECT TO THE TIME VARIABLE T. C F REPRESENTS AN ARBITRARY VECTOR VALUED FUNCTION WHOSE NPDE C COMPONENTS DEFINE THE RESPECTIVE PARTIAL DIFFERENTIAL EQUATIONS OF C THE PDE SYSTEM. SEE SUBROUTINE F DESCRIPTION BELOW. C C BOUNDARY CONDITIONS C C IN THE NEW VERSION OF PDECOL, 2 BOUNDARY CONDITIONS C ARE REQUIRED FOR EACH PDE IN THE SYSTEM. THESE ARE IMPOSED AT XLEFT C AND RIGHT AND EACH MUST BE OF THE FORM.... C C B(U,UX) = Z(T) C C WHERE B AND Z ARE ARBITRARY VECTOR VALUED FUNCTIONS WITH C NPDE COMPONENTS AND U, UX, AND T ARE AS ABOVE. THESE BOUNDARY C CONDITIONS MUST BE CONSISTENT WITH THE INITIAL CONDITIONS WHICH ARE C DESCRIBED NEXT. C C INITIAL CONDITIONS C C EACH SOLUTION COMPONENT U(K) IS ASSUMED TO BE A KNOWN (USER C PROVIDED) FUNCTION OF X AT THE INITIAL TIME T = T0. THE C INITIAL CONDITION FUNCTIONS MUST BE CONSISTENT WITH THE BOUNDARY C CONDITIONS ABOVE, I.E. THE INITIAL CONDITION FUNCTIONS MUST C SATISFY THE BOUNDARY CONDITIONS FOR T = T0. SEE SUBROUTINE UINIT C DESCRIPTION BELOW. C----------------------------------------------------------------------- C C REQUIRED USER SUPPLIED SUBROUTINES C C THE USER IS REQUIRED TO CONSTRUCT THREE SUBPROGRAMS AND A MAIN C PROGRAM WHICH DEFINE THE PDE PROBLEM WHOSE SOLUTION IS TO BE C ATTEMPTED. THE THREE SUBPROGRAMS ARE... C C 1) SUBROUTINE F( T, X, U, UX, UXX, FVAL, NPDE ) C DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE) C THIS ROUTINE DEFINES THE DESIRED PARTIAL DIFFERENTIAL C EQUATIONS TO BE SOLVED. THE PACKAGE PROVIDES VALUES OF THE C INPUT SCALARS T AND X AND INPUT ARRAYS (LENGTH NPDE) U, UX, C AND UXX, AND THE USER MUST CONSTRUCT THIS ROUTINE TO COMPUTE C THE OUTPUT ARRAY FVAL (LENGTH NPDE) WHICH CONTAINS THE C CORRESPONDING VALUES OF THE RIGHT HAND SIDES OF THE DESIRED C PARTIAL DIFFERENTIAL EQUATIONS, I.E. C C FVAL(K) = THE VALUE OF THE RIGHT HAND SIDE OF THE K-TH PDE IN C THE PDE SYSTEM ABOVE, FOR K = 1 TO NPDE. C C THE INCOMING VALUE OF THE SCALAR QUANTITY X WILL BE A C COLLOCATION POINT VALUE (SEE INITAL AND COLPNT) AND THE C INCOMING VALUES IN THE ARRAYS U, UX AND UXX CORRESPOND TO THIS C POINT X AND TIME T. C RETURN C END C C 2) SUBROUTINE BNDRY( T, X, U, UX, DBDU, DBDUX, DZDT, NPDE ) C DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE) C DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE) C THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH NEEDED C INFORMATION ABOUT THE BOUNDARY CONDITION FUNCTIONS B AND Z C ABOVE. THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES C T, X, U, AND UX, AND THE USER IS TO DEFINE THE CORRESPONDING C OUTPUT VALUES OF THE DERIVATIVES OF THE FUNCTIONS B AND Z C WHERE.... C DBDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO C THE J-TH VARIABLE U(J). C DBDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO C THE J-TH VARIABLE UX(J). C DZDT(K) = DERIVATIVE OF THE K-TH COMPONENT OF THE VECTOR C FUNCTION Z(T) ABOVE WITH RESPECT TO THE C VARIABLE T. C NOTE... THE INCOMING VALUE OF X WILL BE EITHER XLEFT OR XRIGHT. C WE REQUIRE A BOUNDARY CONDITION TO BE IMPOSED AT BOTH ENDPOINTS C OF EACH PDE IN THE SYSTEM. THUS, FOR SAY THE K-TH PDE, WE REQUI C AT BOTH OF THE ENDPOINTS, XLEFT AND XRIGHT, THAT ONE OF DBDU(K, C OR DBDUX(K,K) SHOULD BE NON-ZERO WHEN BNDRY IS CALLED. C THIS ROUTINE CAN BE STRUCTURED AS FOLLOWS... C THE COMMON BLOCK /ENDPT/ IS NOT A PART OF PDECOL AND C MUST BE SUPPLIED AND DEFINED BY THE USER. C COMMON /ENDPT/ XLEFT C IF( X .NE. XLEFT ) GO TO 10 C HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J), C AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE LEFT BOUNDARY POINT C X = XLEFT. C RETURN C 10 CONTINUE C HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J), C AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE RIGHT BOUNDARY POINT C X = XRIGHT. C RETURN C END C C 3) SUBROUTINE UINIT( X, U, NPDE ) C DIMENSION U(NPDE) C THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH THE C NEEDED INITIAL CONDITION FUNCTION VALUES. THE PACKAGE C PROVIDES A VALUE OF THE INPUT VARIABLE X, AND THE USER IS TO C DEFINE THE PROPER INITIAL VALUES (AT T = T0) FOR ALL OF THE C PDE COMPONENTS, I.E. C U(K) = DESIRED INITIAL VALUE OF PDE COMPONENT U(K) AT C X AND T = T0 FOR K = 1 TO NPDE. C NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT C VALUE. THE INITIAL CONDITIONS AND BOUNDARY CONDITIONS C MUST BE CONSISTENT (SEE ABOVE). C RETURN C END C----------------------------------------------------------------------- C C OPTIONAL USER SUPPLIED SUBROUTINE C C IF THE USER DESIRES TO USE THE MF = 11 OR 21 OPTION IN ORDER TO SAVE C ABOUT 10-20 PERCENT IN EXECUTION TIME (SEE BELOW), THEN THE USER MUST C PROVIDE THE FOLLOWING SUBROUTINE WHICH PROVIDES INFORMATION ABOUT THE C DERIVATIVES OF THE FUNCTION F ABOVE. THIS PROVIDES FOR MORE EFFICIENT C JACOBIAN MATRIX GENERATION. ON MOST COMPUTER SYSTEMS, THE USER WILL C BE REQUIRED TO SUPPLY THIS SUBROUTINE AS A DUMMY SUBROUTINE IF THE C OPTIONS MF = 12 OR 22 ARE USED (SEE BELOW). C C 1) SUBROUTINE DERIVF( T, X, U, UX, UXX, DFDU, DFDUX, DFDUXX, NPDE ) C DIMENSION U(NPDE), UX(NPDE), UXX(NPDE) C DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE) C THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES T, X, U, UX, C AND UXX, AND THE USER SHOULD CONSTRUCT THIS ROUTINE TO PROVIDE C THE FOLLOWING CORRESPONDING VALUES OF THE OUTPUT ARRAYS C DFDU, DFDUX, AND DFDUXX FOR K,J = 1 TO NPDE... C DFDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C PDE DEFINING FUNCTION F WITH RESPECT TO THE C VARIABLE U(J). C DFDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C PDE DEFINING FUNCTION F WITH RESPECT TO THE C VARIABLE UX(J). C DFDUXX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C PDE DEFINING FUNCTION F WITH RESPECT TO THE C VARIABLE UXX(J). C NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT C VALUE. C RETURN C END C----------------------------------------------------------------------- C C METHODS USED C C THE PACKAGE PDECOL IS BASED ON THE METHOD OF LINES AND USES A C FINITE ELEMENT COLLOCATION PROCEDURE (WITH PIECEWISE POLYNOMIALS C AS THE TRIAL SPACE) FOR THE DISCRETIZATION OF THE SPATIAL VARIABLE C X. THE COLLOCATION PROCEDURE REDUCES THE PDE SYSTEM TO A SEMI- C DISCRETE SYSTEM WHICH THEN DEPENDS ONLY ON THE TIME VARIABLE T. C THE TIME INTEGRATION IS THEN ACCOMPLISHED BY USE OF SLIGHTLY C MODIFIED STANDARD TECHNIQUES (SEE REFS. 1,2). C C PIECEWISE POLYNOMIALS C C THE USER IS REQUIRED TO SELECT THE PIECEWISE POLYNOMIAL SPACE C WHICH IS TO BE USED TO COMPUTE HIS APPROXIMATE SOLUTION. FIRST, THE C ORDER, KORD, OF THE POLYNOMIALS TO BE USED MUST BE SPECIFIED C (KORD = POLYNOMIAL DEGREE + 1). NEXT, THE NUMBER OF PIECES C (INTERVALS), NINT, INTO WHICH THE SPATIAL DOMAIN (XLEFT,XRIGHT) IS C TO BE DIVIDED, IS CHOSEN. THE NINT + 1 DISTINCT BREAKPOINTS OF C THE DOMAIN MUST BE DEFINED AND SET INTO THE ARRAY XBKPT IN C STRICTLY INCREASING ORDER, I.E. C XLEFT=XBKPT(1) .LT. XBKPT(2) .LT. ... .LT. XBKPT(NINT+1)=XRIGHT. C THE APPROXIMATE SOLUTION AT ANY TIME T WILL BE A POLYNOMIAL OF C ORDER KORD OVER EACH SUBINTERVAL (XBKPT(I),XBKPT(I+1)). THE C NUMBER OF CONTINUITY CONDITIONS, NCC, TO BE IMPOSED ACROSS ALL OF C THE BREAKPOINTS IS THE LAST PIECE OF USER SUPPLIED DATA WHICH IS C REQUIRED TO UNIQUELY DETERMINE THE DESIRED PIECEWISE POLYNOMIAL C SPACE. FOR EXAMPLE, NCC = 2 WOULD REQUIRE THAT THE APPROXIMATE C SOLUTION (MADE UP OF THE SEPARATE POLYNOMIAL PIECES) AND ITS FIRST C SPATIAL DERIVATIVE BE CONTINUOUS AT THE BREAKPOINTS AND HENCE ON C THE ENTIRE DOMAIN (XLEFT,XRIGHT). NCC = 3 WOULD REQUIRE THAT THE C APPROXIMATE SOLUTION AND ITS FIRST AND SECOND SPATIAL DERIVATIVES C BE CONTINUOUS AT THE BREAKPOINTS, ETC. THE DIMENSION OF THIS LINEAR C SPACE IS KNOWN AND FINITE AND IS NCPTS = KORD*NINT - NCC*(NINT-1). C THE WELL-KNOWN B-SPLINE BASIS (SEE REF. 3) FOR THIS SPACE IS USED C BY PDECOL AND IT CONSISTS OF NCPTS KNOWN PIECEWISE POLYNOMIAL C FUNCTIONS BF(I,X), FOR I=1 TO NCPTS, WHICH DO NOT DEPEND ON THE C TIME VARIABLE T. WE WISH TO EMPHASIZE THAT THE PIECEWISE POLYNOMIAL C SPACE USED IN PDECOL (WHICH IS SELECTED BY THE USER) WILL DETERMINE C THE MAGNITUDE OF THE SPATIAL DISCRETIZATION ERRORS IN THE COMPUTED C APPROXIMATE SOLUTION. THE PACKAGE HAS NO CONTROL OVER ERRORS C INTRODUCED BY THE USERS CHOICE OF THIS SPACE. SEE INPUT PARAMETERS C BELOW. C C COLLOCATION OVER PIECEWISE POLYNOMIALS C C THE BASIC ASSUMPTION MADE IS THAT THE APPROXIMATE SOLUTION C SATISFIES C NCPTS C U(T,X) = SUM C(I,T) * BF(I,X) C I=1 C C WHERE THE UNKNOWN COEFFICIENTS C DEPEND ONLY ON THE TIME T AND C THE KNOWN BASIS FUNCTIONS DEPEND ONLY ON X (WE HAVE ASSUMED THAT C NPDE = 1 FOR CONVENIENCE). SO, AT ANY GIVEN TIME T THE APPROX- C IMATE SOLUTION IS A PIECEWISE POLYNOMIAL IN THE USER CHOSEN SPACE. C THE SEMI-DISCRETE EQUATIONS (ACTUALLY ORDINARY DIFFERENTIAL C EQUATIONS) WHICH DETERMINE THE COEFFICIENTS C ARE OBTAINED BY C REQUIRING THAT THE ABOVE APPROXIMATE U(T,X) SATISFY THE PDE AND C BOUNDARY CONDITIONS EXACTLY AT A SET OF NCPTS COLLOCATION POINTS C (SEE COLPNT). THUS, PDECOL ACTUALLY COMPUTES THE BASIS FUNCTION C COEFFICIENTS RATHER THAN SPECIFIC APPROXIMATE SOLUTION VALUES. C C----------------------------------------------------------------------- C C THE MODIFICATION MADE TO THE ORIGINAL PDECOL PROGRAM CONSISTS C ENTIRELY OF THE REMOVAL OF THE BAND LINEAR EQUATION SOLVER, DECB C (FOR THE DECOMPOSITION PHASE) AND SOLB (FOR THE SOLUTION PHASE), C AND ITS SUBSTITUTION BY CRDCMP, CRSLVE, THE CONSTITUENT PARTS OF C COLROW [5]. THE CHANGES NECESSITATED IN SUBROUTINE PDECOL BY THIS C SUBSTITUTION ARE LIMITED TO : C C (1) CHANGES TO THE DIMENSION OF PW OR WORK(IW17). C (2) RESULTING CHANGES IN THE CHECKS ON INPUT PARAMETERS. C (3) THE REMOVAL OF SOME VARIABLES USED SOLELY TO PROVIDE C INFORMATION TO DECB, SOLB, AND THE INSERTION OF NEW C VARIABLES TO PASS INFORMATION TO CRDCMP, CRSLVE. C C----------------------------------------------------------------------- C C REFERENCES C C 1. MADSEN, N.K. AND R.F. SINCOVEC, PDECOL - COLLOCATION SOFTWARE C FOR PARTIAL DIFFERENTIAL EQUATIONS, ACM-TOMS, VOL. 5, NO. 3, C 1979, PP. 326-351. C 2. SINCOVEC, R.F. AND N.K. MADSEN, SOFTWARE FOR NONLINEAR PARTIAL C DIFFERENTIAL EQUATIONS, ACM-TOMS, VOL. 1, NO. 3, C SEPTEMBER 1975, PP. 232-260. C 3. HINDMARSH, A.C., PRELIMINARY DOCUMENTATION OF GEARIB.. SOLUTION C OF IMPLICIT SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS WITH C BANDED JACOBIANS, LAWRENCE LIVERMORE LAB, UCID-30130, FEBRUARY C 1976. C 4. DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J. C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472. C 5. DIAZ, J.C., G. FAIRWEATHER, AND P. KEAST, COLROW AND ARCECO : C FORTRAN PACKAGES FOR SOLVING CERTAIN ALMOST BLOCK DIAGONAL C LINEAR SYSTEMS BY MODIFIED ALTERNATE ROW AND COLUMN C ELIMINATION. ACM TRAN. MATH. SOFTW. VOL. 9, NO. 3, SEPT. 1983, C PP. 376,380. C C----------------------------------------------------------------------- C C USE OF PDECOL C C PDECOL IS CALLED ONCE FOR EACH DESIRED OUTPUT VALUE (TOUT) OF THE C TIME T, AND IT IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR, C STIFIB, WHICH ADVANCES THE TIME BY TAKING SINGLE STEPS UNTIL C T .GE. TOUT. INTERPOLATION TO THE EXACT TIME TOUT IS THEN DONE. C SEE TOUT BELOW. C C C SUMMARY OF SUGGESTED INPUT VALUES C C IT IS OF COURSE IMPOSSIBLE TO SUGGEST INPUT PARAMETER VALUES WHICH C ARE APPROPRIATE FOR ALL PROBLEMS. THE FOLLOWING SUGGESTIONS ARE TO C BE USED ONLY IF YOU HAVE NO IDEA OF BETTER VALUES FOR YOUR PROBLEM. C C DT = 1.E-10 C XBKPT = CHOOSE NINT+1 EQUALLY SPACED VALUES SUCH THAT XBKPT(1) = C XLEFT AND XBKPT(NINT+1) = XRIGHT. C EPS = 1.E-4 C NINT = ENOUGH SO THAT ANY FINE STRUCTURE OF THE PROBLEM MAY BE C RESOLVED. C KORD = 4 C NCC = 2 C MF = 22 C INDEX = 1 (ON FIRST CALL ONLY, THEN 0 THEREAFTER). C C C THE INPUT PARAMETERS ARE.. C T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE C (USED ONLY ON FIRST CALL). C TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. SINCE C THE PACKAGE CHOOSES ITS OWN TIME STEP SIZES, THE C INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT C AND THE PACKAGE WILL INTERPOLATE TO T = TOUT. C DT = THE INITIAL STEP SIZE IN T, IF INDEX = 1, OR, THE C MAXIMUM STEP SIZE ALLOWED (MUST BE .GT. 0), IF INDEX = 3. C USED FOR INPUT ONLY WHEN INDEX = 1 OR 3. SEE BELOW. C XBKPT = THE ARRAY OF PIECEWISE POLYNOMIAL BREAKPOINTS. C THE NINT+1 VALUES MUST BE STRICTLY INCREASING WITH C XBKPT(1) = XLEFT AND XBKPT(NINT+1) = XRIGHT (USED ONLY C ON FIRST CALL). C EPS = THE RELATIVE TIME ERROR BOUND (USED ONLY ON THE C FIRST CALL, UNLESS INDEX = 4). SINGLE STEP ERROR C ESTIMATES DIVIDED BY CMAX(I) WILL BE KEPT LESS THAN C EPS IN ROOT-MEAN-SQUARE NORM. THE VECTOR CMAX OF WEIGHTS C IS COMPUTED IN PDECOL. INITIALLY CMAX(I) IS SET TO C ABS(C(I)), WITH A DEFAULT VALUE OF 1 IF ABS(C(I)) .LT. 1. C THEREAFTER, CMAX(I) IS THE LARGEST VALUE C OF ABS(C(I)) SEEN SO FAR, OR THE INITIAL CMAX(I) IF C THAT IS LARGER. TO ALTER EITHER OF THESE, CHANGE THE C APPROPRIATE STATEMENTS IN THE DO-LOOPS ENDING AT C STATEMENTS 50 AND 130 BELOW. THE USER SHOULD EXERCISE C SOME DISCRETION IN CHOOSING EPS. IN GENERAL, THE C OVERALL RUNNING TIME FOR A PROBLEM WILL BE GREATER IF C EPS IS CHOSEN SMALLER. THERE IS USUALLY LITTLE REASON TO C CHOOSE EPS MUCH SMALLER THAN THE ERRORS WHICH ARE BEING C INTRODUCED BY THE USERS CHOICE OF THE POLYNOMIAL SPACE. C SEE RELATED COMMENTS CONCERNING CMAX BELOW STATEMENT 40. C NINT = THE NUMBER OF SUBINTERVALS INTO WHICH THE SPATIAL DOMAIN C (XLEFT,XRIGHT) IS TO BE DIVIDED (MUST BE .GE. 1) C (USED ONLY ON FIRST CALL). C KORD = THE ORDER OF THE PIECEWISE POLYNOMIAL SPACE TO BE USED. C ITS VALUE MUST BE GREATER THAN 2 AND LESS THAN 21. FOR C FIRST ATTEMPTS WE RECOMMEND KORD = 4. IF THE SOLUTION C IS SMOOTH AND MUCH ACCURACY IS DESIRED, HIGHER VALUES C MAY PROVE TO BE MORE EFFICIENT. WE HAVE SELDOM USED C VALUES OF KORD IN EXCESS OF 8 OR 9, THOUGH THEY ARE C AVAILABLE FOR USE IN PDECOL (USED ONLY ON FIRST CALL). C NCC = THE NUMBER OF CONTINUITY CONDITIONS TO BE IMPOSED ON THE C APPROXIMATE SOLUTION AT THE BREAKPOINTS IN XBKPT. C NCC MUST BE GREATER THAN 1 AND LESS THAN KORD. WE C RECOMMEND THE USE OF NCC = 2 (WITH NOGAUS = 0, SEE C BELOW), SINCE THEORY PREDICTS THAT DRAMATICALLY MORE C ACCURATE RESULTS CAN OFTEN BE OBTAINED USING THIS CHOICE C (USED ONLY ON FIRST CALL). C NPDE = THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS IN THE SYSTEM C TO BE SOLVED (USED ONLY ON FIRST CALL). C MF = THE METHOD FLAG (USED ONLY ON FIRST CALL, UNLESS C INDEX = 4). ALLOWED VALUES ARE 11, 12, 21, 22. C FOR FIRST ATTEMPTS WE RECOMMEND THE USE OF MF = 22. C MF HAS TWO DECIMAL DIGITS, METH AND MITER C (MF = 10*METH + MITER). C METH IS THE BASIC METHOD INDICATOR.. C METH = 1 MEANS THE ADAMS METHODS (GENERALIZATIONS OF C CRANK-NICOLSON). C METH = 2 MEANS THE BACKWARD DIFFERENTIATION C FORMULAS (BDF), OR STIFF METHODS OF GEAR. C MITER IS THE ITERATION METHOD INDICATOR C AND DETERMINES HOW THE JACOBIAN MATRIX IS C TO BE COMPUTED.. C MITER = 1 MEANS CHORD METHOD WITH ANALYTIC JACOBIAN. C FOR THIS USER SUPPLIES SUBROUTINE DERIVF. C SEE DESCRIPTION ABOVE. C MITER = 2 MEANS CHORD METHOD WITH JACOBIAN CALCULATED C INTERNALLY BY FINITE DIFFERENCES. SEE C SUBROUTINES PSETIB AND DIFFF. C INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL, C WITH THE FOLLOWING VALUES AND MEANINGS.. C 1 THIS IS THE FIRST CALL FOR THIS PROBLEM. C 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM, C AND INTEGRATION IS TO CONTINUE. C 2 SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT C EXACTLY (NO INTERPOLATION IS DONE). SEE NOTE C BELOW. ASSUMES TOUT .GE. THE CURRENT T. C IF TOUT IS .LT. THE CURRENT TIME, THEN TOUT IS C RESET TO THE CURRENT TIME AND CONTROL IS C RETURNED TO THE USER. A CALL TO VALUES WILL C PRODUCE ANSWERS FOR THE NEW VALUE OF TOUT. C 3 SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING C PROGRAM AFTER ONE STEP. TOUT IS IGNORED AND C DT MUST BE SET .GT. 0 TO A MAXIMUM ALLOWED C DT VALUE. SEE ABOVE. C 4 THIS IS NOT THE FIRST CALL FOR THE PROBLEM, C AND THE USER HAS RESET EPS AND/OR MF. C SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0, C IT NEED NOT BE RESET FOR NORMAL CONTINUATION. C C NOTE.. THE PACKAGE MUST HAVE TAKEN AT LEAST ONE SUCCESSFUL TIME C STEP BEFORE A CALL WITH INDEX = 2 OR 4 IS ALLOWED. C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN. C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. C A CHANGE OF PARAMETERS WITH INDEX = 4 CAN BE MADE AFTER C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN PROVIDED AT LEAST C ONE SUCCESSFUL TIME STEP HAS BEEN MADE. C C WORK = FLOATING POINT WORKING ARRAY FOR PDECOL. WE RECOMMEND C THAT IT BE INITIALIZED TO ZERO BEFORE THE FIRST CALL C TO PDECOL. ITS TOTAL LENGTH MUST BE AT LEAST C C KORD + 4*NPDE + NCPTS*(3*KORD + 2) + C NPDE*NCPTS*(MAXDER + 6) + (NPDE**2)*(13+KORD*(KORD-NCC)*NINT) C C WHERE MAXDER IS DEFINED BELOW (SEE STORAGE C ALLOCATION). C C IWORK = INTEGER WORKING ARRAY FOR PDECOL. THE FIRST TWO C LOCATIONS MUST BE DEFINED AS FOLLOWS... C IWORK(1) = LENGTH OF USERS ARRAY WORK C IWORK(2) = LENGTH OF USERS ARRAY IWORK C THE TOTAL LENGTH OF IWORK MUST BE AT LEAST C NCPTS*(NPDE + 1). C OUTPUT C C THE SOLUTION VALUES ARE NOT RETURNED DIRECTLY TO THE USER BY PDECOL. C THE METHODS USED IN PDECOL COMPUTE BASIS FUNCTION COEFFICIENTS, SO C THE USER (AFTER A RETURN FROM PDECOL) MUST CALL THE PACKAGE ROUTINE C VALUES TO OBTAIN HIS APPROXIMATE SOLUTION VALUES AT ANY DESIRED SPACE C POINTS X AT THE TIME T = TOUT. SEE THE COMMENTS IN SUBROUTINE VALUES C FOR DETAILS ON HOW TO PROPERLY MAKE THE CALL. C C THE COMMON BLOCK /GEAR0/ CAN BE ACCESSED EXTERNALLY BY THE USER C IF DESIRED. IT CONTAINS THE STEP SIZE LAST USED (SUCCESSFULLY), C THE ORDER LAST USED (SUCCESSFULLY), THE NUMBER OF STEPS TAKEN C SO FAR, THE NUMBER OF RESIDUAL EVALUATIONS (RES CALLS) SO FAR, C AND THE NUMBER OF MATRIX EVALUATIONS (PSETIB CALLS) SO FAR. C DIFFUN CALLS ARE COUNTED IN WITH RESIDUAL EVALUATIONS. C C THE OUTPUT PARAMETERS ARE.. C DT = THE STEP SIZE USED LAST, WHETHER SUCCESSFULLY OR NOT. C TOUT = THE OUTPUT VALUE OF T. IF INTEGRATION WAS SUCCESSFUL, C AND THE INPUT VALUE OF INDEX WAS NOT 3, TOUT IS C UNCHANGED FROM ITS INPUT VALUE. OTHERWISE, TOUT C IS THE CURRENT VALUE OF T TO WHICH THE INTEGRATION C HAS BEEN COMPLETED. C INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS, C WITH THE FOLLOWING VALUES AND MEANINGS.. C 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND. C -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE C ERROR TEST EVEN AFTER REDUCING DT BY A FACTOR OF C 1.E10 FROM ITS INITIAL VALUE. C -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS C HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY C A TEST ON EPS. TOO MUCH ACCURACY HAS BEEN REQUESTED. C -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE C CORRECTOR CONVERGENCE EVEN AFTER REDUCING DT BY A C FACTOR OF 1.E10 FROM ITS INITIAL VALUE. C -4 SINGULAR MATRIX ENCOUNTERED. PROBABLY DUE TO STORAGE C OVERWRITES. C -5 INDEX WAS 4 ON INPUT, BUT THE DESIRED CHANGES OF C PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT C WAS NOT BEYOND T. INTERPOLATION TO T = TOUT WAS C PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN, C SIMPLY CALL AGAIN WITH INDEX = 4 AND A NEW TOUT. C -6 ILLEGAL INDEX VALUE. C -7 ILLEGAL EPS VALUE. C -8 AN ATTEMPT TO INTEGRATE IN THE WRONG DIRECTION. THE C SIGN OF DT IS WRONG RELATIVE TO T0 AND TOUT. C -9 DT .EQ. 0.0. C -10 ILLEGAL NINT VALUE. C -11 ILLEGAL KORD VALUE. C -12 ILLEGAL NCC VALUE. C -13 ILLEGAL NPDE VALUE. C -14 ILLEGAL MF VALUE. C -15 ILLEGAL BREAKPOINTS - NOT STRICTLY INCREASING. C -16 INSUFFICIENT STORAGE FOR WORK OR IWORK. C----------------------------------------------------------------------- C C C SUMMARY OF ALL PACKAGE ROUTINES C C PDECOL - STORAGE ALLOCATION, ERROR CHECKING, INITIALIZATION, REPEATED C CALLS TO STIFIB TO ADVANCE THE TIME. C C INTERP - INTERPOLATES COMPUTED BASIS FUNCTION COEFFICIENTS TO THE C DESIRED OUTPUT TIMES, TOUT, FOR USE BY VALUES. C C INITAL - INITIALIZATION, GENERATION AND STORAGE OF PIECEWISE POLY- C NOMIAL SPACE BASIS FUNCTION VALUES AND DERIVATIVES, DET- C ERMINES THE BASIS FUNCTION COEFFICINTS OF THE PIECEWISE C POLYNOMIALS WHICH INTERPOLATE THE USERS INITIAL CONDITIONS. C C COLPNT - GENERATION OF REQUIRED COLLOCATION POINTS. C C BSPLVD - B-SPLINE PACKAGE ROUTINES WHICH ALLOW FOR EVALUATION OF C BSPLVN ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE. C INTERV C C VALUES - GENERATION AT ANY POINT(S) OF VALUES OF THE COMPUTED C APPROXIMATE SOLUTION AND ITS DERIVATIVES WHICH ARE C PIECEWISE POLYNOMIALS. THE SUBROUTINE IS CALLED ONLY BY C THE USER. C C STIFIB - CORE INTEGRATOR, TAKES SINGLE TIME STEPS TO ADVANCE THE C TIME. ASSEMBLES AND SOLVES THE PROPER NONLINEAR EQUATIONS C WHICH ARE RELATED TO USE OF ADAMS OR GEAR TYPE INTEGRATION C FORMULAS. CHOOSES PROPER STEP SIZE AND INTEGRATION FORMULA C ORDER TO MAINTAIN A DESIRED ACCURACY. DESIGNED FOR ODE C PROBLEMS OF THE FORM A * (DY/DT) = G(T,Y). C C COSET - GENERATES INTEGRATION FORMULA AND ERROR CONTROL COEFFICIENTS. C C RES - COMPUTES RESIDUAL VECTORS USED IN SOLVING THE NONLINEAR C EQUATIONS BY A MODIFIED NEWTON METHOD. C C DIFFUN - COMPUTES A**-1 * G(T,Y) WHERE A AND G ARE AS ABOVE (STIFIB). C C ADDA - ADDS THE A MATRIX TO A GIVEN MATRIX IN ALMOST BLOCK DIAGONAL C C EVAL - EVALUATES THE COMPUTED PIECEWISE POLYNOMIAL SOLUTION AND C DERIVATIVES AT COLLOCATION POINTS. C C GFUN - EVALUATES THE FUNCTION G(T,Y) BY CALLING EVAL AND THE USER C SUBROUTINES F AND BNDRY. C C PSETIB - GENERATES PROPER JACOBIAN MATRICES REQUIRED BY THE MODIFIED C NEWTON METHOD. C C DIFFF - PERFORMS SAME ROLE AS THE USER ROUTINE DERIVF. COMPUTES C DERIVATIVE APPROXIMATIONS BY USE OF FINITE DIFFERENCES. C C CRDCMP - PERFORM AN LU DECOMPOSTION OF THE ALMOST BLOCK DIAGONAL C LINEAR SYSTEMS ARISING. C CRSLVE - PERFORM THE RIGHT-HAND SIDE MODIFICATIONS AND THE BACK- C SUBSTITUTIONS REQUIRED TO SOLVE THE ALMOST BLOCK DIAGONAL C SYSTEM AFTER IT HAS BEEN FACTORED. C C----------------------------------------------------------------------- C C C STORAGE ALLOCATION C C SINCE PDECOL IS A DYNAMICALLY DIMENSIONED PROGRAM, MOST OF ITS C WORKING STORAGE IS PROVIDED BY THE USER IN THE ARRAYS WORK AND IWORK. C THE FOLLOWING GIVES A LIST OF THE ARRAYS WHICH MAKE UP THE CONTENTS C WORK AND IWORK, THEIR LENGTHS, AND THEIR USES. WHEN MORE THAN ONE C NAME IS GIVEN, IT INDICATES THAT DIFFERENT NAMES ARE USED FOR THE C SAME ARRAY IN DIFFERENT PARTS OF THE PROGRAM. THE DIFFERENT NAMES C OCCUR BECAUSE PDECOL IS AN AMALGAMATION OF SEVERAL OTHER CODES C WRITTEN BY DIFFERENT PEOPLE AND WE HAVE TRIED TO LEAVE THE SEPARATE C PARTS AS UNCHANGED FROM THEIR ORIGINAL VERSIONS AS POSSIBLE. C C C NAMES LENGTH USE C --------- ------------ ------------------------------------- C C BC 4*NPDE**2 BOUNDARY CONDITION INFORMATION. C WORK C C A 3*KORD*NCPTS BASIS FUNCTION VALUES AT COLLOCATION POINT C WORK(IW1) C C XT NCPTS + KORD BREAKPOINT SEQUENCE FOR GENERATION OF BASI C WORK(IW2) FUNCTION VALUES. C C XC NCPTS COLLOCATION POINTS. C WORK(IW3) C C CMAX NPDE*NCPTS VALUES USED IN ESTIMATING TIME C YMAX INTEGRATION ERRORS. C WORK(IW4) C C ERROR NPDE*NCPTS TIME INTEGRATION ERRORS. C WORK(IW5) C C SAVE1 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION C WORK(IW6) METHOD. C C SAVE2 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION C WORK(IW7) METHOD. C C SAVE3 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION C WORK(IW8) METHOD. C C UVAL 3*NPDE WORKING STORAGE FOR VALUES OF U, UX, AND C WORK(IW9) USS AT ONE POINT. C C C NPDE*NCPTS* CURRENT BASIS FUNCTION COEFFICIENT VALUES C Y (MAXDER+1) AND THEIR SCALED TIME DERIVATIVES. C WORK(IW10) C C DFDU NPDE**2 WORKING STORAGE USED TO COMPUTE THE C WORK(IW11) JACOBIAN MATRIX. C C DFDUX NPDE**2 WORKING STORAGE USED TO COMPUTE THE C WORK(IW12) JACOBIAN MATRIX. C C DFDUXX NPDE**2 WORKING STORAGE USED TO COMPUTE THE C WORK(IW13) JACOBIAN MATRIX. C C DBDU NPDE**2 BOUNDARY CONDITION INFORMATION. C WORK(IW14) C C DBDUX NPDE**2 BOUNDARY CONDITION INFORMATION. C WORK(IW15) C C DZDT NPDE BOUNDARY CONDITION INFORMATION. C WORK(IW16) C C PW NPDE**2*(NCPTS*KORD- STORAGE AND PROCESSING C WORK(IW17) NCC*(KORD-2)) OF THE JACOBIAN MATRIX. C C ILEFT NCPTS POINTERS TO BREAKPOINT SEQUENCE FOR C IWORK GENERATION OF BASIS FUNCTION VALUES. C C IPIV NPDE*NCPTS PIVOT INFORMATION FOR THE LU DECOMPOSED C IWORK(IW18) JACOBIAN MATRIX PW. C C WHERE... C C NCPTS = KORD*NINT - NCC*(NINT-1) C MAXDER = 5 UNLESS OTHERWISE SET BY THE USER INTO /OPTION/. C C IN THE ORIGINAL VERSION OF PDECOL A SWITCH CALLED IQUAD C (IQUAD = 1 IF KORD = 3 AND A NULL BOUNDARY CONDITION EXISTS C IQUAD = 0 OTHERWISE) WAS USED TO HANDLE CERTAIN PROBLEMS. C THE NEW VERSION OF THIS PACKAGE IS NOT INTENDED FOR USE ON C SUCH PROBLEMS AND ALL CODE RELATED TO IQUAD HAS THUS BEEN DELETED C C THE COMMON BLOCK /OPTION/ CONTAINS THE VARIABLES NOGAUS AND MAXDER. C NOGAUS IS SET .EQ. 0 IN THE BLOCK DATA. IT CAN BE CHANGED TO BE C SET .EQ. 1 IF THE GAUSS-LEGENDRE COLLOCATION POINTS ARE NOT C DESIRED WHEN NCC = 2 (SEE ABOVE AND COLPNT). MAXDER IS SET C .EQ. 5 IN THE BLOCK DATA AND ITS VALUE REPRESENTS THE C MAXIMUM ORDER OF TIME INTEGRATION FORMULA ALLOWED. ITS VALUE C AFFECTS THE STORAGE REQUIRED IN WORK AND MAY BE CHANGED IF C DESIRED. SEE COSET FOR RESTRICTIONS. THESE CHANGES MAY BE MADE BY C THE USER BY ACCESSING /OPTION/ IN HIS CALLING PROGRAM (BEFORE THE C FIRST CALL TO PDECOL) OR BY CHANGING THE DATA STATEMENT IN C THE BLOCK DATA. C C----------------------------------------------------------------------- C C C COMMUNICATION C C EACH SUBROUTINE IN THE PACKAGE CONTAINS A COMMUNICATION SUMMARY C AS INDICATED BELOW. C C PACKAGE ROUTINES CALLED.. EVAL,INITAL,INTERP,STIFIB C USER ROUTINES CALLED.. BNDRY C CALLED BY.. USERS MAIN PROGRAM C FORTRAN FUNCTIONS USED.. ABS,MAX,FLOAT,SQRT C----------------------------------------------------------------------- DIMENSION WORK(1),IWORK(1),XBKPT(1) INTEGER BOTPOS,SIZE COMMON /GEAR0/ DTUSED,NQUSED,NSTEP,NFE,NJE COMMON /GEAR1/ T,DTC,DTMN,DTMX,EPSC,UROUND,N,MFC,KFLAG,JSTART COMMON /GEAR9/ EPSJ,R0,IDUM(6) COMMON /OPTION/ NOGAUS,MAXDER COMMON /SIZES/ NIN,KOR,NC,NPD,NCPTS,NEQN,IFLAG COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10, * IW11,IW12,IW13,IW14,IW15,IW16,IW17,IW18 COMMON /IOUNIT/ LOUT IF (INDEX .EQ. 0) GO TO 60 IF (INDEX .EQ. 2) GO TO 70 IF (INDEX .EQ. 4) GO TO 80 IF (INDEX .EQ. 3) GO TO 90 C----------------------------------------------------------------------- C SEVERAL CHECKS ARE MADE HERE TO DETERMINE IF THE INPUT PARAMETERS C HAVE LEGAL VALUES. ERROR CHECKS ARE MADE ON INDEX, EPS, (T0-TOUT)*DT, C DT, NINT, KORD, NCC, NPDE, MF, WHETHER THE BREAKPOINT SEQUENCE IS C STRICTLY INCREASING, AND WHETHER THERE IS SUFFICIENT STORAGE C PROVIDED FOR WORK AND IWORK. PROBLEM DEPENDENT PARAMETERS ARE C CALCULATED AND PLACED IN COMMON. C----------------------------------------------------------------------- IERID = -6 IF (INDEX .NE. 1) GO TO 320 IERID = IERID - 1 IF (EPS .LE. 0.) GO TO 320 IERID = IERID - 1 IF ((T0-TOUT)*DT .GT. 0.) GO TO 320 IERID = IERID - 1 IF (DT .EQ. 0.0) GO TO 320 IERID = IERID - 1 NIN = NINT IF (NIN .LT. 1) GO TO 320 IERID = IERID - 1 KOR = KORD IF (KOR .LT. 3 .OR. KOR .GT. 20) GO TO 320 IERID = IERID - 1 NC = NCC IF (NCC .LT. 2 .OR. NCC .GE. KOR) GO TO 320 IERID = IERID - 1 NPD = NPDE NPDE2 = NPD*NPD IF (NPDE .LT. 1) GO TO 320 IERID = IERID - 1 IF (MF.NE.22.AND.MF.NE.21.AND.MF.NE.12.AND.MF.NE.11) GO TO 320 IERID = IERID - 1 DO 10 K=1,NIN IF(XBKPT(K) .GE. XBKPT(K+1)) GO TO 320 10 CONTINUE NCPTS = KOR + (NIN - 1) * (KOR - NCC) NEQN = NPDE * NCPTS IWSAVE = IWORK(1) IISAVE = IWORK(2) IW1 = 4*NPDE2 + 1 IW2 = IW1 + 3*KORD*NCPTS IW3 = IW2 + NCPTS + KORD IW4 = IW3 + NCPTS IW5 = IW4 + NEQN IW6 = IW5 + NEQN IW7 = IW6 + NEQN IW8 = IW7 + NEQN IW9 = IW8 + NEQN IW10 = IW9 + 3*NPDE IW11 = IW10 + NEQN*(MAXDER+1) IW12 = IW11 + NPDE2 IW13 = IW12 + NPDE2 IW14 = IW13 + NPDE2 IW15 = IW14 + NPDE2 IW16 = IW15 + NPDE2 IW17 = IW16 + NPDE IW18 = NCPTS + 1 IERID = IERID - 1 KRPT = KORD - NCC NROWS = KRPT*NPDE NELS = NROWS*NPDE*KORD NOVRLP = NCC*NPDE NCLBLK = KORD*NPDE NELTOP = NOVRLP*NPDE BOTPOS = NELTOP + NELS*NINT SIZE = BOTPOS + NELTOP IWSTOR = IW17 + SIZE - 1 IISTOR = IW18 + NEQN - 1 IF ( IWSAVE .LT. IWSTOR .OR. IISAVE .LT. IISTOR ) GO TO 335 C----------------------------------------------------------------------- C PERFORM INITIALIZATION TASKS. C----------------------------------------------------------------------- CALL INITAL(KOR,WORK(IW1),WORK(IW6),XBKPT,WORK(IW2),WORK(IW3), * WORK(IW17),IWORK(IW18),IWORK) IF (IFLAG .NE. 0) GO TO 280 C----------------------------------------------------------------------- C IF INITIAL VALUES OF CMAX OTHER THAN THOSE SET BELOW ARE DESIRED, C THEY SHOULD BE SET HERE. ALL CMAX(I) MUST BE POSITIVE. C HAVING PROPER VALUES OF CMAX FOR THE PROBLEM BEING SOLVED IS AS C IMPORTANT AS CHOOSING EPS (SEE ABOVE), SINCE ERRORS ARE C MEASURED RELATIVE TO CMAX. IF VALUES FOR DTMN OR DTMX, THE C BOUNDS ON ABS(DT), OTHER THAN THOSE BELOW ARE DESIRED, THEY C SHOULD BE SET BELOW. C----------------------------------------------------------------------- DO 50 I = 1,NEQN I1 = I - 1 WORK(IW4+I1) = ABS(WORK(IW6+I1)) IF (WORK(IW4+I1) .LT. 1.) WORK(IW4+I1) = 1. 50 WORK(IW10+I1) = WORK(IW6+I1) N = NEQN T = T0 DTC = DT DTMN = ABS(DT) DTUSED = 0. EPSC = EPS MFC = MF JSTART = 0 EPSJ = SQRT(UROUND) NHCUT = 0 KFLAG = 0 TOUTP = T0 IF ( T0 .EQ. TOUT ) GO TO 360 60 DTMX = ABS(TOUT-TOUTP)*10. GO TO 140 C 70 DTMX = ABS(TOUT-TOUTP)*10. IF ((T-TOUT)*DTC .GE. 0.) GO TO 340 GO TO 150 C 80 IF ((T-TOUT)*DTC .GE. 0.) GO TO 300 JSTART = -1 EPSC = EPS MFC = MF GO TO 100 C 90 DTMX = DT 100 IF ((T+DTC) .EQ. T) WRITE(LOUT,110) 110 FORMAT(36H WARNING.. T + DT = T ON NEXT STEP.) C----------------------------------------------------------------------- C TAKE A TIME STEP BY CALLING THE INTEGRATOR. C----------------------------------------------------------------------- CALL STIFIB (NEQN,WORK(IW10),WORK(IW4),WORK(IW5),WORK(IW6), * WORK(IW7),WORK(IW8),WORK(IW17),IWORK(IW18),WORK,IWORK) C KGO = 1 - KFLAG GO TO (120, 160, 220, 260, 280), KGO C KFLAG = 0, -1, -2, -3 -4 C 120 CONTINUE C----------------------------------------------------------------------- C NORMAL RETURN FROM INTEGRATOR. C C THE WEIGHTS CMAX(I) ARE UPDATED. IF DIFFERENT VALUES ARE DESIRED, C THEY SHOULD BE SET HERE. A TEST IS MADE FOR EPS BEING TOO SMALL C FOR THE MACHINE PRECISION. C C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY C STEP SHOULD BE INSERTED HERE. C C IF INDEX = 3, SAVE1 IS SET TO THE CURRENT C VALUES ON RETURN. C IF INDEX = 2, DT IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF C ERROR), AND THEN THE CURRENT C VALUES ARE PUT IN SAVE1 ON RETURN. C FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR C UNLESS TOUT HAS BEEN REACHED. THEN INTERPOLATED VALUES OF C ARE C COMPUTED AND STORED IN SAVE1 ON RETURN. C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 340 INSTEAD OF 360. C----------------------------------------------------------------------- D = 0. DO 130 I = 1,NEQN I1 = I - 1 AYI = ABS(WORK(IW10+I1)) WORK(IW4+I1) = MAX(WORK(IW4+I1), AYI) 130 D = D + (AYI/WORK(IW4+I1))**2 D = D*(UROUND/EPS)**2 IF (D .GT. FLOAT(NEQN)) GO TO 240 IF (INDEX .EQ. 3) GO TO 340 IF (INDEX .EQ. 2) GO TO 150 140 IF ((T-TOUT)*DTC .LT. 0.) GO TO 100 CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6)) GO TO 360 C 150 IF (((T+DTC)-TOUT)*DTC .LE. 0.) GO TO 100 IF (ABS(T-TOUT) .LE. 100.*UROUND*DTMX) GO TO 340 IF ((T-TOUT)*DTC .GE. 0.) GO TO 340 DTC = (TOUT - T)*(1. - 4.*UROUND) JSTART = -1 GO TO 100 C----------------------------------------------------------------------- C ON AN ERROR RETURN FROM INTEGRATOR, AN IMMEDIATE RETURN OCCURS IF C KFLAG = -2 OR -4, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE. C TO RECOVER, DT AND DTMN ARE REDUCED BY A FACTOR OF .1 UP TO 10 C TIMES BEFORE GIVING UP. C----------------------------------------------------------------------- 160 WRITE (LOUT,170) T 170 FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/ * 40H ERROR TEST FAILED WITH ABS(DT) = DTMIN/) 180 IF (NHCUT .EQ. 10) GO TO 200 NHCUT = NHCUT + 1 DTMN = .1*DTMN DTC = .1*DTC WRITE (LOUT,190) DTC 190 FORMAT(25H DT HAS BEEN REDUCED TO ,E16.8, * 26H AND STEP WILL BE RETRIED//) JSTART = -1 GO TO 100 C 200 WRITE (LOUT,210) 210 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//) GO TO 340 C 220 WRITE (LOUT,230) T,DTC 230 FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,6H DT =, * E16.8/52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//) GO TO 340 C 240 WRITE (LOUT,250) T 250 FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/ * 56H EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/) KFLAG = -2 GO TO 340 C 260 WRITE (LOUT,270) T 270 FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/ * 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/) GO TO 180 C 280 WRITE (LOUT,290) 290 FORMAT(//28H SINGULAR MATRIX ENCOUNTERED, * 35H PROBABLY DUE TO STORAGE OVERWRITES//) KFLAG = -4 GO TO 340 C 300 WRITE(LOUT,310) T,TOUT,DTC 310 FORMAT(//45H INDEX = -1 ON INPUT WITH (T-TOUT)*DT .GE. 0./ * 4H T =,E16.8,9H TOUT =,E16.8,8H DTC =,E16.8/ * 44H INTERPOLATION WAS DONE AS ON NORMAL RETURN./ * 41H DESIRED PARAMETER CHANGES WERE NOT MADE.) CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6)) INDEX = -5 RETURN C 320 WRITE(LOUT,330) IERID 330 FORMAT(//24H ILLEGAL INPUT...INDEX= ,I3//) INDEX = IERID RETURN C 335 WRITE(LOUT,336) IWSTOR,IWSAVE,IISTOR,IISAVE 336 FORMAT(//21H INSUFFICIENT STORAGE/24H WORK MUST BE OF LENGTH, * I10,5X,12HYOU PROVIDED,I10/24H IWORK MUST BE OF LENGTH,I10,5X, * 12HYOU PROVIDED,I10//) INDEX = IERID RETURN C 340 TOUT = T DO 350 I = 1,NEQN I1 = I - 1 350 WORK(IW6+I1) = WORK(IW10+I1) 360 INDEX = KFLAG TOUTP = TOUT DT = DTUSED IF (KFLAG .NE. 0) DT = DTC RETURN END * SUBROUTINE ADDA(PW,A,BC,NPDE) C----------------------------------------------------------------------- C SUBROUTINE ADDA ADDS ELEMENTS OF THE MATRIX A(Y,T) TO THE MATRIX C STORED IN PW. THE ELEMENTS OF A(Y,T) ARE STORED IN THE INPUT PARAMETER C MATRIX A(K,I,J). THIS ROUTINE EXTRACTS THE INFORMATION FROM THE LATTER C ADDS IT TO PW IN AN APPROPRIATE FASHION, CREATING THE EFFECT OF ADDING C A(Y,T) TO PW. ALSO, BOUNDARY CONDITION INFORMATION IS ADDED TO THE FI C AND LAST BLOCKS OF PW. C C DOCUMENTATION OF PARAMETERS C C PW - THE JACOBIAN MATRIX TO WHICH THE MATRIX A(Y,T) IS ADDED. C IT IS A VECTOR OF DIMENSION SIZE. PW CONTAINS, COLUMN-WISE, C THE BLOCK, TOP, ASSOCIATED WITH THE LEFT BOUNDARY CONDITIONS, C FOLLOWED BY THE NINT BLOCKS ASSOCIATED WITH COLLOCATION POINTS C ON EACH SUBINTERVAL, FOLLOWED BY THE BLOCK, BOT, ASSOCIATED WIT C THE RIGHT BOUNDARY CONDITIONS. THE TOP BLOCK CONSISTS OF A SING C ROW OF TWO SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE BOT BLOCK C THE SAME STRUCTURE. EACH OF THE OTHER BLOCKS IS A KRPT BY KORD C OF SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE (I,J)-TH SUCH SUBB C CORRESPONDS TO THE J-TH NON-ZERO BASIS FUNCTION OF THE SUBINTE C EVALUATED AT THE I-TH COLLOCATION POINT OF THAT SUBINTERVAL. C A - BASIS FUNCTION VALUES AT COLLOCATION POINTS. A IS USUALLY C A THREE DIMENSIONAL ARRAY OF SIZE (KORD,3,NCPTS), BUT IS C VIEWED HERE AS A ONE DIMENSIONAL ARRAY OF SIZE KORD*NCPTS*3. C A(K,J,I) CONTAINS THE VALUE OF THE (J-1)ST DERIVATIVE (J=1,2,3) C OF THE K-TH NON-ZERO BASIS FUNCTION (K=1,...,KORD) AT THE I-TH C COLLOCATION POINT. C BC - BOUNDARY CONDITION INFORMATION . BC IS USUALLY A THREE DIMENSIO C MATRIX OF SIZE (NPDE,NPDE,4) BUT HERE IT IS VIEWED AS A VECTOR C OF LENGTH 4*NPDE**2. C NPDE - THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. DIFFUN,PSETIB C FORTRAN FUNCTIONS USED.. NONE C * DIMENSION PW(1),A(1),BC(NPDE,NPDE,4) * INTEGER BOTPOS,SIZE * COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IDUM COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT * C DOCUMENTATION OF VARIABLES IN COMMON : C C NINT - NUMBER OF SUBINTERVALS. C KORD - NUMBER OF NON-ZERO BASIS FUNCTIONS PER SUBINTERVAL. C NCC - NUMBER OF CONTINUITY CONDITIONS AT EACH INTERNAL MESH POINT. C NPDE - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS. C NCPTS - NUMBER OF COLLOCATION POINTS. C NEQN - NUMBER OF DISCRETE EQUATIONS, NPDE * NCPTS C C NROWS - THE NUMBER OF ROWS IN EACH BLOCK TO BE PROVIDED TO CRDCMP. C NCLBLK - THE NUMBER OF COLUMNS IN EACH BLOCK TO BE PROVIDED TO CRDCM C NOVRLP - THE NUMBER OF COLUMNS OF OVERLAP BETWEEN THE BLOCKS PROVIDE C TO CRDCMP. C NELS - THE NUMBER OF ELEMENTS IN ONE COLLOCATION BLOCK OF PW. C NELTOP - THE NUMBER OF ELEMENTS IN THE TOP BLOCK STORED IN PW. C BOTPOS - THE POSITION IN PW WHERE STORAGE LOCATIONS FOR THE INFORMAT C FOR THE BOTTOM BLOCK BEGIN. C SIZE - THE LENGTH OF THE VECTOR PW. C KRPT - THE NUMBER OF COLLOCATION POINTS PER SUBINTERVAL. C----------------------------------------------------------------------- C C ADD THE BOUNDARY CONDITION INFORMATION, PROVIDED IN BC, TO THE TOP C AND BOTTOM BLOCK ROWS OF PW. C NPDE2 = NPDE*NPDE * DO 200 J = 1,NPDE * J1TOP = (J-1)*NPDE * C J1TOP+1 IS THE POSITION IN PW OF THE FIRST ROW OF THE FIRST NPDE C SUBBLOCK OF THE TOP BLOCK. * J2TOP = (J-1)*NPDE + NPDE2 * C J2TOP+1 IS THE POSITION IN PW OF THE FIRST ROW OF THE SECOND NPD C NPDE SUBBLOCK OF THE TOP BLOCK. * J1BOT = BOTPOS + J1TOP * C J1BOT+1 IS THE POSITION IN PW OF THE FIRST ROW OF THE FIRST NPDE C SUBBLOCK OF THE BOT BLOCK. * J2BOT = BOTPOS + J2TOP * C J2BOT+1 IS THE POSITION IN PW OF THE FIRST ROW OF THE SECOND NPD C NPDE SUBBLOCK OF THE BOT BLOCK. * DO 100 I = 1,NPDE * PW(J1TOP+I) = PW(J1TOP+I) + BC(I,J,1) PW(J2TOP+I) = PW(J2TOP+I) + BC(I,J,2) PW(J1BOT+I) = PW(J1BOT+I) + BC(I,J,3) PW(J2BOT+I) = PW(J2BOT+I) + BC(I,J,4) * 100 CONTINUE * 200 CONTINUE * C UPDATE THE REMAINING ROWS OF PW BY ADDING THE APPROPRIATE VALUES C IN A TO PW. * DO 600 I = 1,NINT * IPOS = (I-1)*NELS + NELTOP * C IPOS+1 IS THE POSITION IN PW OF THE START OF THE I-TH BLOCK * IC = (I-1)*KRPT * C IC+1 IS THE POSITION IN XC OF THE FIRST COLLOCATION POINT OF T C I-TH SUBINTERVAL * DO 500 IP = 1,KRPT * IPPOS = IPOS + (IP-1)*NPDE * C IPPOS+1 IS THE POSITION IN PW CORRESONDING TO THE ROW OF C SUBBLOCKS ASSOCIATED WITH THE IP-TH COLLOCATION POINT OF C I-TH SUBINTERVAL * IK = (IC+IP)*KORD*3 * C IK IS THE POSITION IN A(K,I,J) OF THE BASIS FUNCTION C INFORMATION CORRESPONDING TO THE IP-TH COLLOCATION POINT C THE I-TH SUBINTERVAL * DO 400 J = 1,KORD * JK = IK + J * C JK IS THE POSITION IN A(K,I,J) OF THE BASIS FUNCTI C INFORMATION CORRESPONDING TO THE J-TH NON-ZERO BAS C FUNCTION ASSOCIATED WITH THE IP-TH COLLOCATION POI C OF THE I-TH SUBINTERVAL * JPOS = IPPOS + (J-1)*NROWS*NPDE * C JPOS+1 IS THE POSITION IN PW OF THE J-TH SUBBLOCK C ROW OF SUBBLOCKS CORREPONDING TO THE IP-TH COLLOCA C POINT OF THE I-TH SUBINTERVAL * DO 300 JJ = 1,NPDE * JJPOS = JPOS + (JJ-1)*NROWS + JJ * C JJPOS IS THE POSITION IN PW OF THE JJ-TH D C ENTRY OF THE J-TH SUBBLOCK ABOVE. * C HERE WE ASSIGN THE VALUE OF THE J-TH NON-Z C BASIS FUNCTION EVALUATED AT THE IP-TH COLL C POINT OF THE I-TH SUBINTERVAL TO THE DIAGO C OF THE (IP,J)-TH SUBBLOCK OF THE COLLOCATI C BLOCK ASSOCIATED WITH THE I-TH SUBINTERVAL * PW(JJPOS) = PW(JJPOS) + A(JK) * 300 CONTINUE * 400 CONTINUE * 500 CONTINUE * 600 CONTINUE * RETURN END SUBROUTINE DIFFUN (N, T, Y, YDOT, IFLAG, PW, IPIV, WORK, IWORK) * * C----------------------------------------------------------------------- C STIFIB FINDS THE SOLUTION TO A SYSTEM OF O.D.E.'S OF THE FORM C A(Y,T) (DY/DT) = G(Y,T). THIS REQUIRES THE CALCULATION OF THE QUANTITY C A(Y,T)**-1 * G(Y,T). DIFFUN COMPUTES YDOT = A(Y,T)**-1 * G(Y,T) BY C SETTING UP AND SOLVING THE SYSTEM A(Y,T) * YDOT = G(Y,T) FOR YDOT. C C DOCUMENTATION OF PARAMETERS C C N - THE DIMENSION OF Y. (N = NEQN). C T - CURRENT VALUE OF THE INDEPENDENT VARIABLE T. C Y - THE VECTOR OF UNKNOWN BASIS FUNCTION COEFFICIENTS. C IT CONTAINS THE CURRENT APPROXIMATIONS OF THESE COEFFICIENTS. C THIS VECTOR IS ALSO KNOWN AS C IN OTHER PARTS OF THE PACKAGE. C YDOT - TO BE COMPUTED AS INDICATED ABOVE. C IFLAG- ERROR INDICATOR FROM THE MATRIX DECOMPOSITION ROUTINE. C IFLAG <> 0 IMPLIES FAILURE. C PW - COEFFICIENT MATRIX FOR DETERMINATION OF YDOT. C IPIV - PIVOTING INFORMATION FROM THE MATRIX DECOMPOSITION ROUTINE. C WORK - WORKING STORAGE VECTOR CONTAINING MANY OF THE IMPORTANT ARRAYS C USED IN THIS PACKAGE. COMPLETE DOCUMENTATION OF WORK IS GIVEN C IN THE PDECOL ROUTINE. HERE WE MENTION ONLY THOSE COMPONENTS C OF WORK THAT ARE USED IN THIS ROUTINE. EACH COMPONENT WILL BE C DESIGNATED BY ITS STARTING LOCATION IN WORK AND ALSO IT USUAL C NAME THROUGHOUT THE PACKAGE. C C WORK(IW1), A - BASIS FUNCTION VALUES AT COLLOCATION POINTS C A IS USUALLY A THREE DIMENSIONAL ARRAY OF S C (KORD,3,NCPTS), BUT IS VIEWED HERE AS A ONE C DIMENSIONAL ARRAY OF SIZE KORD*NCPTS*3. A(K C CONTAINS THE VALUE OF THE (J-1)ST DERIVATIV C (J=1,2,3) OF THE K-TH NON-ZERO BASIS FUNCTI C (K=1,...,KORD) AT THE I-TH COLLOCATION POIN C WORK(IW3), XC - COLLOCATION POINTS.THIS VECTOR IS OF LENGTH C WORK(IW9), UVAL - STORAGE FOR VALUES OF THE CURRENT SOLUTION C APPROXIMATION, U(X), AND ITS FIRST AND SECO C DERIVATIVES AT A SINGLE POINT, X. C THIS IS A TWO DIMENSIONAL ARRAY OF SIZE (NP C WORK(IW14), DBDU - BOUNDARY CONDITION INFORMATION. THIS TWO C DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AN C CONTAINS THE FIRST DERIVATIVE WITH RESPECT C OF THE LEFT HAND SIDE OF THE BOUNDARY CONDI C EQUATIONS. C C WORK(IW15), DBDUX - BOUNDARY CONDITION INFORMATION. THIS TWO C DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AN C CONTAINS THE FIRST DERIVATIVE WITH RESPECT C UX OF THE LEFT HAND SIDE OF THE BOUNDARY C CONDITION EQUATIONS. C WORK(IW16), DZDT - BOUNDARY CONDITION INFORMATION. THIS VECTOR C LENGTH NPDE AND CONTAINS THE DERIVATIVE WIT C RESPECT TO T OF THE RIGHT HAND SIDE OF THE C BOUNDARY CONDITION EQUATIONS. C C IWORK - INTEGER WORKING STORAGE VECTOR. IT IS DOCUMENTED IN PDECOL. C HERE WE MENTION ONLY THE VECTOR IT CONTAINS THAT IS RELEVANT C TO THIS ROUTINE. C C IWORK(1), ILEFT - VECTOR OF BREAKPOINT SEQUENCE INFORMATION. C THIS VECTOR IS OF LENGTH NCPTS. C C PACKAGE ROUTINES CALLED.. ADDA,DECB,GFUN,SOLB C USER ROUTINES CALLED.. NONE C CALLED BY.. STIFIB C FORTRAN FUNCTIONS USED.. NONE * DIMENSION Y(N),YDOT(N),PW(1),IPIV(1),WORK(1),IWORK(1) * INTEGER BOTPOS,SIZE,NELTP1,BOTPS1 * COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IDUM0 COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT COMMON /ISTART/ IW1,IDUM1,IW3,IDUM2(5),IW9,IDUM3(4),IW14, * IW15,IW16,IDUM4(2) * C DOCUMENTATION OF VARIABLES IN COMMON : C C NINT - NUMBER OF SUBINTERVALS. C KORD - NUMBER OF NON-ZERO BASIS FUNCTIONS PER SUBINTERVAL. C NCC - NUMBER OF CONTINUITY CONDITIONS AT EACH INTERNAL MESH POINT. C NPDE - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS. C NCPTS - NUMBER OF COLLOCATION POINTS. C NEQN - NUMBER OF DISCRETE EQUATIONS, NPDE * NCPTS C C NROWS - THE NUMBER OF ROWS IN EACH BLOCK TO BE PROVIDED TO CRDCMP. C NCLBLK - THE NUMBER OF COLUMNS IN EACH BLOCK TO BE PROVIDED TO CRDCMP. C NOVRLP - THE NUMBER OF COLUMNS OF OVERLAP BETWEEN THE BLOCKS PROVIDED C TO CRDCMP. C NELS - THE NUMBER OF ELEMENTS IN ONE COLLOCATION BLOCK OF PW. C NELTOP - THE NUMBER OF ELEMENTS IN THE TOP BLOCK STORED IN PW. C BOTPOS - THE POSITION IN PW WHERE THE STORAGE LOCATIONS FOR THE INFORM C FOR THE BOTTOM BLOCK BEGIN. C SIZE - THE LENGTH OF THE VECTOR PW. C KRPT - THE NUMBER OF COLLOCATION POINTS PER SUBINTERVAL. C C THE COMMON BLOCK /ISTART/ CONTAINS THE STARTING LOCATIONS, IN WORK C AND IWORK OF THE COMPONENT VECTORS CONTAINING THE ARRAYS OF C INFORMATION USED THROUGHTOUT THIS PACKAGE. * C----------------------------------------------------------------------- * * C EVALUATE THE FUNCTION G(Y,T). IT IS RETURNED IN THE VECTOR YDOT, C AND IS USED AS THE RIGHT HAND SIDE OF THE NEWTON SYSTEM. BOUNDARY C CONDITION INFORMATION IS ALSO UPDATED BY THIS CALL TO GFUN. * CALL GFUN (T, Y, YDOT, NPDE, NCPTS, WORK(IW1), WORK, WORK(IW14), * WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9), IWORK) * C THE FOLLOWING CALL TO GFUN IS EQUVIVALENT TO THE ABOVE CALL WITH THE C ACTUAL PARAMETER NAMES SUBSTITUED FOR THE WORK VECTOR ADDRESSES. C C CALL GFUN (T, Y, YDOT, NPDE, NCPTS, A, BC, DBDU, C * DBDUX, DZDT, XC, UVAL, ILEFT) * C INITIALIZE PW. * DO 10 I = 1,SIZE * PW(I) = 0.0E+0 * 10 CONTINUE * C ADD INFORMATION FROM THE MATRIX A(K,I,J) TO THE MATRIX PW. C PW WILL THEN CONTAIN THE MATRIX A(Y,T). * CALL ADDA (PW, WORK(IW1), WORK, NPDE) * C THE FOLLOWING CALL TO ADDA IS EQUVIVALENT TO THE ABOVE CALL WITH THE C ACTUAL PARAMETER NAMES SUBSTITUED FOR THE WORK VECTOR ADDRESSES. C C CALL ADDA (PW, A, BC, NPDE) * C FACTOR THE MATRIX PW. * BOTPS1 = BOTPOS + 1 NELTP1 = NELTOP + 1 IFLAG = 0 * CALL CRDCMP(NEQN,PW,NPDE,NOVRLP,PW(NELTP1),NROWS,NCLBLK,NINT, * PW(BOTPS1),NPDE,IPIV,IFLAG) * IF (IFLAG .NE. 0) GO TO 1021 * C SOLVE THE TRIANGULAR SYSTEMS RETURNING THE SOLUTION IN YDOT. * CALL CRSLVE(PW,NPDE,NOVRLP,PW(NELTP1),NROWS,NCLBLK,NINT, * PW(BOTPS1),NPDE,IPIV,YDOT,YDOT) * 1021 CONTINUE RETURN END SUBROUTINE INITAL(K,A,RHS,X,XT,XC,PW,IPIV,ILEFT) * * C----------------------------------------------------------------------- C INITAL IS CALLED ONLY ONCE BY PDECOL TO PERFORM INITIALIZATION TASKS. C THESE TASKS INCLUDE - 1) DEFINING THE PIECEWISE POLYNOMIAL SPACE C BREAKPOINT SEQUENCE, 2) CALLING THE SUBROUTINE COLPNT TO DEFINE THE C REQUIRED COLLOCATION POINTS, 3) DEFINING THE PIECEWISE POLYNOMIAL SPAC C BASIS FUNCTION VALUES (PLUS FIRST AND SECOND DERIVATIVE VALUES) AT C THE COLLOCATION POINTS, AND 4) DEFINING THE INITIAL BASIS FUNCTION C COEFFICIENTS WHICH DETERMINE THE PIECEWISE POLYNOMIAL WHICH C INTERPOLATES THE USER SUPPLIED INITIAL CONDITIONS AT THE C COLLOCATION POINTS. THIS LAST STEP INVOLVES THE CONSTRUCTION OF A LINE C SYSTEM OF THE FORM PW*C = RHS, WHERE C IS THE VECTOR OF C LENGTH NPDE*NCPTS CONTAINING THE UNKNOWN COEFFICIENTS TO BE DETERMINED C BY THIS ROUTINE. RHS CONSISTS OF NCPTS VECTORS EACH OF LENGTH NPDE. C THE I-TH SUCH VECTOR IS THE RIGHT HAND SIDE OF THE INITIAL CONDITIONS C EVALUATED AT THE I-TH COLLOCATION POINT. THE MATRIX PW CONSISTS OF AN C NCPTS BY NCPTS TWO DIMENSIONAL ARRAY OF BLOCKS, WITH EACH BLOCK BEING C SIZE NPDE BY NPDE. THE (I,J)-TH SUCH BLOCK CONSISTS OF THE IDENTITY MA C MULTIPLIED BY THE J-TH BASIS FUNCTION EVALUATED AT THE I-TH COLLOCATIO C POINT. BECAUSE OF SPECIAL PROPERTIES OF THE BASIS FUNCTIONS, ONLY KORD C OF THEM ARE NON-ZERO WHEN EVALUATED AT THE SET OF COLLOCATION POINTS C ASSOCIATED WITH A SINGLE SUBINTERVAL. THIS PROPERTY RESULTS IN THE PW C HAVING AN ALMOST BLOCK DIAGONAL STRUCTURE. C C DOCUMENTATION OF PARAMETERS C C K - EQUAL TO KORD, THE NUMBER OF NON-ZERO BASIS FUNCTIONS C PER SUBINTERVAL; USED FOR DIMENSIONING PURPOSES. C A - BASIS FUNCTION VALUES AT COLLOCATION POINTS. A IS C A THREE DIMENSIONAL ARRAY OF SIZE (KORD,3,NCPTS). C A(K,J,I) CONTAINS THE VALUE OF THE (J-1)ST DERIVATIVE (J=1,2,3 C OF THE K-TH NON-ZERO BASIS FUNCTION (K=1,...,KORD) AT THE I-TH C COLLOCATION POINT. C RHS - IN THIS ROUTINE THE RIGHT HAND SIDE OF THE INITAL CONDITIONS C EVALUATED AT EACH OF THE NCPTS COLLOCATION POINTS IS ASSIGNED C TO RHS. ON OUTPUT THIS VECTOR CONTAINS THE COEFFICIENTS WHICH C DETERMINE THE PIECE-WISE POLYNOMIAL INTERPOLATING THE INITIAL C CONDITIONS AT THE COLLOCATION POINTS. IT IS OF SIZE NPDE*NCPT C X - THE MESHPOINT SEQUENCE. IT IS OF SIZE NINT + 1. C XT - THE BREAKPOINT SEQUENCE. IT IS OF SIZE NCPTS + KORD. C XC - THE COLLOCATION POINT SEQUENCE. IT IS OF SIZE NCPTS. C PW - IN INITAL, THIS VECTOR, WHICH IS OF DIMENSION SIZE, WILL BE C ASSIGNED VALUES AS INDICATED ABOVE. PW CONTAINS, COLUMN-WISE, C THE BLOCK, TOP, ASSOCIATED WITH THE LEFT BOUNDARY CONDITIONS, C FOLLOWED BY THE NINT BLOCKS ASSOCIATED WITH COLLOCATION POINTS C ON EACH SUBINTERVAL, FOLLOWED BY THE BLOCK, BOT, ASSOCIATED WI C THE RIGHT BOUNDARY CONDITIONS. THE TOP BLOCK CONSISTS OF A SIN C ROW OF TWO SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE BOT BLOCK C THE SAME STRUCTURE. EACH OF THE OTHER BLOCKS IS A KRPT BY KORD C OF SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE (I,J)-TH SUCH SUB C CORRESPONDS TO THE J-TH NON-ZERO BASIS FUNCTION OF THE SUBINT C EVALUATED AT THE I-TH COLLOCATION POINT OF THAT SUBINTERVAL. C ON OUTPUT PW CONTAINS THE FACTORED COEFFICIENT MATRIX C USED TO DETERMINE THE COMPONENTS OF RHS. C IPIV - PIVOTING INFORMATION FROM THE FACTORIZATION OF PW. C IT IS OF SIZE NPDE*NCPTS. C ILEFT - BREAKPOINT INFORMATION. IT IS OF SIZE NCPTS. C C PACKAGE ROUTINES CALLED.. BSPLVD,COLPNT,CRDCMP,INTERV,CRSLVE C USER ROUTINES CALLED.. UINIT C CALLED BY.. PDECOL C FORTRAN FUNCTIONS USED.. MAX0,MIN0 * DIMENSION A(K,3,1),RHS(1),X(1),XT(1),XC(1),PW(1),IPIV(1),ILEFT(1) INTEGER BOTPOS,SIZE,BOTPS1 * * COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IFLAG COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT * C DOCUMENTATION OF VARIABLES IN COMMON : C C NINT - NUMBER OF SUBINTERVALS. C KORD - NUMBER OF NON-ZERO BASIS FUNCTIONS PER SUBINTERVAL. C NCC - NUMBER OF CONTINUITY CONDITIONS AT EACH INTERNAL MESH POINT. C NPDE - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS. C NCPTS - NUMBER OF COLLOCATION POINTS. C NEQN - NUMBER OF DISCRETE EQUATIONS, NPDE * NCPTS C IFLAG - ERROR INDICATOR FOR LINEAR SYSTEM FACTORIZATION ROUTINE. C C NROWS - THE NUMBER OF ROWS IN EACH BLOCK TO BE PROVIDED TO CRDCMP. C NCLBLK - THE NUMBER OF COLUMNS IN EACH BLOCK TO BE PROVIDED TO CRDCM C NOVRLP - THE NUMBER OF COLUMNS OF OVERLAP BETWEEN THE BLOCKS PROVIDE C TO CRDCMP. C NELS - THE NUMBER OF ELEMENTS IN ONE COLLOCATION BLOCK OF PW. C NELTOP - THE NUMBER OF ELEMENTS IN THE TOP BLOCK STORED IN PW. C BOTPOS - THE POSITION IN PW WHERE STORAGE LOCATIONS FOR THE INFORMAT C FOR THE BOTTOM BLOCK BEGIN. C SIZE - THE LENGTH OF THE VECTOR PW. C KRPT - THE NUMBER OF COLLOCATION POINTS PER SUBINTERVAL. C C----------------------------------------------------------------------- * C SET UP THE PIECEWISE POLYNOMIAL SPACE BREAKPOINT SEQUENCE. C FIRST SET UP THE SEQUENCE AT EACH ENDPOINT. DO 10 I=1,KORD * XT(NCPTS+I) = X(NINT+1) XT(I) = X(1) * 10 CONTINUE * C SET UP THE REMAINING POINTS OF THE SEQUENCE. * DO 30 I=2,NINT * IK = (I-2)*KRPT + KORD * C IK+1 IS THE POSITION IN XT OF THE START OF THE BREAKPOINT INFORM C ASSOCIATED WITH THE KRPT COLLOCATION POINTS LOCATED ON THE I-TH C SUBINTERVAL. * DO 20 J=1,KRPT * XT(IK+J) = X(I) * 20 CONTINUE * 30 CONTINUE * C SET UP COLLOCATION POINTS ARRAY XC. * CALL COLPNT(X, XC, XT) * C GENERATE THE ILEFT ARRAY. STORE THE BASIS FUNCTION VALUES IN THE C ARRAY A. THE ARRAY A IS DIMENSIONED A(KORD,3,NCPTS) AND A(K,J,I) C CONTAINS THE VALUE OF THE (J-1)-ST DERIVATIVE (J = 1,2,3) OF THE K-TH C NONZERO BASIS FUNCTION (K = 1, ... ,KORD) AT THE I-TH COLLOCATION C POINT (I = 1, ... ,NCPTS). SET UP RHS FOR INTERPOLATING THE INITIAL C CONDITIONS AT THE COLLOCATION POINTS. ASSIGN VALUES FROM THE ABOVE MA C A, TO THE MATRIX PW. C C INITIALIZE PW TO ZERO. C DO 1130 I = 1,SIZE * PW(I) = 0.0E+0 * 1130 CONTINUE * C FIRST, SET UP THE TOP BLOCK OF THE MATRIX, MAKING USE OF THE FACT C THAT ONLY THE FIRST B-SPLINE HAS A NON-ZERO VALUE AT THE LEFT END C POINT. ALSO, SET UP THE FIRST NPDE COMPONENTS OF THE RIGHT HAND C SIDE. * MFLAG = -2 CALL INTERV(XT,NCPTS,XC(1),ILEFT(1),MFLAG) * C HERE INTERV COMPUTES THE VALUE OF ILEFT(1). MFLAG IS INITIALIZED T C -2 TO INSTRUCT INTERV TO BEGIN ITS SEARCH FOR THE APPROPRIATE ILEF C WITH THE VALUE 1. * CALL BSPLVD(XT,KORD,XC(1),ILEFT(1),A(1,1,1),3) C C BSPLVD IS CALLED TO COMPUTE THE COMPONENTS OF A(K,I,J) ASSOCIATED C THE FIRST COLLOCATION POINT. C CALL UINIT(XC(1),RHS(1),NPDE) C C UINIT IS CALLED TO EVALUATE THE FIRST NPDE COMPONENTS OF THE RHS V C THE RIGHT HAND SIDE OF THE INITIAL CONDITIONS, EVALUATED AT THE FI C COLLOCATION POINT IS ASSIGNED TO THE RHS VECTOR. * DO 4040 N = 1,NPDE * NN = (N-1)*NPDE + N PW(NN) = A(1,1,1) * C HERE WE ASSIGN A(1,1,1), THE FIRST BASIS FUNCTION EVALUATED AT C FIRST COLLOCATION POINT, TO THE DIAGONAL OF THE FIRST (I.E TOP C BLOCK OF PW. THE INDEX NN MOVES THROUGH THE LOCATIONS IN PW C CORRESPONDING TO THIS DIAGONAL. * 4040 CONTINUE * C THE NINT BLOCKS OF THE MAIN MATRIX WILL NOW BE SET UP. EACH SUCH BLOCK C CONSISTS OF A KRPT BY KORD TWO DIMENSIONAL ARRAY OF SUBBLOCKS, EACH OF C SIZE NPDE BY NPDE. THE CORRESPONDING COMPONENTS OF THE RIGHT HAND SID C VECTOR ARE ALSO ASSIGNED. * * DO 1040 I = 1,NINT * IPOS = (I-1)*NELS + NELTOP * C IPOS+1 IS THE POSITION IN PW OF THE START OF THE I-TH BLOCK * IC = (I-1)*KRPT * C IC+1 IS THE POSITION IN XC OF THE FIRST COLLOCATION POINT OF THE C I-TH SUBINTERVAL * C FOR EACH COLLOCATION POINT OF THE I-TH SUBINTERVAL, C CONSTRUCT THE CORRESPONDING SUBBLOCK ROW OF THE I-TH BLOCK. * DO 1030 IP = 1,KRPT * IPPOS = IPOS + (IP-1)*NPDE * C IPPOS+1 IS THE POSITION IN PW CORRESONDING TO THE ROW OF SUBBL C ASSOCIATED WITH THE IP-TH COLLOCATION POINT OF THE I-TH SUBINT * ICC = IC + IP + 1 * C ICC IS THE ABSOLUTE POSITION IN XC OF THE IP-TH COLLKCATION PO C OF THE I-TH SUBINTERVAL * IN = (ICC-1)*NPDE * C IN+1 IS THE POSITION IN THE RHS VECTOR WHERE THE VALUES FOR TH C RIGHT HAND SIDE OF THE INITAL CONDITIONS, EVALUATED AT THE ICC C COLLOCATION POINT ARE STORED. * C COMPUTE INFORMATION FOR ICC-TH COLLOCATION POINT * CALL INTERV(XT,NCPTS,XC(ICC),ILEFT(ICC),MFLAG) CALL BSPLVD(XT,KORD,XC(ICC),ILEFT(ICC),A(1,1,ICC),3) CALL UINIT(XC(ICC),RHS(IN+1),NPDE) * C THE ABOVE THREE CALLS PERFORM THE SAME TASKS AS MENTIONED ABOV C EXCEPT THAT ALL EVAULATIONS ARE DONE AT THE ICC-TH COLLOCATION C POINT INSTEAD OF AT THE FIRST. * DO 1010 J = 1,KORD * C FOR J RANGING OVER THE KORD SUBBLOCKS OF THE IP-TH SUBBLOCK C OF THE I-TH SUBINTERVAL CONSTRUCT THE CORRESPONDING COMPONEN C PW. * JPOS = IPPOS + (J-1)*NROWS*NPDE * C JPOS+1 IS THE POSITION IN PW OF THE J-TH SUBBLOCK OF THE ROW C SUBBLOCKS CORREPONDING TO THE IP-TH COLLOCATION POINT OF THE C SUBINTERVAL * DO 1005 L = 1,NPDE * LPOS = JPOS + (L-1)*NROWS + L * C LPOS IS THE POSITION IN PW OF THE (L,L)-TH ENTRY OF THE C J-TH SUBBLOCK ABOVE. * PW(LPOS) = A(J,1,ICC) C print *,'LPOS = ',lpos,' J = ',j,' ICC = ',icc * C HERE WE ASSIGN THE J-TH NON-ZERO BASIS FUNCTION OF THE I C SUBINTERVAL, EVALUATED AT THE ICC-TH COLLOCATION POINT, C DIAGONAL OF THE (IP,J)-TH NPDE BY NPDE SUBBLOCK OF THE I C COLLOCATION BLOCK. * 1005 CONTINUE * 1010 CONTINUE * 1030 CONTINUE * 1040 CONTINUE * C NOW, SET UP THE BOTTOM BLOCK, USING THE FACT THAT ONLY THE LAST C B-SPLINE BASIS FUNCTION IS NON-ZERO AT THE RIGHT END POINT. C SIMULTANEOUSLY, SET UP THE CORRESPONDING PART OF THE RIGHT HAND SIDE. * CALL INTERV(XT,NCPTS,XC(NCPTS),ILEFT(NCPTS),MFLAG) CALL BSPLVD(XT,KORD,XC(NCPTS),ILEFT(NCPTS),A(1,1,NCPTS),3) IN = NEQN - NPDE +1 CALL UINIT(XC(NCPTS),RHS(IN),NPDE) * DO 3040 N = 1,NPDE * NN = ((N-1) + (NCC-1)*NPDE)*NPDE + N + BOTPOS PW(NN) = A(KORD,1,NCPTS) * 3040 CONTINUE * C LU DECOMPOSE THE MATRIX PW. * NELTP1 = NELTOP + 1 BOTPS1 = BOTPOS + 1 IFLAG = 0 * CALL CRDCMP(NEQN,PW,NPDE,NOVRLP,PW(NELTP1),NROWS,NCLBLK,NINT, * PW(BOTPS1),NPDE,IPIV,IFLAG) * * IF (IFLAG.NE.0) GO TO 1021 * C SOLVE THE LINEAR SYSTEM PW*Z = RHS. THIS GIVES THE BASIS FUNCTION C COEFFICIENTS FOR THE INITIAL CONDITIONS. THE SOLUTION, Z, IS RETURNED C IN THE VECTOR RHS. * * CALL CRSLVE(PW,NPDE,NOVRLP,PW(NELTP1),NROWS,NCLBLK,NINT, * PW(BOTPS1),NPDE,IPIV,RHS,RHS) * * 1021 CONTINUE * RETURN END SUBROUTINE PSETIB (C, PW, N0, CON, MITER, IFLAG, A, ILEFT, XC, * UVAL, SAVE2,IPIV,CMAX,DFDU,DFDUX,DFDUXX,DZDT,DBDU,DBDUX,BC, * NPDE) * * C----------------------------------------------------------------------- C STIFIB IS USED TO SOLVE A SYSTEM OF ODE'S OF THE FORM A(C,T)*DC/DT = G C WHEN NEWTON'S METHOD IS APPLIED TO SOLVE THE DISCRETE NON-LINEAR SYSTE C ARISING FROM THE ABOVE SYSTEM OF ODE'S, THE JACOBIAN MUST BE COMPUTED. C THE MATRIX THAT MUST BE COMPUTED HAS THE FORM A(C,T) + CON*(DG/DC), C WHERE THE CONSTANT, CON, IS DEFINED BELOW. PSETIB IS CALLED BY STIFIB C CONSTRUCT AND FACTOR THIS MATRIX, WHICH IS TO BE STORED IN PW. C DG/DC IS COMPUTED ACCORDING TO THE VALUE OF MITER DEFINED BELOW. C FINALLY, PW IS SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER C SOLUTION OF LINEAR SYSTEMS WITH PW AS COEFFICIENT MATRIX. C C DOCUMENTATION OF PARAMETERS C C C - THE VECTOR OF UNKNOWN COEFFICIENTS; THE VECTOR C C CONTAINS THE CURRENT APPROXIMATE VALUES OF THESE COEFFICIENTS. C THIS VECTOR IS OF SIZE NCPTS*NPDE. C C PW - IN PSETIB, THIS VECTOR, WHICH IS OF DIMENSION SIZE, WILL BE C ASSIGNED VALUES AS INDICATED ABOVE. PW CONTAINS, COLUMN-WISE, C THE BLOCK, TOP, ASSOCIATED WITH THE LEFT BOUNDARY CONDITIONS, C FOLLOWED BY THE NINT BLOCKS ASSOCIATED WITH COLLOCATION POINTS C ON EACH SUBINTERVAL, FOLLOWED BY THE BLOCK, BOT, ASSOCIATED WI C THE RIGHT BOUNDARY CONDITIONS. THE TOP BLOCK CONSISTS OF A SIN C ROW OF TWO SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE BOT BLOCK C THE SAME STRUCTURE. EACH OF THE OTHER BLOCKS IS A KRPT BY KORD C OF SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE (I,J)-TH SUCH SUB C CORRESPONDS TO THE J-TH NON-ZERO BASIS FUNCTION OF THE SUBINT C EVALUATED AT THE I-TH COLLOCATION POINT OF THAT SUBINTERVAL. C ON OUTPUT PW CONTAINS THE FACTORED COEFFICIENT MATRIX. C C NO - EQUAL TO NEQN. USED TO PROVIDE RUN-TIME C DIMENSIONING FOR SEVERAL ARRAYS. C C CON - THIS CONSTANT CONTAINS THE PRODUCT OF TWO VALUES MENTIONED C ABOVE IN REFERENCE TO THE DEFINITION OF PW. IT EQUALS C -H * EL(1), WHERE EL(1) IS A COEFFICIENT FROM THE NUMERICAL C FORMULA USED TO DISCRETIZE THE SYSTEM OF ODE'S. C C MITER - THIS FLAG INDICATES WHETHER OR NOT AN ANALYTIC JACOBIAN IS C TO BE PROVIDED BY THE USER; MITER = 1 IMPLIES AN ANALYTIC C JACOBIAN WILL BE PROVIDED THROUGH THE ROUTINE DERIVF; C MITER = 2 IMPLIES THAT A DIVIDED-DIFFERENCE APPROXIMATION C TO THE JACOBIAN MUST BE COMPUTED. C C IFLAG - ERROR INDICATOR FROM THE MATRIX DECOMPOSITION ROUTINE. C IFLAG <> 0 IMPLIES A FAILURE. C C A - BASIS FUNCTION VALUES AT COLLOCATION POINTS. A IS USUALLY C A THREE DIMENSIONAL ARRAY OF SIZE (KORD,3,NCPTS), BUT IS C VIEWED HERE AS A ONE DIMENSIONAL ARRAY OF SIZE KORD*NCPTS*3. C A(K,J,I) CONTAINS THE VALUE OF THE (J-1)ST DERIVATIVE (J=1,2,3 C OF THE K-TH NON-ZERO BASIS FUNCTION (K=1,...,KORD) AT THE I-TH C COLLOCATION POINT. C C ILEFT - POINTERS TO BREAKPOINT SEQUENCE. THIS VECTOR IS OF SIZE NCPTS. C C XC - COLLOCATION POINTS. THIS VECTOR IS OF LENGTH NCPTS. C C UVAL - WORKING STORAGE FOR VALUES OF U, UX, AND UXX AT ONE POINT. C THIS IS A TWO DIMENSIONAL ARRAY OF SIZE (NPDE,3). C C SAVE2 - WORKING STORAGE FOR THE TIME INTEGRATION METHOD. C THIS VECTOR IS OF LENGTH NPDE*NCPTS. IT CONTAINS NCPTS VECTORS C EACH OF LENGTH NPDE. EACH OF THESE VECTORS CORRESPONDS TO AN C EVALUATION OF THE RIGHT HAND SIDE OF THE SYSTEM OF PDE'S AT C A COLLOCATION POINT. C C IPIV - PIVOT INFORMATION FOR THE LU DECOMPOSED JACOBIAN MATRIX PW. C THIS VECTOR IS OF LENGTH NPDE*NCPTS. C C CMAX - VALUES USED IN ESTIMATING TIME INTEGRATION ERRORS. C THIS INFORMATION IS USED BY THE ROUTINE DIFFF WHEN A NUMERICAL C APPROXIMATION TO THE JACOBIAN MUST BE COMPUTED. C C DFDU - WORKING STORAGE USED TO COMPUTE THE JACOBIAN MATRIX. C THIS TWO DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AND CONTAINS C THE DERIVATIVE WITH RESPECT TO U OF THE RIGHT HAND SIDE OF THE C SYSTEM OF PDE'S, EVALUATED AT A COLLOCATION POINT. C C DFDUX - WORKING STORAGE USED TO COMPUTE THE JACOBIAN MATRIX. C THIS TWO DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AND CONTAINS C THE DERIVATIVE WITH RESPECT TO UX OF THE RIGHT HAND SIDE OF TH C SYSTEM OF PDE'S, EVALUATED AT A COLLOCATION POINT. C C DFDUXX- WORKING STORAGE USED TO COMPUTE THE JACOBIAN MATRIX. C THIS TWO DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AND CONTAINS C THE DERIVATIVE WITH RESPECT TO UXX OF THE RIGHT HAND SIDE OF T C SYSTEM OF PDE'S, EVALUATED AT A COLLOCATION POINT. C C DZDT - BOUNDARY CONDITION INFORMATION. THIS VECTOR IS OF LENGTH NPDE C AND CONTAINS THE DEPIVATIVE WITH RESPECT TO T OF THE RIGHT HAN C SIDE OF THE BOUNDARY CONDITION EQUATIONS. C C DBDU - BOUNDARY CONDITION INFORMATION. THIS TWO DIMENSIONAL ARRAY IS C SIZE (NPDE,NPDE) AND CONTAINS THE FIRST DERIVATIVE WITH RESPEC C OF THE LEFT HAND SIDE OF THE BOUNDARY CONDITION EQUATIONS. C C DBDUX - BOUNDARY CONDITION INFORMATION. THIS TWO DIMENSIONAL ARRAY IS C SIZE (NPDE,NPDE) AND CONTAINS THE FIRST DERIVATIVE WITH RESPEC C OF THE LEFT HAND SIDE OF THE BOUNDARY CONDITION EQUATIONS. C C BC - BOUNDARY CONDITION INFORMATION. BC IS USUALLY A THREE DIMENSIO C MATRIX OF SIZE (NPDE,NPDE,4) BUT HERE IT IS VIEWED AS A VECTOR C LENGTH 4*NPDE**2. C C C NPDE - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS C C PACKAGE ROUTINES CALLED.. ADDA,CRDCMP,DIFFF,EVAL,GFUN C USER ROUTINES CALLED.. BNDRY,DERIVF C CALLED BY.. STIFIB C FORTRAN FUNCTIONS USED.. ABS,FLOAT,MAX,MIN,SQRT C INTEGER BOTPOS,SIZE,NELTP1,BOTPS1 DIMENSION PW(SIZE),C(1),CMAX(1) DIMENSION A(1),ILEFT(1),BC(1),XC(1),UVAL(NPDE,3),SAVE2(1),IPIV(1) DIMENSION DFDU(NPDE,NPDE),DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE) DIMENSION DZDT(NPDE),DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE) * COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IDUM COMMON /GEAR1/ T,H,DUM1(3),UROUND,N,IDUM1(3) COMMON /GEAR9/ DUM2,R0,IDUM2(6) COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT * C DOCUMENTATION OF VARIABLES IN COMMON : C C NINT - NUMBER OF SUBINTERVALS. C KORD - NUMBER OF NON-ZERO BASIS FUNCTIONS PER SUBINTERVAL. C NCC - NUMBER OF CONTINUITY CONDITIONS AT EACH INTERNAL MESH POINT. C NPDE - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS. C NCPTS - NUMBER OF COLLOCATION POINTS. C NEQN - NUMBER OF DISCRETE EQUATIONS, NPDE * NCPTS C C T - THE INDEPENDENT VARIABLE REPRESENTING THE TIME DIMENSION. C H - THE CURRENT STEPSIZE BEING USED IN THE STIFIB ROUTINE. C UROUND - THE UNIT ROUND-OFF OF THE MACHINE. C N - THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS C IN THE SYSTEMS TO BE SOLVED BY STIFIB. C C R0 - SCALING FACTOR - USED BY THE DIFFF ROUTINE AND CALCULATED IN C PSETIB. C C NROWS - THE NUMBER OF ROWS IN EACH BLOCK TO BE PROVIDED TO CRDCMP. C NCLBLK - THE NUMBER OF COLUMNS IN EACH BLOCK TO BE PROVIDED TO CRDCMP. C NOVRLP - THE NUMBER OF COLUMNS OF OVERLAP BETWEEN THE BLOCKS PROVIDED C TO CRDCMP. C NELS - THE NUMBER OF ELEMENTS IN ONE COLLOCATION BLOCK OF PW. C NELTOP - THE NUMBER OF ELEMENTS IN THE TOP BLOCK STORED IN PW. C BOTPOS - THE POSITION IN PW WHERE THE INFORMATION FOR THE BOTTOM C BLOCK IS STORED. C SIZE - THE LENGTH OF THE VECTOR PW. C KRPT - THE NUMBER OF COLLOCATION POINTS PER SUBINTERVAL. C C C----------------------------------------------------------------------- * C INITIALIZE PW. * DO 10 I=1,SIZE * PW(I) = 0.0E+0 * 10 CONTINUE * C IF THE USER DOES NOT PROVIDE THE ANALYTIC JACOBIAN, (I.E. C MITER IS NOT EQUAL TO 1), THEN, BY CALLING GFUN, COMPUTE C THE REQUIRED INFORMATION TO BE USED LATER, DURING THE CALL C TO DIFFF, IN THE CALCULATION OF AN APPROXIMATION TO THE JACOBIAN. * IF ( MITER .EQ. 1 ) GO TO 25 * CALL GFUN (T, C, SAVE2, NPDE, NCPTS,A,BC,DBDU,DBDUX,DZDT,XC, * UVAL,ILEFT) * C COMPUTE R0 FOR USE IN THE DIFFF ROUTINE. * D = 0.0E+0 * DO 20 I = 1,N * D = D + SAVE2(I)**2 * 20 CONTINUE * R0 = ABS(H)* SQRT(D/FLOAT(N0))*1.0D+3*UROUND * 25 CONTINUE * C FOR EACH OF THE NINT SUBINTERVALS COMPUTE THE CORRESPONDING C BLOCK OF THE COEFFICIENT MATRIX PW. * DO 300 I=1,NINT * IPOS = (I-1)*NELS + NELTOP * C IPOS+1 IS THE POSITION IN PW OF THE START OF THE I-TH BLOCK * IC = (I-1)*KRPT * C IC+1 IS THE POSITION IN XC OF THE FIRST COLLOCATION POINT OF THE C I-TH SUBINTERVAL * C FOR EACH COLLOCATION POINT OF THE I-TH SUBINTERVAL, C CONSTRUCT THE CORRESPONDING SUBBLOCK ROW OF THE I-TH BLOCK. * DO 250 IP = 1,KRPT * IPPOS = IPOS + (IP-1)*NPDE * C IPPOS+1 IS THE POSITION IN PW CORRESONDING TO THE START OF THE C OF SUBBLOCKS ASSOCIATED WITH THE IP-TH COLLOCATION POINT OF TH C SUBINTERVAL * IK = (IC+IP)*KORD*3 * C IK IS THE POSITION IN A OF THE BASIS FUNCTIGN INFORMATION CORR C ING TO THE IP-TH COLLOCATION POINT OF THE I-TH SUBINTERVAL * ICC = IC + IP + 1 * C ICC IS THE ABSOLUTE POSITION IN XC OF THE IP-TH COLLKCATION PO C OF THE I-TH SUBINTERVAL * CALL EVAL(ICC,NPDE,C,UVAL,A,ILEFT) * C EVAL RETURNS AN EVALUATION OF THE SOLUTION COMPONENTS AND THEI C FIRST AND SECOND DERIVATIVES IN THE MATRIX UVAL, WITH THE EVAL C UATIONS BEING DONE AT THE ICC-TH COLLOCATION POINT. * IF ( MITER .EQ. 1 ) * CALL DERIVF(T,XC(ICC),UVAL,UVAL(1,2),UVAL(1,3), * DFDU,DFDUX,DFDUXX,NPDE) * C DERIVF RETURNS AN EVALUATION OF THE JACOBIAN OF F AS WELL C AS THE DERIVATIVES WITH RESPECT TO THE FIRST AND SECOND C DERIVATIVES OF THE SOLUTION IN THE MATRICES DFDU, DFDUX, C AND DFDUXX. THE EVALUATION IS DONE AT THE ICC-TH COLLOCATION C POINT. * IF ( MITER .EQ. 2 ) * CALL DIFFF(T,XC(ICC),ICC,UVAL,UVAL(1,2),UVAL(1,3), * DFDU,DFDUX,DFDUXX,NPDE,CMAX,SAVE2) * C IN THIS CASE THE JACOBIAN IS CALCULATED NUMERICALLY BY THE DIF C ROUTINE * DO 200 J = 1,KORD * C FOR J RANGING OVER THE KORD SUBBLOCKS OF THE IP-TH SUBBLOCK C OF THE I-TH SUBINTERVAL CONSTRUCT THE CORRESPONDING COMPONEN C PW. (THESE CORRESPOND TO THE KORD NON-ZERO BASIS FUNCTIONS O C I-TH SUBINTERVAL). * JK = IK + J * C JK IS THE POSITION IN A OF THE BASIS FUNCTION INFORMATION C CORRESPONDING TO THE J-TH NON-ZERO BASIS FUNCTION ASSOCIATED C WITH THE IP-TH COLLOCATION POINT OF THE I-TH SUBINTERVAL * JK2 = JK + KORD * C JK2 IS THE POSITION IN A OF THE FIRST DERIVATIVE INFORMATION C FOR THE ABOVE BASIS FUNCTION * JK3 = JK2 + KORD * C JK3 IS THE POSITION IN A OF THE SECOND DERIVATIVE INFORMATIO C FOR THE ABOVE BASIS FUNCTION * JPOS = IPPOS + (J-1)*NROWS*NPDE * C JPOS+1 IS THE POSITION IN PW OF THE START OF THE J-TH SUBBLO C THE ROW OF SUBBLOCKS CORREPONDING TO THE IP-TH COLLOCATION P C OF THE I-TH SUBINTERVAL * DO 150 JJ = 1,NPDE * DO 100 KK = 1,NPDE * JJPOS = JPOS + (JJ-1)*NROWS + KK * C JJPOS IS THE POSITION IN PW OF THE (JJ,KK)-TH ENTRY OF T C J-TH SUBBLOCK ABOVE. * C WE NOW ASSIGN THE APPROPRIATE COMPONENT OF (DG/DC), IN T C OF DF/DU, DF/DUX, AND DF/DUXX AND A(K,I,J) TO PW. * PW(JJPOS) = DFDU(KK,JJ)*A(JK) + DFDUX(KK,JJ)*A(JK2) * + DFDUXX(KK,JJ)*A(JK3) * 100 CONTINUE * 150 CONTINUE * 200 CONTINUE * 250 CONTINUE * 300 CONTINUE * C THE NEW VERSION DOES NOT ALLOW COLLOCATION AT THE ENDPOINTS SINCE C THE SIZES OF THE TOP AND BOTTOM BLOCKS WOULD MAKE THEM UNSUITABLE C FOR USE IN COLROW. THUS AT THIS POINT THERE IS NO NEED TO MAKE C ADJUSTMENTS TO THE ENTRIES OF PW BASED ON WHETHER OR NOT SOME C COMPONENTS OF THE SYSTEM OF PDE'S HAVE NO BOUNDARY CONDITIONS IMPO C ON THEM, AS IS DONE IS THE ORIGINAL VERSION. * C AT THIS POINT PW CONTAINS (DG/DC). THE REMAINING STEPS ARE TO MULT C BY THE CONSTANT, CON, AND THEN ADD ON THE MATRIX A(Y,T) * DO 85 I = 1,SIZE * PW(I) = PW(I)*CON * 85 CONTINUE * C ADD MATRIX A(C,T) TO PW. * CALL ADDA (PW, A, BC, NPDE) * C DO LU DECOMPOSITION ON PW. * NELTP1 = NELTOP + 1 BOTPS1 = BOTPOS + 1 IFLAG = 0 * CALL CRDCMP (NEQN, PW, NPDE, NOVRLP, PW(NELTP1), NROWS, NCLBLK, * NINT,PW(BOTPS1), NPDE, IPIV, IFLAG) * RETURN END SUBROUTINE STIFIB (N0,Y,YMAX,ERROR,SAVE1,SAVE2,SAVE3, * PW,IPIV,WORK,IWORK) C----------------------------------------------------------------------- C STIFIB PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM, C A(Y,T)*(DY/DT) = G(Y,T), WHERE Y = (Y(1),Y(2), ... ,Y(N)). C STIFIB IS FOR USE WHEN THE MATRICES A AND DG/DY HAVE AN ALMOST BLOCK C DIAGONAL FORM. THE DEPENDENCE OF A(Y,T) ON Y IS ASSUMED TO BE WEAK. C C--------------------------------------------------------------------- C THE MODIFICATION TO THE ORIGINAL STIFIB PROGRAM CONSISTS ENTIRELY OF T C REMOVAL OF THE BAND LINEAR EQUATION SOLVER, SOLB, AND ITS SUBSTITUTION C BY CRSLVE, THE ALMOST BLOCK DIAGONAL LINEAR EQUATION SOLVER. THE CHANG C IN THE STIFIB ROUTINE ARE LIMITED TO: C C (1) REPLACEMENT OF THE CALL TO SOLV BY A CALL TO CRSLVE C (2) THE REMOVAL OF SOME VARIABLES USED SOLELY TO PROVIDE C INFORMATION TO SOLB AND THE INSERTION OF NEW VARIABLES C TO PASS INFORMATION TO CRSLVE. C----------------------------------------------------------------------- C C REFERENCE C C HINDMARSH, A.C., PRELIMINARY DOCUMENTATION OF GEARIB.. SOLUTION C OF IMPLICIT SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS WITH C BANDED JACOBIANS, LAWRENCE LIVERMORE LAB, UCID-30130, FEBRUARY C 1976. C C COMMUNICATION WITH STIFIB IS DONE WITH THE FOLLOWING VARIABLES.. C C Y AN N0 BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES C AND THEIR SCALED DERIVATIVES. LMAX IS 13 FOR THE ADAMS C METHODS AND 6 FOR THE GEAR METHODS. LMAX - 1 = MAXDER C IS THE MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET. C Y(I,J+1) CONTAINS THE J-TH DERIVATIVE OF Y(I), SCALED BY C H**J/FACTORIAL(J) (J = 0,1,...,NQ). C N0 A CONSTANT INTEGER .GE. N, USED FOR DIMENSIONING PURPOSES. C T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN. C H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. C HMIN, THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEP SIZE C HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY C TIME, BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE. C EPS THE RELATIVE ERROR BOUND. SEE DESCRIPTION IN PDECOL. C UROUND THE UNIT ROUNDOFF OF THE MACHINE. C N THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS. C MF THE METHOD FLAG. SEE DESCRIPTION IN PDECOL. C KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. C 0 THE STEP WAS SUCCESFUL. C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED C WITH ABS(H) = HMIN. C -2 THE REQUESTED ERROR IS SMALLER THAN CAN C BE HANDLED FOR THIS PROBLEM. C -3 CORRECTOR CONVERGENCE COULD NOT BE C ACHIEVED FOR ABS(H) = HMIN. C -4 SINGULAR A-MATRIX ENCOUNTERED. C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND C THE Y ARRAY ARE AS OF THE BEGINNING OF THE LAST C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED. C JSTART AN INTEGER USED ON INPUT AND OUTPUT. C ON INPUT, IT HAS THE FOLLOWING VALUES AND MEANINGS.. C 0 PERFORM THE FIRST STEP. C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST. C .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF C H, EPS, N, AND/OR MF. C ON EXIT, JSTART IS NQ, THE CURRENT ORDER OF THE METHOD. C YMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL C ERRORS IN Y ARE COMPARED. C ERROR AN ARRAY OF N ELEMENTS. ERROR(I)/TQ(2) IS THE ESTIMATED C ONE-STEP ERROR IN Y(I). C SAVE1,SAVE2,SAVE3 THREE WORKING STORAGE ARRAYS, EACH OF LENGTH N. C PW A BLOCK OF LOCATIONS USED FOR THE CHORD ITERATION C MATRIX. SEE DESCRIPTION IN PDECOL. C IPIV AN INTEGER ARRAY OF LENGTH N FOR PIVOT INFORMATION. C WORK,IWORK WORKING ARRAYS WHICH ARE USED TO PASS APPROPRIATE C ARRAYS TO OTHER SUBROUTINES. SEE DESCRIPTION IN C PDECOL. C C PACKAGE ROUTINES CALLED.. COSET,DIFFUN,PSETIB,RES,CRSLVE C USER ROUTINES CALLED.. NONE C CALLED BY.. PDECOL C FORTRAN FUNCTIONS USED.. ABS,MAX,MIN,FLOAT C----------------------------------------------------------------------- DIMENSION Y(N0,1),YMAX(N0),ERROR(N0),SAVE1(N0),SAVE2(N0), * SAVE3(N0),PW(1),IPIV(1),WORK(1),IWORK(1) * COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11, * IW12,IW13,IW14,IW15,IW16,IW17,IW18 COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,N,MF,KFLAG,JSTART COMMON /GEAR0/ HUSED,NQUSED,NSTEP,NFE,NJE COMMON /OPTION/ NOGAUS,MAXDER COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT * INTEGER BOTPOS,SIZE,NELTP1,BOTPS1 DIMENSION EL(13),TQ(4) DATA EL(2)/1./, OLDL0/1./, TQ(1)/0./, IER/0/ KFLAG = 0 TOLD = T IF (JSTART .GT. 0) GO TO 200 IF (JSTART .NE. 0) GO TO 120 C----------------------------------------------------------------------- C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT IS C CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2 C FOR THE NEXT INCREASE. C----------------------------------------------------------------------- NQ = 1 IER = 0 CALL DIFFUN (N, T, Y, SAVE1, IER, PW, IPIV, WORK, IWORK) IF ( IER .NE. 0 ) GO TO 685 DO 110 I = 1,N 110 Y(I,2) = H*SAVE1(I) METH = MF/10 MITER = MF - 10*METH L = 2 IDOUB = 3 RMAX = 1.D+4 RC = 0. CRATE = 1. EPSOLD = EPS HOLD = H MFOLD = MF NOLD = N NSTEP = 0 NSTEPJ = 0 NFE = 0 NJE = 1 IRET = 3 GO TO 130 C----------------------------------------------------------------------- C IF THE CALLER HAS CHANGED METH, COSET IS CALLED TO SET C THE COEFFICIENTS OF THE METHOD. IF THE CALLER HAS CHANGED C N, EPS, OR METH, THE CONSTANTS E, EDN, EUP, AND BND MUST BE RESET. C E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS C TO TEST FOR INCREASING THE ORDER, EDN FOR DECREASING THE ORDER. C BND IS USED TO TEST FOR CONVERGENCE OF THE CORRECTOR ITERATES. C IF THE CALLER HAS CHANGED H, Y MUST BE RESCALED. C IF H OR METH HAS BEEN CHANGED, IDOUB IS RESET TO L + 1 TO PREVENT C FURTHER CHANGES IN H FOR THAT MANY STEPS. C----------------------------------------------------------------------- 120 IF (MF .EQ. MFOLD) GO TO 150 MEO = METH MIO = MITER METH = MF/10 MITER = MF - 10*METH MFOLD = MF IF (MITER .NE. MIO) IWEVAL = MITER IF (METH .EQ. MEO) GO TO 150 IDOUB = L + 1 IRET = 1 130 CALL COSET (METH, NQ, EL, TQ) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDL0 OLDL0 = EL(1) 140 FN = FLOAT(N) EDN = FN*(TQ(1)*EPS)**2 E = FN*(TQ(2)*EPS)**2 EUP = FN*(TQ(3)*EPS)**2 BND = FN*(TQ(4)*EPS)**2 GO TO (160, 170, 200), IRET 150 IF ((EPS .EQ. EPSOLD) .AND. (N .EQ. NOLD)) GO TO 160 EPSOLD = EPS NOLD = N IRET = 1 GO TO 140 160 IF (H .EQ. HOLD) GO TO 200 RH = H/HOLD H = HOLD IREDO = 3 GO TO 175 170 RH = MAX(RH,HMIN/ ABS(H)) 175 RH = MIN(RH,HMAX/ ABS(H),RMAX) R1 = 1. DO 180 J = 2,L R1 = R1*RH DO 180 I = 1,N 180 Y(I,J) = Y(I,J)*R1 H = H*RH RC = RC*RH IDOUB = L + 1 IF (IREDO .EQ. 0) GO TO 690 C----------------------------------------------------------------------- C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY C MULTIPLYING THE Y ARRAY BY THE PASCAL TRIANGLE MATRIX. C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1). C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, OR THE CALLER HAS C CHANGED MITER, IWEVAL IS SET TO MITER TO FORCE PW TO BE UPDATED. C IN ANY CASE, PW IS UPDATED AT LEAST EVERY 40-TH STEP. C PW IS THE CHORD ITERATION MATRIX A - H*EL(1)*(DG/DY). C----------------------------------------------------------------------- 200 IF ( ABS(RC-1.) .GT. 0.3) IWEVAL = MITER IF (NSTEP .GE. NSTEPJ+40) IWEVAL = MITER T = T + H DO 210 J1 = 1,NQ DO 210 J2 = J1,NQ J = (NQ + J1) - J2 DO 210 I = 1,N 210 Y(I,J) = Y(I,J) + Y(I,J+1) C----------------------------------------------------------------------- C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS C MADE ON THE R.M.S. NORM OF EACH CORRECTION, USING BND, WHICH C IS DEPENDENT ON EPS. THE SUM OF THE CORRECTIONS IS ACCUMULATED C IN THE VECTOR ERROR(I). THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR C LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1. C THE UPDATED H*YDOT IS STORED IN SAVE2. C----------------------------------------------------------------------- 220 DO 230 I = 1,N SAVE2(I) = Y(I,2) 230 ERROR(I) = 0. M = 0 CALL RES (T, H, Y, SAVE2, SAVE3, NPDE, NCPTS, WORK(IW1), IWORK, * WORK, WORK(IW14), WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9)) NFE = NFE + 1 IF (IWEVAL .LE. 0) GO TO 350 C----------------------------------------------------------------------- C IF INDICATED, THE MATRIX PW IS REEVALUATED BEFORE STARTING THE C CORRECTOR ITERATION. IWEVAL IS SET TO 0 AS AN INDICATOR C THAT THIS HAS BEEN DONE. PW IS COMPUTED AND PROCESSED IN PSETIB. C----------------------------------------------------------------------- IWEVAL = 0 RC = 1. NJE = NJE + 1 NSTEPJ = NSTEP CON = -H*EL(1) CALL PSETIB (Y, PW, N0, CON, MITER, IER, WORK(IW1), IWORK, * WORK(IW3),WORK(IW9),SAVE2,IPIV,YMAX,WORK(IW11),WORK(IW12), * WORK(IW13),WORK(IW16),WORK(IW14),WORK(IW15),WORK,NPDE) IF (IER .NE. 0) GO TO 420 C----------------------------------------------------------------------- C COMPUTE THE CORRECTOR ERROR, R SUB M, AND SOLVE THE LINEAR SYSTEM C WITH THAT AS RIGHT-HAND SIDE AND PW AS COEFFICIENT MATRIX, C USING THE LU DECOMPOSITION OF PW. C----------------------------------------------------------------------- 350 CONTINUE NELTP1 = NELTOP + 1 BOTPS1 = BOTPOS + 1 * CALL CRSLVE( PW, NPDE, NOVRLP, PW(NELTP1), NROWS, NCLBLK, * NINT, PW(BOTPS1), NPDE, IPIV, SAVE3, SAVE3) * 370 D = 0. DO 380 I = 1,N ERROR(I) = ERROR(I) + SAVE3(I) D = D + (SAVE3(I)/YMAX(I))**2 SAVE1(I) = Y(I,1) + EL(1)*ERROR(I) 380 SAVE2(I) = Y(I,2) + ERROR(I) C----------------------------------------------------------------------- C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. C----------------------------------------------------------------------- 400 IF (M .NE. 0) CRATE = MAX(.9*CRATE,D/D1) IF ((D*MIN(1.E-0,2.*CRATE)) .LE. BND) GO TO 450 D1 = D M = M + 1 IF (M .EQ. 3) GO TO 410 CALL RES(T, H, SAVE1, SAVE2, SAVE3, NPDE, NCPTS, WORK(IW1), IWORK, * WORK, WORK(IW14), WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9)) GO TO 350 C----------------------------------------------------------------------- C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. C IF THE MATRIX PW IS NOT UP TO DATE, IT IS REEVALUATED FOR THE C NEXT TRY. OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF NOT, A C NO-CONVERGENCE EXIT IS TAKEN. C----------------------------------------------------------------------- 410 NFE = NFE + 2 IF (IWEVAL .EQ. -1) GO TO 440 420 T = TOLD RMAX = 2. DO 430 J1 = 1,NQ DO 430 J2 = J1,NQ J = (NQ + J1) - J2 DO 430 I = 1,N 430 Y(I,J) = Y(I,J) - Y(I,J+1) IF ( ABS(H) .LE. HMIN*1.00001) GO TO 680 RH = .25 IREDO = 1 GO TO 170 440 IWEVAL = MITER GO TO 220 C----------------------------------------------------------------------- C THE CORRECTOR HAS CONVERGED. IWEVAL IS SET TO -1 TO SIGNAL C THAT PW MAY NEED UPDATING ON SUBSEQUENT STEPS. THE ERROR TEST C IS MADE AND CONTROL PASSES TO STATEMENT 500 IF IT FAILS. C----------------------------------------------------------------------- 450 IWEVAL = -1 NFE = NFE + M D = 0. DO 460 I = 1,N 460 D = D + (ERROR(I)/YMAX(I))**2 IF (D .GT. E) GO TO 500 C----------------------------------------------------------------------- C AFTER A SUCCESSFUL STEP, UPDATE THE Y ARRAY. C CONSIDER CHANGING H IF IDOUB = 1. OTHERWISE DECREASE IDOUB BY 1. C IF IDOUB IS THEN 1 AND NQ .LT. MAXDER, THEN ERROR IS SAVED FOR C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A C FACTOR OF AT LEAST 1.1. IF NOT, IDOUB IS SET TO 10 TO PREVENT C TESTING FOR THAT MANY STEPS. C----------------------------------------------------------------------- KFLAG = 0 IREDO = 0 NSTEP = NSTEP + 1 HUSED = H NQUSED = NQ DO 470 J = 1,L DO 470 I = 1,N 470 Y(I,J) = Y(I,J) + EL(J)*ERROR(I) IF (IDOUB .EQ. 1) GO TO 520 IDOUB = IDOUB - 1 IF (IDOUB .GT. 1) GO TO 700 IF (L .EQ. LMAX) GO TO 700 DO 490 I = 1,N 490 Y(I,LMAX) = ERROR(I) GO TO 700 C----------------------------------------------------------------------- C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES. C RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR C ONE LOWER ORDER. C----------------------------------------------------------------------- 500 KFLAG = KFLAG - 1 T = TOLD DO 510 J1 = 1,NQ DO 510 J2 = J1,NQ J = (NQ + J1) - J2 DO 510 I = 1,N 510 Y(I,J) = Y(I,J) - Y(I,J+1) RMAX = 2. IF ( ABS(H) .LE. HMIN*1.00001) GO TO 660 IF (KFLAG .LE. -3) GO TO 640 IREDO = 2 PR3 = 1.E+20 GO TO 540 C----------------------------------------------------------------------- C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS C PR1, PR2, AND PR3 ARE COMPUTED, BY WHICH H COULD BE DIVIDED C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. C IN THE CASE OF FAILURE, PR3 = 1.E20 TO AVOID AN ORDER INCREASE. C THE SMALLEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE C ADDITIONAL SCALED DERIVATIVE. C----------------------------------------------------------------------- 520 PR3 = 1.E+20 IF (L .EQ. LMAX) GO TO 540 D1 = 0. DO 530 I = 1,N 530 D1 = D1 + ((ERROR(I) - Y(I,LMAX))/YMAX(I))**2 ENQ3 = .5/ FLOAT(L+1) PR3 = ((D1/EUP)**ENQ3)*1.4 + 1.4E-6 540 ENQ2 = .5/ FLOAT(L) PR2 = ((D/E)**ENQ2)*1.2 + 1.2E-6 PR1 = 1.E+20 IF (NQ .EQ. 1) GO TO 560 D = 0. DO 550 I = 1,N 550 D = D + (Y(I,L)/YMAX(I))**2 ENQ1 = .5/ FLOAT(NQ) PR1 = ((D/EDN)**ENQ1)*1.3 + 1.3E-6 560 IF (PR2 .LE. PR3) GO TO 570 IF (PR3 .LT. PR1) GO TO 590 GO TO 580 570 IF (PR2 .GT. PR1) GO TO 580 NEWQ = NQ RH = 1./PR2 GO TO 620 580 NEWQ = NQ - 1 RH = 1./PR1 GO TO 620 590 NEWQ = L RH = 1./PR3 IF (RH .LT. 1.1) GO TO 610 DO 600 I = 1,N 600 Y(I,NEWQ+1) = ERROR(I)*EL(L)/ FLOAT(L) GO TO 630 610 IDOUB = 10 GO TO 700 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1)) GO TO 610 C----------------------------------------------------------------------- C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. C IN ANY CASE H IS RESET ACCORDING TO RH AND THE Y ARRAY IS RESCALED. C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE. C----------------------------------------------------------------------- IF (NEWQ .EQ. NQ) GO TO 170 630 NQ = NEWQ L = NQ + 1 IRET = 2 GO TO 130 C----------------------------------------------------------------------- C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED. C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE C Y ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED. C AFTER A TOTAL OF 7 FAILURES, AN EXIT IS TAKEN WITH KFLAG = -2. C----------------------------------------------------------------------- 640 IF (KFLAG .EQ. -7) GO TO 670 RH = .1 RH = MAX(HMIN/ ABS(H),RH) H = H*RH IER = 0 CALL DIFFUN (N, T, Y, SAVE1, IER, PW, IPIV, WORK, IWORK) IF ( IER .NE. 0 ) GO TO 685 NJE = NJE + 1 DO 650 I = 1,N 650 Y(I,2) = H*SAVE1(I) IWEVAL = MITER IDOUB = 10 IF (NQ .EQ. 1) GO TO 200 NQ = 1 L = 2 IRET = 3 GO TO 130 C----------------------------------------------------------------------- C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. C----------------------------------------------------------------------- 660 KFLAG = -1 GO TO 700 670 KFLAG = -2 GO TO 700 680 KFLAG = -3 GO TO 700 685 KFLAG = -4 GO TO 700 690 RMAX = 10. 700 HOLD = H JSTART = NQ RETURN END BLOCK DATA C----------------------------------------------------------------------- C IN THE FOLLOWING DATA STATEMENT, SET.. C LOUT = THE LOGICAL UNIT NUMBER FOR THE OUTPUT OF MESSAGES DURING C THE INTEGRATION. C NOGAUS = SET .EQ. 1 IF THE GAUSS-LEGENDRE COLLOCATION POINTS ARE C NOT DESIRED WHEN NCC = 2 (SEE PDECOL AND COLPNT). C MAXDER = SET .EQ. 5. ITS VALUE REPRESENTS THE MAXIMUM ORDER OF C THE TIME INTEGRATION ALLOWED. ITS VALUE AFFECTS THE STOR- C AGE REQUIRED IN /HISTRY/ AND MAY BE CHANGED IF DESIRED C (SEE COSET FOR RESTRICTIONS). C UROUND = THE UNIT ROUNDOFF OF THE MACHINE, I.E. THE SMALLEST C POSITIVE U SUCH THAT 1. + U .NE. 1. ON THE MACHINE. C----------------------------------------------------------------------- COMMON /GEAR1/ DUM(5),UROUND,IDUM(4) COMMON /OPTION/ NOGAUS,MAXDER COMMON /IOUNIT/ LOUT DATA LOUT,NOGAUS,MAXDER,UROUND/6,0,5,2.22E-16/ END SUBROUTINE CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK, * NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,IFLAG) C C*************************************************************** C C C R D C M P DECOMPOSES THE ALMOST BLOCK DIAGONAL MATRIX A C USING MODIFIED ALTERNATE ROW AND COLUMN ELIMINATION WITH C PARTIAL PIVOTING. THE MATRIX A IS STORED IN THE ARRAYS C TOPBLK, ARRAY, AND BOTBLK. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM, C GIVEN BY NBLOKS*NRWBLK + NOVRLP C C TOPBLK - REAL(NRWTOP,NOVRLP) C THE FIRST BLOCK OF THE ALMOST BLOCK C DIAGONAL MATRIX A TO BE DECOMPOSED C C NRWTOP - INTEGER C NUMBER OF ROWS IN THE BLOCK TOPBLK C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN WHICH SUCC- C ESSIVE BLOCKS OVERLAP, WHERE C NOVRLP = NRWTOP + NRWBOT C C ARRAY - REAL(NRWBLK,NCLBLK,NBLOKS) C ARRAY(,,K) CONTAINS THE K-TH NRWBLK C BY NCLBLK BLOCK OF THE MATRIX A C C NRWBLK - INTEGER C NUMBER OF ROWS IN K-TH BLOCK C C NCLBLK - INTEGER C NUMBER OF COLUMNS IN K-TH BLOCK C C NBLOKS - INTEGER C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN C THE MATRIX A C C BOTBLK - REAL(NRWBOT,NOVRLP) C THE LAST BLOCK OF THE MATRIX A C C NRWBOT - INTEGER C NUMBER OF ROWS IN THE BLOCK BOTBLK C C PIVOT - INTEGER(N) C WORK SPACE C C *** ON RETURN ... C C TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE C DESIRED DECOMPOSITION OF THE MATRIX A C (IF IFLAG = 0) C C PIVOT - INTEGER(N) C RECORDS THE PIVOTING INDICES DETER- C MINED IN THE DECOMPOSITION C C IFLAG - INTEGER C = 1, IF INPUT PARAMETERS ARE INVALID C = -1, IF MATRIX IS SINGULAR C = 0, OTHERWISE C C*************************************************************** C REAL TOPBLK,ARRAY,BOTBLK REAL ROWMAX,ROWPIV,ROWMLT,COLMAX,COLPIV REAL SWAP,COLMLT,PIVMAX,ZERO,TEMPIV INTEGER PIVOT(1) DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1), * BOTBLK(NRWBOT,1) DATA ZERO/0.0D0/ C C*************************************************************** C C **** DEFINE THE CONSTANTS USED THROUGHOUT **** C C*************************************************************** C IFLAG = 0 PIVMAX = ZERO NRWTP1 = NRWTOP+1 NROWEL = NRWBLK-NRWTOP NRWEL1 = NROWEL+1 NVRLP0 = NOVRLP-1 C C*************************************************************** C C **** CHECK VALIDITY OF THE INPUT PARAMETERS.... C C IF PARAMETERS ARE INVALID THEN TERMINATE AT 10; C ELSE CONTINUE AT 100. C C*************************************************************** C IF(N.NE.NBLOKS*NRWBLK+NOVRLP)GO TO 10 IF(NOVRLP.NE.NRWTOP+NRWBOT)GO TO 10 IF(NCLBLK.NE.NOVRLP+NRWBLK)GO TO 10 IF(NOVRLP.GT.NRWBLK)GO TO 10 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PARAMETERS ARE ACCEPTABLE - CONTINUE AT 100. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C GO TO 100 10 CONTINUE C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PARAMETERS ARE INVALID. SET IFLAG = 1, AND TERMINATE C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IFLAG = 1 RETURN 100 CONTINUE C C*************************************************************** C C **** FIRST, IN TOPBLK.... C C*************************************************************** C C *** APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN C PIVOTING .... C C*************************************************************** C DO 190 I = 1,NRWTOP IPLUS1 = I+1 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DETERMINE COLUMN PIVOT AND PIVOT INDEX C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IPVT = I COLMAX = ABS(TOPBLK(I,I)) DO 110 J = IPLUS1,NOVRLP TEMPIV = ABS(TOPBLK(I,J)) IF(TEMPIV.LE.COLMAX)GO TO 110 IPVT = J COLMAX = TEMPIV 110 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TEST FOR SINGULARITY: C C IF SINGULAR THEN TERMINATE AT 1000; C ELSE CONTINUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF(PIVMAX+COLMAX.EQ.PIVMAX)GO TO 1000 PIVMAX = MAX(COLMAX,PIVMAX) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF NECESSARY INTERCHANGE COLUMNS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PIVOT(I) = IPVT IF(IPVT.EQ.I)GO TO 140 DO 120 L = I,NRWTOP SWAP = TOPBLK(L,IPVT) TOPBLK(L,IPVT) = TOPBLK(L,I) TOPBLK(L,I) = SWAP 120 CONTINUE DO 130 L = 1,NRWBLK SWAP = ARRAY(L,IPVT,1) ARRAY(L,IPVT,1) = ARRAY(L,I,1) ARRAY(L,I,1) = SWAP 130 CONTINUE 140 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COMPUTE MULTIPLIERS AND PERFORM COLUMN C ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COLPIV = TOPBLK(I,I) DO 180 J = IPLUS1,NOVRLP COLMLT = TOPBLK(I,J)/COLPIV TOPBLK(I,J) = COLMLT IF(IPLUS1.GT.NRWTOP)GO TO 160 DO 150 L = IPLUS1,NRWTOP TOPBLK(L,J) = TOPBLK(L,J)-COLMLT*TOPBLK(L,I) 150 CONTINUE 160 CONTINUE DO 170 L = 1,NRWBLK ARRAY(L,J,1) = ARRAY(L,J,1)-COLMLT*ARRAY(L,I,1) 170 CONTINUE 180 CONTINUE 190 CONTINUE C C*************************************************************** C C **** IN EACH BLOCK ARRAY(,,K).... C C*************************************************************** C INCR = 0 DO 395 K = 1,NBLOKS KPLUS1 = K+1 C C ***************************************************** C C *** FIRST APPLY NRWBLK-NRWTOP ROW ELIMINATIONS WITH C ROW PIVOTING.... C C ***************************************************** C DO 270 J = NRWTP1,NRWBLK JPLUS1 = J+1 JMINN = J-NRWTOP C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DETERMINE ROW PIVOT AND PIVOT INDEX C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IPVT = JMINN ROWMAX = ABS(ARRAY(JMINN,J,K)) LOOP = JMINN+1 DO 210 I = LOOP,NRWBLK TEMPIV = ABS(ARRAY(I,J,K)) IF(TEMPIV.LE.ROWMAX)GO TO 210 IPVT = I ROWMAX = TEMPIV 210 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TEST FOR SINGULARITY: C C IF SINGULAR THEN TERMINATE AT 1000; C ELSE CONTINUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF(PIVMAX+ROWMAX.EQ.PIVMAX)GO TO 1000 PIVMAX = MAX(ROWMAX,PIVMAX) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF NECESSARY INTERCHANGE ROWS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCRJ = INCR+J PIVOT(INCRJ) = INCR+IPVT+NRWTOP IF(IPVT.EQ.JMINN)GO TO 230 DO 220 L = J,NCLBLK SWAP = ARRAY(IPVT,L,K) ARRAY(IPVT,L,K) = ARRAY(JMINN,L,K) ARRAY(JMINN,L,K) = SWAP 220 CONTINUE 230 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COMPUTE MULTIPLERS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ROWPIV = ARRAY(JMINN,J,K) DO 240 I = LOOP,NRWBLK ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV 240 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PERFORM ROW ELIMINATION WITH COLUMN INDEXING C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 260 L = JPLUS1,NCLBLK ROWMLT = ARRAY(JMINN,L,K) DO 250 I = LOOP,NRWBLK ARRAY(I,L,K) = ARRAY(I,L,K) * -ROWMLT*ARRAY(I,J,K) 250 CONTINUE 260 CONTINUE 270 CONTINUE C C ***************************************************** C C *** NOW APPLY NRWTOP COLUMN ELIMINATIONS WITH C COLUMN PIVOTING.... C C ***************************************************** C DO 390 I = NRWEL1,NRWBLK IPLUSN = I+NRWTOP IPLUS1 = I+1 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DETERMINE COLUMN PIVOT AND PIVOT INDEX C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IPVT = IPLUSN COLMAX = ABS(ARRAY(I,IPVT,K)) LOOP = IPLUSN+1 DO 310 J = LOOP,NCLBLK TEMPIV = ABS(ARRAY(I,J,K)) IF(TEMPIV.LE.COLMAX)GO TO 310 IPVT = J COLMAX = TEMPIV 310 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TEST FOR SINGULARITY: C C IF SINGULAR THEN TERMINATE AT 1000; C ELSE CONTINUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF(PIVMAX+COLMAX.EQ.PIVMAX)GO TO 1000 PIVMAX = MAX(COLMAX,PIVMAX) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF NECESSARY INTERCHANGE COLUMNS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCRN = INCR+IPLUSN PIVOT(INCRN) = INCR+IPVT IRWBLK = IPLUSN-NRWBLK IF(IPVT.EQ.IPLUSN)GO TO 340 DO 315 L = I,NRWBLK SWAP = ARRAY(L,IPVT,K) ARRAY(L,IPVT,K) = ARRAY(L,IPLUSN,K) ARRAY(L,IPLUSN,K) = SWAP 315 CONTINUE IPVBLK = IPVT-NRWBLK IF(K.EQ.NBLOKS)GO TO 330 DO 320 L = 1,NRWBLK SWAP = ARRAY(L,IPVBLK,KPLUS1) ARRAY(L,IPVBLK,KPLUS1) * = ARRAY(L,IRWBLK,KPLUS1) ARRAY(L,IRWBLK,KPLUS1) = SWAP 320 CONTINUE GO TO 340 330 CONTINUE DO 335 L = 1,NRWBOT SWAP = BOTBLK(L,IPVBLK) BOTBLK(L,IPVBLK) = BOTBLK(L,IRWBLK) BOTBLK(L,IRWBLK) = SWAP 335 CONTINUE 340 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COMPUTE MULTIPLIERS AND PERFORM COLUMN C ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COLPIV = ARRAY(I,IPLUSN,K) DO 380 J = LOOP,NCLBLK COLMLT = ARRAY(I,J,K)/COLPIV ARRAY(I,J,K) = COLMLT IF(I.EQ.NRWBLK)GO TO 350 DO 345 L = IPLUS1,NRWBLK ARRAY(L,J,K) = ARRAY(L,J,K) * -COLMLT*ARRAY(L,IPLUSN,K) 345 CONTINUE 350 CONTINUE JRWBLK = J-NRWBLK IF(K.EQ.NBLOKS)GO TO 370 DO 360 L = 1,NRWBLK ARRAY(L,JRWBLK,KPLUS1) = * ARRAY(L,JRWBLK,KPLUS1) * -COLMLT*ARRAY(L,IRWBLK,KPLUS1) 360 CONTINUE GO TO 380 370 CONTINUE DO 375 L = 1,NRWBOT BOTBLK(L,JRWBLK) = BOTBLK(L,JRWBLK) * -COLMLT*BOTBLK(L,IRWBLK) 375 CONTINUE 380 CONTINUE 390 CONTINUE INCR = INCR + NRWBLK 395 CONTINUE C C*************************************************************** C C **** FINALLY, IN BOTBLK.... C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C *** APPLY NRWBOT ROW ELIMINATIONS WITH ROW C PIVOTING.... C C IF BOT HAS JUST ONE ROW GO TO 500 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF(NRWBOT.EQ.1)GO TO 500 DO 470 J = NRWTP1,NVRLP0 JPLUS1 = J+1 JMINN = J-NRWTOP C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DETERMINE ROW PIVOT AND PIVOT INDEX C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IPVT = JMINN ROWMAX = ABS(BOTBLK(JMINN,J)) LOOP = JMINN+1 DO 410 I = LOOP,NRWBOT TEMPIV = ABS(BOTBLK(I,J)) IF(TEMPIV.LE.ROWMAX) GO TO 410 IPVT = I ROWMAX = TEMPIV 410 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TEST FOR SINGULARITY: C C IF SINGULAR THEN TERMINATE AT 1000; C ELSE CONTINUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF(PIVMAX+ROWMAX.EQ.PIVMAX)GO TO 1000 PIVMAX = MAX(ROWMAX,PIVMAX) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF NECESSARY INTERCHANGE ROWS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCRJ = INCR+J PIVOT(INCRJ) = INCR+IPVT+NRWTOP IF(IPVT.EQ.JMINN)GO TO 430 DO 420 L = J,NOVRLP SWAP = BOTBLK(IPVT,L) BOTBLK(IPVT,L) = BOTBLK(JMINN,L) BOTBLK(JMINN,L) = SWAP 420 CONTINUE 430 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COMPUTE MULTIPLIERS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ROWPIV = BOTBLK(JMINN,J) DO 440 I = LOOP,NRWBOT BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV 440 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PERFORM ROW ELIMINATION WITH COLUMN INDEXING C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 460 L = JPLUS1,NOVRLP ROWMLT = BOTBLK(JMINN,L) DO 450 I = LOOP,NRWBOT BOTBLK(I,L) = BOTBLK(I,L)-ROWMLT*BOTBLK(I,J) 450 CONTINUE 460 CONTINUE 470 CONTINUE 500 CONTINUE C C*************************************************************** C C DONE PROVIDED THE LAST ELEMENT IS NOT ZERO C C*************************************************************** C IF(PIVMAX+ABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN C C*************************************************************** C C **** MATRIX IS SINGULAR - SET IFLAG = - 1. C TERMINATE AT 1000. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 1000 CONTINUE IFLAG = -1 RETURN END SUBROUTINE CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK, * NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,B,X) C C*************************************************************** C C C R S L V E SOLVES THE LINEAR SYSTEM C A*X = B C USING THE DECOMPOSITION ALREADY GENERATED IN C R D C M P. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C TOPBLK - REAL(NRWTOP,NOVRLP) C OUTPUT FROM C R D C M P C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN WHICH SUCC- C ESSIVE BLOCKS OVERLAP, WHERE C NOVRLP = NRWTOP + NRWBOT C C NRWTOP - INTEGER C NUMBER OF ROWS IN THE BLOCK TOPBLK C C ARRAY - REAL(NRWBLK,NCLBLK,NBLOKS) C OUTPUT FROM C R D C M P C C NRWBLK - INTEGER C NUMBER OF ROWS IN K-TH BLOCK C C NCLBLK - INTEGER C NUMBER OF COLUMNS IN K-TH BLOCK C C NBLOKS - INTEGER C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN C THE MATRIX A C C BOTBLK - REAL(NRWBOT,NOVRLP) C OUTPUT FROM C R D C M P C C NRWBOT - INTEGER C NUMBER OF ROWS IN THE BLOCK BOTBLK C C PIVOT - INTEGER(N) C THE PIVOT VECTOR FROM C R D C M P C C B - REAL(N) C THE RIGHT HAND SIDE VECTOR C C X - REAL(N) C WORK SPACE C C *** ON RETURN ... C C X - REAL(N) C THE SOLUTION VECTOR C C*************************************************************** C REAL TOPBLK,ARRAY,BOTBLK,X,B REAL DOTPRD,XJ,XINCRJ,BINCRJ,SWAP INTEGER PIVOT(1) DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1), * BOTBLK(NRWBOT,1),B(1),X(1) C C*************************************************************** C C **** DEFINE THE CONSTANTS USED THROUGHOUT **** C C*************************************************************** C * NRWTP1 = NRWTOP+1 NRWBK1 = NRWBLK+1 NVRLP1 = NOVRLP+1 NRWTP0 = NRWTOP-1 NRWBT1 = NRWBOT+1 NROWEL = NRWBLK-NRWTOP NRWEL1 = NROWEL+1 NVRLP0 = NOVRLP-1 NBLKS1 = NBLOKS+1 NBKTOP = NRWBLK+NRWTOP C C*************************************************************** C C **** FORWARD RECURSION **** C C*************************************************************** C C *** FIRST, IN TOPBLK.... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD SOLUTION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 130 J = 1,NRWTOP X(J) = B(J)/TOPBLK(J,J) IF(J.EQ.NRWTOP)GO TO 120 XJ = -X(J) LOOP = J+1 DO 110 I = LOOP,NRWTOP B(I) = B(I)+TOPBLK(I,J)*XJ 110 CONTINUE 120 CONTINUE 130 CONTINUE C C ******************************************************** C C *** IN EACH BLOCK ARRAY(,,K).... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCR = 0 DO 280 K = 1,NBLOKS INCRTP = INCR+NRWTOP C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD MODIFICATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 220 J = 1,NRWTOP INCRJ = INCR+J XINCRJ = -X(INCRJ) DO 210 I = 1,NRWBLK B(INCRTP+I) = B(INCRTP+I)+ARRAY(I,J,K)*XINCRJ 210 CONTINUE 220 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 240 J = NRWTP1,NRWBLK INCRJ = INCR+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.EQ.INCRJ)GO TO 225 SWAP = B(INCRJ) B(INCRJ) = B(JPIVOT) B(JPIVOT) = SWAP 225 CONTINUE BINCRJ = -B(INCRJ) LOOP = J-NRWTP0 DO 230 I = LOOP,NRWBLK B(INCRTP+I) = B(INCRTP+I)+ARRAY(I,J,K)*BINCRJ 230 CONTINUE 240 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C