*DECK RUMACH REAL FUNCTION RUMACH () C***BEGIN PROLOGUE RUMACH C***PURPOSE Compute the unit roundoff of the machine. C***CATEGORY R1 C***TYPE SINGLE PRECISION (RUMACH-S, DUMACH-D) C***KEYWORDS MACHINE CONSTANTS C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C *Usage: C REAL A, RUMACH C A = RUMACH() C C *Function Return Values: C A : the unit roundoff of the machine. C C *Description: C The unit roundoff is defined as the smallest positive machine C number u such that 1.0 + u .ne. 1.0. This is computed by RUMACH C in a machine-independent manner. C C***REFERENCES (NONE) C***ROUTINES CALLED RUMSUM C***REVISION HISTORY (YYYYMMDD) C 19930216 DATE WRITTEN C 19930818 Added SLATEC-format prologue. (FNF) C 20030707 Added RUMSUM to force normal storage of COMP. (ACH) C***END PROLOGUE RUMACH C REAL U, COMP C***FIRST EXECUTABLE STATEMENT RUMACH U = 1.0E0 10 U = U*0.5E0 CALL RUMSUM(1.0E0, U, COMP) IF (COMP .NE. 1.0E0) GO TO 10 RUMACH = U*2.0E0 RETURN C----------------------- End of Function RUMACH ------------------------ END SUBROUTINE RUMSUM(A,B,C) C Routine to force normal storing of A + B, for RUMACH. REAL A, B, C C = A + B RETURN END *DECK SCFODE SUBROUTINE SCFODE (METH, ELCO, TESCO) C***BEGIN PROLOGUE SCFODE C***SUBSIDIARY C***PURPOSE Set ODE integrator coefficients. C***TYPE SINGLE PRECISION (SCFODE-S, DCFODE-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C SCFODE is called by the integrator routine to set coefficients C needed there. The coefficients for the current method, as C given by the value of METH, are set for all orders and saved. C The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2. C (A smaller value of the maximum order is also allowed.) C SCFODE is called once at the beginning of the problem, C and is not called again unless and until METH is changed. C C The ELCO array contains the basic method coefficients. C The coefficients el(i), 1 .le. i .le. nq+1, for the method of C order nq are stored in ELCO(i,nq). They are given by a genetrating C polynomial, i.e., C l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq. C For the implicit Adams methods, l(x) is given by C dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0. C For the BDF methods, l(x) is given by C l(x) = (x+1)*(x+2)* ... *(x+nq)/K, C where K = factorial(nq)*(1 + 1/2 + ... + 1/nq). C C The TESCO array contains test constants used for the C local error test and the selection of step size and/or order. C At order nq, TESCO(k,nq) is used for the selection of step C size at order nq - 1 if k = 1, at order nq if k = 2, and at order C nq + 1 if k = 3. C C***SEE ALSO SLSODE C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890501 Modified prologue to SLATEC/LDOC format. (FNF) C 890503 Minor cosmetic changes. (FNF) C 930809 Renamed to allow single/double precision versions. (ACH) C***END PROLOGUE SCFODE C**End INTEGER METH INTEGER I, IB, NQ, NQM1, NQP1 REAL ELCO, TESCO REAL AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ, 1 RQFAC, RQ1FAC, TSIGN, XPIN DIMENSION ELCO(13,12), TESCO(3,12) DIMENSION PC(12) C C***FIRST EXECUTABLE STATEMENT SCFODE GO TO (100, 200), METH C 100 ELCO(1,1) = 1.0E0 ELCO(2,1) = 1.0E0 TESCO(1,1) = 0.0E0 TESCO(2,1) = 2.0E0 TESCO(1,2) = 1.0E0 TESCO(3,12) = 0.0E0 PC(1) = 1.0E0 RQFAC = 1.0E0 DO 140 NQ = 2,12 C----------------------------------------------------------------------- C The PC array will contain the coefficients of the polynomial C p(x) = (x+1)*(x+2)*...*(x+nq-1). C Initially, p(x) = 1. C----------------------------------------------------------------------- RQ1FAC = RQFAC RQFAC = RQFAC/NQ NQM1 = NQ - 1 FNQM1 = NQM1 NQP1 = NQ + 1 C Form coefficients of p(x)*(x+nq-1). ---------------------------------- PC(NQ) = 0.0E0 DO 110 IB = 1,NQM1 I = NQP1 - IB 110 PC(I) = PC(I-1) + FNQM1*PC(I) PC(1) = FNQM1*PC(1) C Compute integral, -1 to 0, of p(x) and x*p(x). ----------------------- PINT = PC(1) XPIN = PC(1)/2.0E0 TSIGN = 1.0E0 DO 120 I = 2,NQ TSIGN = -TSIGN PINT = PINT + TSIGN*PC(I)/I 120 XPIN = XPIN + TSIGN*PC(I)/(I+1) C Store coefficients in ELCO and TESCO. -------------------------------- ELCO(1,NQ) = PINT*RQ1FAC ELCO(2,NQ) = 1.0E0 DO 130 I = 2,NQ 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/I AGAMQ = RQFAC*XPIN RAGQ = 1.0E0/AGAMQ TESCO(2,NQ) = RAGQ IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1 TESCO(3,NQM1) = RAGQ 140 CONTINUE RETURN C 200 PC(1) = 1.0E0 RQ1FAC = 1.0E0 DO 230 NQ = 1,5 C----------------------------------------------------------------------- C The PC array will contain the coefficients of the polynomial C p(x) = (x+1)*(x+2)*...*(x+nq). C Initially, p(x) = 1. C----------------------------------------------------------------------- FNQ = NQ NQP1 = NQ + 1 C Form coefficients of p(x)*(x+nq). ------------------------------------ PC(NQP1) = 0.0E0 DO 210 IB = 1,NQ I = NQ + 2 - IB 210 PC(I) = PC(I-1) + FNQ*PC(I) PC(1) = FNQ*PC(1) C Store coefficients in ELCO and TESCO. -------------------------------- DO 220 I = 1,NQP1 220 ELCO(I,NQ) = PC(I)/PC(2) ELCO(2,NQ) = 1.0E0 TESCO(1,NQ) = RQ1FAC TESCO(2,NQ) = NQP1/ELCO(1,NQ) TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ) RQ1FAC = RQ1FAC/FNQ 230 CONTINUE RETURN C----------------------- END OF SUBROUTINE SCFODE ---------------------- END *DECK SINTDY SUBROUTINE SINTDY (T, K, YH, NYH, DKY, IFLAG) C***BEGIN PROLOGUE SINTDY C***SUBSIDIARY C***PURPOSE Interpolate solution derivatives. C***TYPE SINGLE PRECISION (SINTDY-S, DINTDY-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C SINTDY computes interpolated values of the K-th derivative of the C dependent variable vector y, and stores it in DKY. This routine C is called within the package with K = 0 and T = TOUT, but may C also be called by the user for any K up to the current order. C (See detailed instructions in the usage documentation.) C C The computed values in DKY are gotten by interpolation using the C Nordsieck history array YH. This array corresponds uniquely to a C vector-valued polynomial of degree NQCUR or less, and DKY is set C to the K-th derivative of this polynomial at T. C The formula for DKY is: C q C DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1) C j=K C where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR. C The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are C communicated by COMMON. The above sum is done in reverse order. C IFLAG is returned negative if either K or T is out of bounds. C C***SEE ALSO SLSODE C***ROUTINES CALLED XERRWV C***COMMON BLOCKS SLS001 C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890501 Modified prologue to SLATEC/LDOC format. (FNF) C 890503 Minor cosmetic changes. (FNF) C 930809 Renamed to allow single/double precision versions. (ACH) C 010412 Reduced size of Common block /SLS001/. (ACH) C 031105 Restored 'own' variables to Common block /SLS001/, to C enable interrupt/restart feature. (ACH) C 050427 Corrected roundoff decrement in TP. (ACH) C***END PROLOGUE SINTDY C**End INTEGER K, NYH, IFLAG REAL T, YH, DKY DIMENSION YH(NYH,*), DKY(*) INTEGER IOWND, IOWNS, 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU REAL ROWNS, 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND COMMON /SLS001/ ROWNS(209), 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 2 IOWND(6), IOWNS(6), 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1 REAL C, R, S, TP CHARACTER*80 MSG C C***FIRST EXECUTABLE STATEMENT SINTDY IFLAG = 0 IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80 TP = TN - HU - 100.0E0*UROUND*SIGN(ABS(TN) + ABS(HU), HU) IF ((T-TP)*(T-TN) .GT. 0.0E0) GO TO 90 C S = (T - TN)/H IC = 1 IF (K .EQ. 0) GO TO 15 JJ1 = L - K DO 10 JJ = JJ1,NQ 10 IC = IC*JJ 15 C = IC DO 20 I = 1,N 20 DKY(I) = C*YH(I,L) IF (K .EQ. NQ) GO TO 55 JB2 = NQ - K DO 50 JB = 1,JB2 J = NQ - JB JP1 = J + 1 IC = 1 IF (K .EQ. 0) GO TO 35 JJ1 = JP1 - K DO 30 JJ = JJ1,J 30 IC = IC*JJ 35 C = IC DO 40 I = 1,N 40 DKY(I) = C*YH(I,JP1) + S*DKY(I) 50 CONTINUE IF (K .EQ. 0) RETURN 55 R = H**(-K) DO 60 I = 1,N 60 DKY(I) = R*DKY(I) RETURN C 80 MSG = 'SINTDY- K (=I1) illegal ' CALL XERRWV (MSG, 30, 51, 0, 1, K, 0, 0, 0.0E0, 0.0E0) IFLAG = -1 RETURN 90 MSG = 'SINTDY- T (=R1) illegal ' CALL XERRWV (MSG, 30, 52, 0, 0, 0, 0, 1, T, 0.0E0) MSG=' T not in interval TCUR - HU (= R1) to TCUR (=R2) ' CALL XERRWV (MSG, 60, 52, 0, 0, 0, 0, 2, TP, TN) IFLAG = -2 RETURN C----------------------- END OF SUBROUTINE SINTDY ---------------------- END *DECK SPREPJ SUBROUTINE SPREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, 1 F, JAC) C***BEGIN PROLOGUE SPREPJ C***SUBSIDIARY C***PURPOSE Compute and process Newton iteration matrix. C***TYPE SINGLE PRECISION (SPREPJ-S, DPREPJ-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C SPREPJ is called by SSTODE to compute and process the matrix C P = I - h*el(1)*J , where J is an approximation to the Jacobian. C Here J is computed by the user-supplied routine JAC if C MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5. C If MITER = 3, a diagonal approximation to J is used. C J is stored in WM and replaced by P. If MITER .ne. 3, P is then C subjected to LU decomposition in preparation for later solution C of linear systems with P as coefficient matrix. This is done C by SGEFA if MITER = 1 or 2, and by SGBFA if MITER = 4 or 5. C C In addition to variables described in SSTODE and SLSODE prologues, C communication with SPREPJ uses the following: C Y = array containing predicted values on entry. C FTEM = work array of length N (ACOR in SSTODE). C SAVF = array containing f evaluated at predicted y. C WM = real work space for matrices. On output it contains the C inverse diagonal matrix if MITER = 3 and the LU decomposition C of P if MITER is 1, 2 , 4, or 5. C Storage of matrix elements starts at WM(3). C WM also contains the following matrix-related data: C WM(1) = SQRT(UROUND), used in numerical Jacobian increments. C WM(2) = H*EL0, saved for later use if MITER = 3. C IWM = integer work space containing pivot information, starting at C IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band C parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. C EL0 = EL(1) (input). C IERPJ = output error flag, = 0 if no trouble, .gt. 0 if C P matrix found to be singular. C JCUR = output flag = 1 to indicate that the Jacobian matrix C (or approximation) is now current. C This routine also uses the COMMON variables EL0, H, TN, UROUND, C MITER, N, NFE, and NJE. C C***SEE ALSO SLSODE C***ROUTINES CALLED SGBFA, SGEFA, SVNORM C***COMMON BLOCKS SLS001 C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890501 Modified prologue to SLATEC/LDOC format. (FNF) C 890504 Minor cosmetic changes. (FNF) C 930809 Renamed to allow single/double precision versions. (ACH) C 010412 Reduced size of Common block /SLS001/. (ACH) C 031105 Restored 'own' variables to Common block /SLS001/, to C enable interrupt/restart feature. (ACH) C***END PROLOGUE SPREPJ C**End EXTERNAL F, JAC INTEGER NEQ, NYH, IWM REAL Y, YH, EWT, FTEM, SAVF, WM DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*), 1 WM(*), IWM(*) INTEGER IOWND, IOWNS, 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU REAL ROWNS, 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND COMMON /SLS001/ ROWNS(209), 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 2 IOWND(6), IOWNS(6), 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP, 1 MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1 REAL CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, 1 SVNORM C C***FIRST EXECUTABLE STATEMENT SPREPJ NJE = NJE + 1 IERPJ = 0 JCUR = 1 HL0 = H*EL0 GO TO (100, 200, 300, 400, 500), MITER C If MITER = 1, call JAC and multiply by scalar. ----------------------- 100 LENP = N*N DO 110 I = 1,LENP 110 WM(I+2) = 0.0E0 CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N) CON = -HL0 DO 120 I = 1,LENP 120 WM(I+2) = WM(I+2)*CON GO TO 240 C If MITER = 2, make N calls to F to approximate J. -------------------- 200 FAC = SVNORM (N, SAVF, EWT) R0 = 1000.0E0*ABS(H)*UROUND*N*FAC IF (R0 .EQ. 0.0E0) R0 = 1.0E0 SRUR = WM(1) J1 = 2 DO 230 J = 1,N YJ = Y(J) R = MAX(SRUR*ABS(YJ),R0/EWT(J)) Y(J) = Y(J) + R FAC = -HL0/R CALL F (NEQ, TN, Y, FTEM) DO 220 I = 1,N 220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC Y(J) = YJ J1 = J1 + N 230 CONTINUE NFE = NFE + N C Add identity matrix. ------------------------------------------------- 240 J = 3 NP1 = N + 1 DO 250 I = 1,N WM(J) = WM(J) + 1.0E0 250 J = J + NP1 C Do LU decomposition on P. -------------------------------------------- CALL SGEFA (WM(3), N, N, IWM(21), IER) IF (IER .NE. 0) IERPJ = 1 RETURN C If MITER = 3, construct a diagonal approximation to J and P. --------- 300 WM(2) = HL0 R = EL0*0.1E0 DO 310 I = 1,N 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) CALL F (NEQ, TN, Y, WM(3)) NFE = NFE + 1 DO 320 I = 1,N R0 = H*SAVF(I) - YH(I,2) DI = 0.1E0*R0 - H*(WM(I+2) - SAVF(I)) WM(I+2) = 1.0E0 IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320 IF (ABS(DI) .EQ. 0.0E0) GO TO 330 WM(I+2) = 0.1E0*R0/DI 320 CONTINUE RETURN 330 IERPJ = 1 RETURN C If MITER = 4, call JAC and multiply by scalar. ----------------------- 400 ML = IWM(1) MU = IWM(2) ML3 = ML + 3 MBAND = ML + MU + 1 MEBAND = MBAND + ML LENP = MEBAND*N DO 410 I = 1,LENP 410 WM(I+2) = 0.0E0 CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND) CON = -HL0 DO 420 I = 1,LENP 420 WM(I+2) = WM(I+2)*CON GO TO 570 C If MITER = 5, make MBAND calls to F to approximate J. ---------------- 500 ML = IWM(1) MU = IWM(2) MBAND = ML + MU + 1 MBA = MIN(MBAND,N) MEBAND = MBAND + ML MEB1 = MEBAND - 1 SRUR = WM(1) FAC = SVNORM (N, SAVF, EWT) R0 = 1000.0E0*ABS(H)*UROUND*N*FAC IF (R0 .EQ. 0.0E0) R0 = 1.0E0 DO 560 J = 1,MBA DO 530 I = J,N,MBAND YI = Y(I) R = MAX(SRUR*ABS(YI),R0/EWT(I)) 530 Y(I) = Y(I) + R CALL F (NEQ, TN, Y, FTEM) DO 550 JJ = J,N,MBAND Y(JJ) = YH(JJ,1) YJJ = Y(JJ) R = MAX(SRUR*ABS(YJJ),R0/EWT(JJ)) FAC = -HL0/R I1 = MAX(JJ-MU,1) I2 = MIN(JJ+ML,N) II = JJ*MEB1 - ML + 2 DO 540 I = I1,I2 540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC 550 CONTINUE 560 CONTINUE NFE = NFE + MBA C Add identity matrix. ------------------------------------------------- 570 II = MBAND + 2 DO 580 I = 1,N WM(II) = WM(II) + 1.0E0 580 II = II + MEBAND C Do LU decomposition of P. -------------------------------------------- CALL SGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER) IF (IER .NE. 0) IERPJ = 1 RETURN C----------------------- END OF SUBROUTINE SPREPJ ---------------------- END *DECK SSOLSY SUBROUTINE SSOLSY (WM, IWM, X, TEM) C***BEGIN PROLOGUE SSOLSY C***SUBSIDIARY C***PURPOSE ODEPACK linear system solver. C***TYPE SINGLE PRECISION (SSOLSY-S, DSOLSY-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C This routine manages the solution of the linear system arising from C a chord iteration. It is called if MITER .ne. 0. C If MITER is 1 or 2, it calls SGESL to accomplish this. C If MITER = 3 it updates the coefficient h*EL0 in the diagonal C matrix, and then computes the solution. C If MITER is 4 or 5, it calls SGBSL. C Communication with SSOLSY uses the following variables: C WM = real work space containing the inverse diagonal matrix if C MITER = 3 and the LU decomposition of the matrix otherwise. C Storage of matrix elements starts at WM(3). C WM also contains the following matrix-related data: C WM(1) = SQRT(UROUND) (not used here), C WM(2) = HL0, the previous value of h*EL0, used if MITER = 3. C IWM = integer work space containing pivot information, starting at C IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band C parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. C X = the right-hand side vector on input, and the solution vector C on output, of length N. C TEM = vector of work space of length N, not used in this version. C IERSL = output flag (in COMMON). IERSL = 0 if no trouble occurred. C IERSL = 1 if a singular matrix arose with MITER = 3. C This routine also uses the COMMON variables EL0, H, MITER, and N. C C***SEE ALSO SLSODE C***ROUTINES CALLED SGBSL, SGESL C***COMMON BLOCKS SLS001 C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890501 Modified prologue to SLATEC/LDOC format. (FNF) C 890503 Minor cosmetic changes. (FNF) C 930809 Renamed to allow single/double precision versions. (ACH) C 010412 Reduced size of Common block /SLS001/. (ACH) C 031105 Restored 'own' variables to Common block /SLS001/, to C enable interrupt/restart feature. (ACH) C***END PROLOGUE SSOLSY C**End INTEGER IWM REAL WM, X, TEM DIMENSION WM(*), IWM(*), X(*), TEM(*) INTEGER IOWND, IOWNS, 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU REAL ROWNS, 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND COMMON /SLS001/ ROWNS(209), 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 2 IOWND(6), IOWNS(6), 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU INTEGER I, MEBAND, ML, MU REAL DI, HL0, PHL0, R C C***FIRST EXECUTABLE STATEMENT SSOLSY IERSL = 0 GO TO (100, 100, 300, 400, 400), MITER 100 CALL SGESL (WM(3), N, N, IWM(21), X, 0) RETURN C 300 PHL0 = WM(2) HL0 = H*EL0 WM(2) = HL0 IF (HL0 .EQ. PHL0) GO TO 330 R = HL0/PHL0 DO 320 I = 1,N DI = 1.0E0 - R*(1.0E0 - 1.0E0/WM(I+2)) IF (ABS(DI) .EQ. 0.0E0) GO TO 390 320 WM(I+2) = 1.0E0/DI 330 DO 340 I = 1,N 340 X(I) = WM(I+2)*X(I) RETURN 390 IERSL = 1 RETURN C 400 ML = IWM(1) MU = IWM(2) MEBAND = 2*ML + MU + 1 CALL SGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0) RETURN C----------------------- END OF SUBROUTINE SSOLSY ---------------------- END *DECK SSRCOM SUBROUTINE SSRCOM (RSAV, ISAV, JOB) C***BEGIN PROLOGUE SSRCOM C***SUBSIDIARY C***PURPOSE Save/restore ODEPACK COMMON blocks. C***TYPE SINGLE PRECISION (SSRCOM-S, DSRCOM-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C This routine saves or restores (depending on JOB) the contents of C the COMMON block SLS001, which is used internally C by one or more ODEPACK solvers. C C RSAV = real array of length 218 or more. C ISAV = integer array of length 37 or more. C JOB = flag indicating to save or restore the COMMON blocks: C JOB = 1 if COMMON is to be saved (written to RSAV/ISAV) C JOB = 2 if COMMON is to be restored (read from RSAV/ISAV) C A call with JOB = 2 presumes a prior call with JOB = 1. C C***SEE ALSO SLSODE C***ROUTINES CALLED (NONE) C***COMMON BLOCKS SLS001 C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890501 Modified prologue to SLATEC/LDOC format. (FNF) C 890503 Minor cosmetic changes. (FNF) C 921116 Deleted treatment of block /EH0001/. (ACH) C 930801 Reduced Common block length by 2. (ACH) C 930809 Renamed to allow single/double precision versions. (ACH) C 010412 Reduced Common block length by 209+12. (ACH) C 031105 Restored 'own' variables to Common block /SLS001/, to C enable interrupt/restart feature. (ACH) C 031112 Added SAVE statement for data-loaded constants. C***END PROLOGUE SSRCOM C**End INTEGER ISAV, JOB INTEGER ILS INTEGER I, LENILS, LENRLS REAL RSAV, RLS DIMENSION RSAV(*), ISAV(*) SAVE LENRLS, LENILS COMMON /SLS001/ RLS(218), ILS(37) DATA LENRLS/218/, LENILS/37/ C C***FIRST EXECUTABLE STATEMENT SSRCOM IF (JOB .EQ. 2) GO TO 100 C DO 10 I = 1,LENRLS 10 RSAV(I) = RLS(I) DO 20 I = 1,LENILS 20 ISAV(I) = ILS(I) RETURN C 100 CONTINUE DO 110 I = 1,LENRLS 110 RLS(I) = RSAV(I) DO 120 I = 1,LENILS 120 ILS(I) = ISAV(I) RETURN C----------------------- END OF SUBROUTINE SSRCOM ---------------------- END *DECK SSTODE SUBROUTINE SSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, 1 WM, IWM, F, JAC, PJAC, SLVS) C***BEGIN PROLOGUE SSTODE C***SUBSIDIARY C***PURPOSE Performs one step of an ODEPACK integration. C***TYPE SINGLE PRECISION (SSTODE-S, DSTODE-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C SSTODE performs one step of the integration of an initial value C problem for a system of ordinary differential equations. C Note: SSTODE is independent of the value of the iteration method C indicator MITER, when this is .ne. 0, and hence is independent C of the type of chord method used, or the Jacobian structure. C Communication with SSTODE is done with the following variables: C C NEQ = integer array containing problem size in NEQ(1), and C passed as the NEQ argument in all calls to F and JAC. C Y = an array of length .ge. N used as the Y argument in C all calls to F and JAC. C YH = an NYH by LMAX array containing the dependent variables C and their approximate scaled derivatives, where C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate C j-th derivative of y(i), scaled by h**j/factorial(j) C (j = 0,1,...,NQ). on entry for the first step, the first C two columns of YH must be set from the initial values. C NYH = a constant integer .ge. N, the first dimension of YH. C YH1 = a one-dimensional array occupying the same space as YH. C EWT = an array of length N containing multiplicative weights C for local error measurements. Local errors in Y(i) are C compared to 1.0/EWT(i) in various error tests. C SAVF = an array of working storage, of length N. C Also used for input of YH(*,MAXORD+2) when JSTART = -1 C and MAXORD .lt. the current order NQ. C ACOR = a work array of length N, used for the accumulated C corrections. On a successful return, ACOR(i) contains C the estimated one-step local error in Y(i). C WM,IWM = real and integer work arrays associated with matrix C operations in chord iteration (MITER .ne. 0). C PJAC = name of routine to evaluate and preprocess Jacobian matrix C and P = I - h*el0*JAC, if a chord method is being used. C SLVS = name of routine to solve linear system in chord iteration. C CCMAX = maximum relative change in h*el0 before PJAC is called. 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 absolute value of the step size h to be used. C HMXI = inverse of the maximum absolute value of h to be used. C HMXI = 0.0 is allowed and corresponds to an infinite hmax. C HMIN and HMXI may be changed at any time, but will not C take effect until the next change of h is considered. C TN = the independent variable. TN is updated on each step taken. C JSTART = an integer used for input only, with the following C values and meanings: C 0 perform the first step. C .gt.0 take a new step continuing from the last. C -1 take the next step with a new value of H, MAXORD, C N, METH, MITER, and/or matrix parameters. C -2 take the next step with a new value of H, C but with other inputs unchanged. C On return, JSTART is set to 1 to facilitate continuation. 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 -2 corrector convergence could not be achieved. C -3 fatal error in PJAC or SLVS. C A return with KFLAG = -1 or -2 means either C abs(H) = HMIN or 10 consecutive failures occurred. C On a return with KFLAG negative, the values of TN and C the YH array are as of the beginning of the last C step, and H is the last step size attempted. C MAXORD = the maximum order of integration method to be allowed. C MAXCOR = the maximum number of corrector iterations allowed. C MSBP = maximum number of steps between PJAC calls (MITER .gt. 0). C MXNCF = maximum number of convergence failures allowed. C METH/MITER = the method flags. See description in driver. C N = the number of first-order differential equations. C The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD, C MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON. C C***SEE ALSO SLSODE C***ROUTINES CALLED SCFODE, SVNORM C***COMMON BLOCKS SLS001 C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890501 Modified prologue to SLATEC/LDOC format. (FNF) C 890503 Minor cosmetic changes. (FNF) C 930809 Renamed to allow single/double precision versions. (ACH) C 010413 Reduced size of Common block /SLS001/. (ACH) C 031105 Restored 'own' variables to Common block /SLS001/, to C enable interrupt/restart feature. (ACH) C***END PROLOGUE SSTODE C**End EXTERNAL F, JAC, PJAC, SLVS INTEGER NEQ, NYH, IWM REAL Y, YH, YH1, EWT, SAVF, ACOR, WM DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*), 1 ACOR(*), WM(*), IWM(*) INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, 1 R, RH, RHDN, RHSM, RHUP, TOLD, SVNORM COMMON /SLS001/ CONIT, CRATE, EL(13), ELCO(13,12), 1 HOLD, RMAX, TESCO(3,12), 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 3 IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU C C***FIRST EXECUTABLE STATEMENT SSTODE KFLAG = 0 TOLD = TN NCF = 0 IERPJ = 0 IERSL = 0 JCUR = 0 ICF = 0 DELP = 0.0E0 IF (JSTART .GT. 0) GO TO 200 IF (JSTART .EQ. -1) GO TO 100 IF (JSTART .EQ. -2) GO TO 160 C----------------------------------------------------------------------- C On the first call, the order is set to 1, and other variables are C initialized. 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 to 2 C for the next increase. C----------------------------------------------------------------------- LMAX = MAXORD + 1 NQ = 1 L = 2 IALTH = 2 RMAX = 10000.0E0 RC = 0.0E0 EL0 = 1.0E0 CRATE = 0.7E0 HOLD = H MEO = METH NSLP = 0 IPUP = MITER IRET = 3 GO TO 140 C----------------------------------------------------------------------- C The following block handles preliminaries needed when JSTART = -1. C IPUP is set to MITER to force a matrix update. C If an order increase is about to be considered (IALTH = 1), C IALTH is reset to 2 to postpone consideration one more step. C If the caller has changed METH, SCFODE is called to reset C the coefficients of the method. C If the caller has changed MAXORD to a value less than the current C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly. C If H is to be changed, YH must be rescaled. C If H or METH is being changed, IALTH is reset to L = NQ + 1 C to prevent further changes in H for that many steps. C----------------------------------------------------------------------- 100 IPUP = MITER LMAX = MAXORD + 1 IF (IALTH .EQ. 1) IALTH = 2 IF (METH .EQ. MEO) GO TO 110 CALL SCFODE (METH, ELCO, TESCO) MEO = METH IF (NQ .GT. MAXORD) GO TO 120 IALTH = L IRET = 1 GO TO 150 110 IF (NQ .LE. MAXORD) GO TO 160 120 NQ = MAXORD L = LMAX DO 125 I = 1,L 125 EL(I) = ELCO(I,NQ) NQNYH = NQ*NYH RC = RC*EL(1)/EL0 EL0 = EL(1) CONIT = 0.5E0/(NQ+2) DDN = SVNORM (N, SAVF, EWT)/TESCO(1,L) EXDN = 1.0E0/L RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) RH = MIN(RHDN,1.0E0) IREDO = 3 IF (H .EQ. HOLD) GO TO 170 RH = MIN(RH,ABS(H/HOLD)) H = HOLD GO TO 175 C----------------------------------------------------------------------- C SCFODE is called to get all the integration coefficients for the C current METH. Then the EL vector and related constants are reset C whenever the order NQ is changed, or at the start of the problem. C----------------------------------------------------------------------- 140 CALL SCFODE (METH, ELCO, TESCO) 150 DO 155 I = 1,L 155 EL(I) = ELCO(I,NQ) NQNYH = NQ*NYH RC = RC*EL(1)/EL0 EL0 = EL(1) CONIT = 0.5E0/(NQ+2) GO TO (160, 170, 200), IRET C----------------------------------------------------------------------- C If H is being changed, the H ratio RH is checked against C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to C L = NQ + 1 to prevent a change of H for that many steps, unless C forced by a convergence or error test failure. C----------------------------------------------------------------------- 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,RMAX) RH = RH/MAX(1.0E0,ABS(H)*HMXI*RH) R = 1.0E0 DO 180 J = 2,L R = R*RH DO 180 I = 1,N 180 YH(I,J) = YH(I,J)*R H = H*RH RC = RC*RH IALTH = L IF (IREDO .EQ. 0) GO TO 690 C----------------------------------------------------------------------- C This section computes the predicted values by effectively C multiplying the YH 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 CCMAX, IPUP is set to MITER C to force PJAC to be called, if a Jacobian is involved. C In any case, PJAC is called at least every MSBP steps. C----------------------------------------------------------------------- 200 IF (ABS(RC-1.0E0) .GT. CCMAX) IPUP = MITER IF (NST .GE. NSLP+MSBP) IPUP = MITER TN = TN + H I1 = NQNYH + 1 DO 215 JB = 1,NQ I1 = I1 - NYH Cdir$ ivdep DO 210 I = I1,NQNYH 210 YH1(I) = YH1(I) + YH1(I+NYH) 215 CONTINUE C----------------------------------------------------------------------- C Up to MAXCOR corrector iterations are taken. A convergence test is C made on the R.M.S. norm of each correction, weighted by the error C weight vector EWT. The sum of the corrections is accumulated in the C vector ACOR(i). The YH array is not altered in the corrector loop. C----------------------------------------------------------------------- 220 M = 0 DO 230 I = 1,N 230 Y(I) = YH(I,1) CALL F (NEQ, TN, Y, SAVF) NFE = NFE + 1 IF (IPUP .LE. 0) GO TO 250 C----------------------------------------------------------------------- C If indicated, the matrix P = I - h*el(1)*J is reevaluated and C preprocessed before starting the corrector iteration. IPUP is set C to 0 as an indicator that this has been done. C----------------------------------------------------------------------- CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC) IPUP = 0 RC = 1.0E0 NSLP = NST CRATE = 0.7E0 IF (IERPJ .NE. 0) GO TO 430 250 DO 260 I = 1,N 260 ACOR(I) = 0.0E0 270 IF (MITER .NE. 0) GO TO 350 C----------------------------------------------------------------------- C In the case of functional iteration, update Y directly from C the result of the last function evaluation. C----------------------------------------------------------------------- DO 290 I = 1,N SAVF(I) = H*SAVF(I) - YH(I,2) 290 Y(I) = SAVF(I) - ACOR(I) DEL = SVNORM (N, Y, EWT) DO 300 I = 1,N Y(I) = YH(I,1) + EL(1)*SAVF(I) 300 ACOR(I) = SAVF(I) GO TO 400 C----------------------------------------------------------------------- C In the case of the chord method, compute the corrector error, C and solve the linear system with that as right-hand side and C P as coefficient matrix. C----------------------------------------------------------------------- 350 DO 360 I = 1,N 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) CALL SLVS (WM, IWM, Y, SAVF) IF (IERSL .LT. 0) GO TO 430 IF (IERSL .GT. 0) GO TO 410 DEL = SVNORM (N, Y, EWT) DO 380 I = 1,N ACOR(I) = ACOR(I) + Y(I) 380 Y(I) = YH(I,1) + EL(1)*ACOR(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(0.2E0*CRATE,DEL/DELP) DCON = DEL*MIN(1.0E0,1.5E0*CRATE)/(TESCO(2,NQ)*CONIT) IF (DCON .LE. 1.0E0) GO TO 450 M = M + 1 IF (M .EQ. MAXCOR) GO TO 410 IF (M .GE. 2 .AND. DEL .GT. 2.0E0*DELP) GO TO 410 DELP = DEL CALL F (NEQ, TN, Y, SAVF) NFE = NFE + 1 GO TO 270 C----------------------------------------------------------------------- C The corrector iteration failed to converge. C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for C the next try. Otherwise the YH array is retracted to its values C before prediction, and H is reduced, if possible. If H cannot be C reduced or MXNCF failures have occurred, exit with KFLAG = -2. C----------------------------------------------------------------------- 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430 ICF = 1 IPUP = MITER GO TO 220 430 ICF = 2 NCF = NCF + 1 RMAX = 2.0E0 TN = TOLD I1 = NQNYH + 1 DO 445 JB = 1,NQ I1 = I1 - NYH Cdir$ ivdep DO 440 I = I1,NQNYH 440 YH1(I) = YH1(I) - YH1(I+NYH) 445 CONTINUE IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 670 IF (NCF .EQ. MXNCF) GO TO 670 RH = 0.25E0 IPUP = MITER IREDO = 1 GO TO 170 C----------------------------------------------------------------------- C The corrector has converged. JCUR is set to 0 C to signal that the Jacobian involved may need updating later. C The local error test is made and control passes to statement 500 C if it fails. C----------------------------------------------------------------------- 450 JCUR = 0 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) IF (M .GT. 0) DSM = SVNORM (N, ACOR, EWT)/TESCO(2,NQ) IF (DSM .GT. 1.0E0) GO TO 500 C----------------------------------------------------------------------- C After a successful step, update the YH array. C Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1. C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR 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, IALTH is set to 3 to prevent C testing for that many steps. C----------------------------------------------------------------------- KFLAG = 0 IREDO = 0 NST = NST + 1 HU = H NQU = NQ DO 470 J = 1,L DO 470 I = 1,N 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) IALTH = IALTH - 1 IF (IALTH .EQ. 0) GO TO 520 IF (IALTH .GT. 1) GO TO 700 IF (L .EQ. LMAX) GO TO 700 DO 490 I = 1,N 490 YH(I,LMAX) = ACOR(I) GO TO 700 C----------------------------------------------------------------------- C The error test failed. KFLAG keeps track of multiple failures. C Restore TN and the YH 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. After 2 or more failures, H is forced to decrease C by a factor of 0.2 or less. C----------------------------------------------------------------------- 500 KFLAG = KFLAG - 1 TN = TOLD I1 = NQNYH + 1 DO 515 JB = 1,NQ I1 = I1 - NYH Cdir$ ivdep DO 510 I = I1,NQNYH 510 YH1(I) = YH1(I) - YH1(I+NYH) 515 CONTINUE RMAX = 2.0E0 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 660 IF (KFLAG .LE. -3) GO TO 640 IREDO = 2 RHUP = 0.0E0 GO TO 540 C----------------------------------------------------------------------- C Regardless of the success or failure of the step, factors C RHDN, RHSM, and RHUP are computed, by which H could be multiplied C at order NQ - 1, order NQ, or order NQ + 1, respectively. C In the case of failure, RHUP = 0.0 to avoid an order increase. C The largest 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 RHUP = 0.0E0 IF (L .EQ. LMAX) GO TO 540 DO 530 I = 1,N 530 SAVF(I) = ACOR(I) - YH(I,LMAX) DUP = SVNORM (N, SAVF, EWT)/TESCO(3,NQ) EXUP = 1.0E0/(L+1) RHUP = 1.0E0/(1.4E0*DUP**EXUP + 0.0000014E0) 540 EXSM = 1.0E0/L RHSM = 1.0E0/(1.2E0*DSM**EXSM + 0.0000012E0) RHDN = 0.0E0 IF (NQ .EQ. 1) GO TO 560 DDN = SVNORM (N, YH(1,L), EWT)/TESCO(1,NQ) EXDN = 1.0E0/NQ RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) 560 IF (RHSM .GE. RHUP) GO TO 570 IF (RHUP .GT. RHDN) GO TO 590 GO TO 580 570 IF (RHSM .LT. RHDN) GO TO 580 NEWQ = NQ RH = RHSM GO TO 620 580 NEWQ = NQ - 1 RH = RHDN IF (KFLAG .LT. 0 .AND. RH .GT. 1.0E0) RH = 1.0E0 GO TO 620 590 NEWQ = L RH = RHUP IF (RH .LT. 1.1E0) GO TO 610 R = EL(L)/L DO 600 I = 1,N 600 YH(I,NEWQ+1) = ACOR(I)*R GO TO 630 610 IALTH = 3 GO TO 700 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1E0)) GO TO 610 IF (KFLAG .LE. -2) RH = MIN(RH,0.2E0) 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 YH 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 150 C----------------------------------------------------------------------- C Control reaches this section if 3 or more failures have occured. C If 10 failures have occurred, exit with KFLAG = -1. C It is assumed that the derivatives that have accumulated in the C YH 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 until it succeeds or H reaches HMIN. C----------------------------------------------------------------------- 640 IF (KFLAG .EQ. -10) GO TO 660 RH = 0.1E0 RH = MAX(HMIN/ABS(H),RH) H = H*RH DO 645 I = 1,N 645 Y(I) = YH(I,1) CALL F (NEQ, TN, Y, SAVF) NFE = NFE + 1 DO 650 I = 1,N 650 YH(I,2) = H*SAVF(I) IPUP = MITER IALTH = 5 IF (NQ .EQ. 1) GO TO 200 NQ = 1 L = 2 IRET = 3 GO TO 150 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 720 670 KFLAG = -2 GO TO 720 680 KFLAG = -3 GO TO 720 690 RMAX = 10.0E0 700 R = 1.0E0/TESCO(2,NQU) DO 710 I = 1,N 710 ACOR(I) = ACOR(I)*R 720 HOLD = H JSTART = 1 RETURN C----------------------- END OF SUBROUTINE SSTODE ---------------------- END *DECK SEWSET SUBROUTINE SEWSET (N, ITOL, RTOL, ATOL, YCUR, EWT) C***BEGIN PROLOGUE SEWSET C***SUBSIDIARY C***PURPOSE Set error weight vector. C***TYPE SINGLE PRECISION (SEWSET-S, DEWSET-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C This subroutine sets the error weight vector EWT according to C EWT(i) = RTOL(i)*ABS(YCUR(i)) + ATOL(i), i = 1,...,N, C with the subscript on RTOL and/or ATOL possibly replaced by 1 above, C depending on the value of ITOL. C C***SEE ALSO SLSODE C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890501 Modified prologue to SLATEC/LDOC format. (FNF) C 890503 Minor cosmetic changes. (FNF) C 930809 Renamed to allow single/double precision versions. (ACH) C***END PROLOGUE SEWSET C**End INTEGER N, ITOL INTEGER I REAL RTOL, ATOL, YCUR, EWT DIMENSION RTOL(*), ATOL(*), YCUR(N), EWT(N) C C***FIRST EXECUTABLE STATEMENT SEWSET GO TO (10, 20, 30, 40), ITOL 10 CONTINUE DO 15 I = 1,N 15 EWT(I) = RTOL(1)*ABS(YCUR(I)) + ATOL(1) RETURN 20 CONTINUE DO 25 I = 1,N 25 EWT(I) = RTOL(1)*ABS(YCUR(I)) + ATOL(I) RETURN 30 CONTINUE DO 35 I = 1,N 35 EWT(I) = RTOL(I)*ABS(YCUR(I)) + ATOL(1) RETURN 40 CONTINUE DO 45 I = 1,N 45 EWT(I) = RTOL(I)*ABS(YCUR(I)) + ATOL(I) RETURN C----------------------- END OF SUBROUTINE SEWSET ---------------------- END *DECK SVNORM REAL FUNCTION SVNORM (N, V, W) C***BEGIN PROLOGUE SVNORM C***SUBSIDIARY C***PURPOSE Weighted root-mean-square vector norm. C***TYPE SINGLE PRECISION (SVNORM-S, DVNORM-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C This function routine computes the weighted root-mean-square norm C of the vector of length N contained in the array V, with weights C contained in the array W of length N: C SVNORM = SQRT( (1/N) * SUM( V(i)*W(i) )**2 ) C C***SEE ALSO SLSODE C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890501 Modified prologue to SLATEC/LDOC format. (FNF) C 890503 Minor cosmetic changes. (FNF) C 930809 Renamed to allow single/double precision versions. (ACH) C***END PROLOGUE SVNORM C**End INTEGER N, I REAL V, W, SUM DIMENSION V(N), W(N) C C***FIRST EXECUTABLE STATEMENT SVNORM SUM = 0.0E0 DO 10 I = 1,N 10 SUM = SUM + (V(I)*W(I))**2 SVNORM = SQRT(SUM/N) RETURN C----------------------- END OF FUNCTION SVNORM ------------------------ END *DECK SIPREP SUBROUTINE SIPREP (NEQ, Y, RWORK, IA, JA, IPFLAG, F, JAC) EXTERNAL F, JAC INTEGER NEQ, IA, JA, IPFLAG REAL Y, RWORK DIMENSION NEQ(*), Y(*), RWORK(*), IA(*), JA(*) INTEGER IOWND, IOWNS, 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP, 1 IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA, 2 LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ, 3 NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU REAL ROWNS, 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND REAL RLSS COMMON /SLS001/ ROWNS(209), 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 2 IOWND(6), IOWNS(6), 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU COMMON /SLSS01/ RLSS(6), 1 IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP, 2 IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA, 3 LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ, 4 NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU INTEGER I, IMAX, LEWTN, LYHD, LYHN C----------------------------------------------------------------------- C This routine serves as an interface between the driver and C Subroutine SPREP. It is called only if MITER is 1 or 2. C Tasks performed here are: C * call SPREP, C * reset the required WM segment length LENWK, C * move YH back to its final location (following WM in RWORK), C * reset pointers for YH, SAVF, EWT, and ACOR, and C * move EWT to its new position if ISTATE = 1. C IPFLAG is an output error indication flag. IPFLAG = 0 if there was C no trouble, and IPFLAG is the value of the SPREP error flag IPPER C if there was trouble in Subroutine SPREP. C----------------------------------------------------------------------- IPFLAG = 0 C Call SPREP to do matrix preprocessing operations. -------------------- CALL SPREP (NEQ, Y, RWORK(LYH), RWORK(LSAVF), RWORK(LEWT), 1 RWORK(LACOR), IA, JA, RWORK(LWM), RWORK(LWM), IPFLAG, F, JAC) LENWK = MAX(LREQ,LWMIN) IF (IPFLAG .LT. 0) RETURN C If SPREP was successful, move YH to end of required space for WM. ---- LYHN = LWM + LENWK IF (LYHN .GT. LYH) RETURN LYHD = LYH - LYHN IF (LYHD .EQ. 0) GO TO 20 IMAX = LYHN - 1 + LENYHM DO 10 I = LYHN,IMAX 10 RWORK(I) = RWORK(I+LYHD) LYH = LYHN C Reset pointers for SAVF, EWT, and ACOR. ------------------------------ 20 LSAVF = LYH + LENYH LEWTN = LSAVF + N LACOR = LEWTN + N IF (ISTATC .EQ. 3) GO TO 40 C If ISTATE = 1, move EWT (left) to its new position. ------------------ IF (LEWTN .GT. LEWT) RETURN DO 30 I = 1,N 30 RWORK(I+LEWTN-1) = RWORK(I+LEWT-1) 40 LEWT = LEWTN RETURN C----------------------- End of Subroutine SIPREP ---------------------- END *DECK SPREP SUBROUTINE SPREP (NEQ, Y, YH, SAVF, EWT, FTEM, IA, JA, 1 WK, IWK, IPPER, F, JAC) EXTERNAL F,JAC INTEGER NEQ, IA, JA, IWK, IPPER REAL Y, YH, SAVF, EWT, FTEM, WK DIMENSION NEQ(*), Y(*), YH(*), SAVF(*), EWT(*), FTEM(*), 1 IA(*), JA(*), WK(*), IWK(*) INTEGER IOWND, IOWNS, 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP, 1 IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA, 2 LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ, 3 NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU REAL ROWNS, 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND REAL CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH COMMON /SLS001/ ROWNS(209), 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 2 IOWND(6), IOWNS(6), 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU COMMON /SLSS01/ CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH, 1 IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP, 2 IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA, 3 LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ, 4 NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU INTEGER I, IBR, IER, IPIL, IPIU, IPTT1, IPTT2, J, JFOUND, K, 1 KNEW, KMAX, KMIN, LDIF, LENIGP, LIWK, MAXG, NP1, NZSUT REAL DQ, DYJ, ERWT, FAC, YJ C----------------------------------------------------------------------- C This routine performs preprocessing related to the sparse linear C systems that must be solved if MITER = 1 or 2. C The operations that are performed here are: C * compute sparseness structure of Jacobian according to MOSS, C * compute grouping of column indices (MITER = 2), C * compute a new ordering of rows and columns of the matrix, C * reorder JA corresponding to the new ordering, C * perform a symbolic LU factorization of the matrix, and C * set pointers for segments of the IWK/WK array. C In addition to variables described previously, SPREP uses the C following for communication: C YH = the history array. Only the first column, containing the C current Y vector, is used. Used only if MOSS .ne. 0. C SAVF = a work array of length NEQ, used only if MOSS .ne. 0. C EWT = array of length NEQ containing (inverted) error weights. C Used only if MOSS = 2 or if ISTATE = MOSS = 1. C FTEM = a work array of length NEQ, identical to ACOR in the driver, C used only if MOSS = 2. C WK = a real work array of length LENWK, identical to WM in C the driver. C IWK = integer work array, assumed to occupy the same space as WK. C LENWK = the length of the work arrays WK and IWK. C ISTATC = a copy of the driver input argument ISTATE (= 1 on the C first call, = 3 on a continuation call). C IYS = flag value from SODRV or SCDRV. C IPPER = output error flag with the following values and meanings: C 0 no error. C -1 insufficient storage for internal structure pointers. C -2 insufficient storage for JGROUP. C -3 insufficient storage for SODRV. C -4 other error flag from SODRV (should never occur). C -5 insufficient storage for SCDRV. C -6 other error flag from SCDRV. C----------------------------------------------------------------------- IBIAN = LRAT*2 IPIAN = IBIAN + 1 NP1 = N + 1 IPJAN = IPIAN + NP1 IBJAN = IPJAN - 1 LIWK = LENWK*LRAT IF (IPJAN+N-1 .GT. LIWK) GO TO 210 IF (MOSS .EQ. 0) GO TO 30 C IF (ISTATC .EQ. 3) GO TO 20 C ISTATE = 1 and MOSS .ne. 0. Perturb Y for structure determination. -- DO 10 I = 1,N ERWT = 1.0E0/EWT(I) FAC = 1.0E0 + 1.0E0/(I + 1.0E0) Y(I) = Y(I) + FAC*SIGN(ERWT,Y(I)) 10 CONTINUE GO TO (70, 100), MOSS C 20 CONTINUE C ISTATE = 3 and MOSS .ne. 0. Load Y from YH(*,1). -------------------- DO 25 I = 1,N 25 Y(I) = YH(I) GO TO (70, 100), MOSS C C MOSS = 0. Process user's IA,JA. Add diagonal entries if necessary. - 30 KNEW = IPJAN KMIN = IA(1) IWK(IPIAN) = 1 DO 60 J = 1,N JFOUND = 0 KMAX = IA(J+1) - 1 IF (KMIN .GT. KMAX) GO TO 45 DO 40 K = KMIN,KMAX I = JA(K) IF (I .EQ. J) JFOUND = 1 IF (KNEW .GT. LIWK) GO TO 210 IWK(KNEW) = I KNEW = KNEW + 1 40 CONTINUE IF (JFOUND .EQ. 1) GO TO 50 45 IF (KNEW .GT. LIWK) GO TO 210 IWK(KNEW) = J KNEW = KNEW + 1 50 IWK(IPIAN+J) = KNEW + 1 - IPJAN KMIN = KMAX + 1 60 CONTINUE GO TO 140 C C MOSS = 1. Compute structure from user-supplied Jacobian routine JAC. 70 CONTINUE C A dummy call to F allows user to create temporaries for use in JAC. -- CALL F (NEQ, TN, Y, SAVF) K = IPJAN IWK(IPIAN) = 1 DO 90 J = 1,N IF (K .GT. LIWK) GO TO 210 IWK(K) = J K = K + 1 DO 75 I = 1,N 75 SAVF(I) = 0.0E0 CALL JAC (NEQ, TN, Y, J, IWK(IPIAN), IWK(IPJAN), SAVF) DO 80 I = 1,N IF (ABS(SAVF(I)) .LE. SETH) GO TO 80 IF (I .EQ. J) GO TO 80 IF (K .GT. LIWK) GO TO 210 IWK(K) = I K = K + 1 80 CONTINUE IWK(IPIAN+J) = K + 1 - IPJAN 90 CONTINUE GO TO 140 C C MOSS = 2. Compute structure from results of N + 1 calls to F. ------- 100 K = IPJAN IWK(IPIAN) = 1 CALL F (NEQ, TN, Y, SAVF) DO 120 J = 1,N IF (K .GT. LIWK) GO TO 210 IWK(K) = J K = K + 1 YJ = Y(J) ERWT = 1.0E0/EWT(J) DYJ = SIGN(ERWT,YJ) Y(J) = YJ + DYJ CALL F (NEQ, TN, Y, FTEM) Y(J) = YJ DO 110 I = 1,N DQ = (FTEM(I) - SAVF(I))/DYJ IF (ABS(DQ) .LE. SETH) GO TO 110 IF (I .EQ. J) GO TO 110 IF (K .GT. LIWK) GO TO 210 IWK(K) = I K = K + 1 110 CONTINUE IWK(IPIAN+J) = K + 1 - IPJAN 120 CONTINUE C 140 CONTINUE IF (MOSS .EQ. 0 .OR. ISTATC .NE. 1) GO TO 150 C If ISTATE = 1 and MOSS .ne. 0, restore Y from YH. -------------------- DO 145 I = 1,N 145 Y(I) = YH(I) 150 NNZ = IWK(IPIAN+N) - 1 LENIGP = 0 IPIGP = IPJAN + NNZ IF (MITER .NE. 2) GO TO 160 C C Compute grouping of column indices (MITER = 2). ---------------------- MAXG = NP1 IPJGP = IPJAN + NNZ IBJGP = IPJGP - 1 IPIGP = IPJGP + N IPTT1 = IPIGP + NP1 IPTT2 = IPTT1 + N LREQ = IPTT2 + N - 1 IF (LREQ .GT. LIWK) GO TO 220 CALL JGROUP (N, IWK(IPIAN), IWK(IPJAN), MAXG, NGP, IWK(IPIGP), 1 IWK(IPJGP), IWK(IPTT1), IWK(IPTT2), IER) IF (IER .NE. 0) GO TO 220 LENIGP = NGP + 1 C C Compute new ordering of rows/columns of Jacobian. -------------------- 160 IPR = IPIGP + LENIGP IPC = IPR IPIC = IPC + N IPISP = IPIC + N IPRSP = (IPISP - 2)/LRAT + 2 IESP = LENWK + 1 - IPRSP IF (IESP .LT. 0) GO TO 230 IBR = IPR - 1 DO 170 I = 1,N 170 IWK(IBR+I) = I NSP = LIWK + 1 - IPISP CALL SODRV (N, IWK(IPIAN), IWK(IPJAN), WK, IWK(IPR), IWK(IPIC), 1 NSP, IWK(IPISP), 1, IYS) IF (IYS .EQ. 11*N+1) GO TO 240 IF (IYS .NE. 0) GO TO 230 C C Reorder JAN and do symbolic LU factorization of matrix. -------------- IPA = LENWK + 1 - NNZ NSP = IPA - IPRSP LREQ = MAX(12*N/LRAT, 6*N/LRAT+2*N+NNZ) + 3 LREQ = LREQ + IPRSP - 1 + NNZ IF (LREQ .GT. LENWK) GO TO 250 IBA = IPA - 1 DO 180 I = 1,NNZ 180 WK(IBA+I) = 0.0E0 IPISP = LRAT*(IPRSP - 1) + 1 CALL SCDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN), 1 WK(IPA),WK(IPA),WK(IPA),NSP,IWK(IPISP),WK(IPRSP),IESP,5,IYS) LREQ = LENWK - IESP IF (IYS .EQ. 10*N+1) GO TO 250 IF (IYS .NE. 0) GO TO 260 IPIL = IPISP IPIU = IPIL + 2*N + 1 NZU = IWK(IPIL+N) - IWK(IPIL) NZL = IWK(IPIU+N) - IWK(IPIU) IF (LRAT .GT. 1) GO TO 190 CALL ADJLR (N, IWK(IPISP), LDIF) LREQ = LREQ + LDIF 190 CONTINUE IF (LRAT .EQ. 2 .AND. NNZ .EQ. N) LREQ = LREQ + 1 NSP = NSP + LREQ - LENWK IPA = LREQ + 1 - NNZ IBA = IPA - 1 IPPER = 0 RETURN C 210 IPPER = -1 LREQ = 2 + (2*N + 1)/LRAT LREQ = MAX(LENWK+1,LREQ) RETURN C 220 IPPER = -2 LREQ = (LREQ - 1)/LRAT + 1 RETURN C 230 IPPER = -3 CALL CNTNZU (N, IWK(IPIAN), IWK(IPJAN), NZSUT) LREQ = LENWK - IESP + (3*N + 4*NZSUT - 1)/LRAT + 1 RETURN C 240 IPPER = -4 RETURN C 250 IPPER = -5 RETURN C 260 IPPER = -6 LREQ = LENWK RETURN C----------------------- End of Subroutine SPREP ----------------------- END *DECK JGROUP SUBROUTINE JGROUP (N,IA,JA,MAXG,NGRP,IGP,JGP,INCL,JDONE,IER) INTEGER N, IA, JA, MAXG, NGRP, IGP, JGP, INCL, JDONE, IER DIMENSION IA(*), JA(*), IGP(*), JGP(*), INCL(*), JDONE(*) C----------------------------------------------------------------------- C This subroutine constructs groupings of the column indices of C the Jacobian matrix, used in the numerical evaluation of the C Jacobian by finite differences. C C Input: C N = the order of the matrix. C IA,JA = sparse structure descriptors of the matrix by rows. C MAXG = length of available storage in the IGP array. C C Output: C NGRP = number of groups. C JGP = array of length N containing the column indices by groups. C IGP = pointer array of length NGRP + 1 to the locations in JGP C of the beginning of each group. C IER = error indicator. IER = 0 if no error occurred, or 1 if C MAXG was insufficient. C C INCL and JDONE are working arrays of length N. C----------------------------------------------------------------------- INTEGER I, J, K, KMIN, KMAX, NCOL, NG C IER = 0 DO 10 J = 1,N 10 JDONE(J) = 0 NCOL = 1 DO 60 NG = 1,MAXG IGP(NG) = NCOL DO 20 I = 1,N 20 INCL(I) = 0 DO 50 J = 1,N C Reject column J if it is already in a group.-------------------------- IF (JDONE(J) .EQ. 1) GO TO 50 KMIN = IA(J) KMAX = IA(J+1) - 1 DO 30 K = KMIN,KMAX C Reject column J if it overlaps any column already in this group.------ I = JA(K) IF (INCL(I) .EQ. 1) GO TO 50 30 CONTINUE C Accept column J into group NG.---------------------------------------- JGP(NCOL) = J NCOL = NCOL + 1 JDONE(J) = 1 DO 40 K = KMIN,KMAX I = JA(K) 40 INCL(I) = 1 50 CONTINUE C Stop if this group is empty (grouping is complete).------------------- IF (NCOL .EQ. IGP(NG)) GO TO 70 60 CONTINUE C Error return if not all columns were chosen (MAXG too small).--------- IF (NCOL .LE. N) GO TO 80 NG = MAXG 70 NGRP = NG - 1 RETURN 80 IER = 1 RETURN C----------------------- End of Subroutine JGROUP ---------------------- END *DECK ADJLR SUBROUTINE ADJLR (N, ISP, LDIF) INTEGER N, ISP, LDIF DIMENSION ISP(*) C----------------------------------------------------------------------- C This routine computes an adjustment, LDIF, to the required C integer storage space in IWK (sparse matrix work space). C It is called only if the word length ratio is LRAT = 1. C This is to account for the possibility that the symbolic LU phase C may require more storage than the numerical LU and solution phases. C----------------------------------------------------------------------- INTEGER IP, JLMAX, JUMAX, LNFC, LSFC, NZLU C IP = 2*N + 1 C Get JLMAX = IJL(N) and JUMAX = IJU(N) (sizes of JL and JU). ---------- JLMAX = ISP(IP) JUMAX = ISP(IP+IP) C NZLU = (size of L) + (size of U) = (IL(N+1)-IL(1)) + (IU(N+1)-IU(1)). NZLU = ISP(N+1) - ISP(1) + ISP(IP+N+1) - ISP(IP+1) LSFC = 12*N + 3 + 2*MAX(JLMAX,JUMAX) LNFC = 9*N + 2 + JLMAX + JUMAX + NZLU LDIF = MAX(0, LSFC - LNFC) RETURN C----------------------- End of Subroutine ADJLR ----------------------- END *DECK CNTNZU SUBROUTINE CNTNZU (N, IA, JA, NZSUT) INTEGER N, IA, JA, NZSUT DIMENSION IA(*), JA(*) C----------------------------------------------------------------------- C This routine counts the number of nonzero elements in the strict C upper triangle of the matrix M + M(transpose), where the sparsity C structure of M is given by pointer arrays IA and JA. C This is needed to compute the storage requirements for the C sparse matrix reordering operation in SODRV. C----------------------------------------------------------------------- INTEGER II, JJ, J, JMIN, JMAX, K, KMIN, KMAX, NUM C NUM = 0 DO 50 II = 1,N JMIN = IA(II) JMAX = IA(II+1) - 1 IF (JMIN .GT. JMAX) GO TO 50 DO 40 J = JMIN,JMAX IF (JA(J) - II) 10, 40, 30 10 JJ =JA(J) KMIN = IA(JJ) KMAX = IA(JJ+1) - 1 IF (KMIN .GT. KMAX) GO TO 30 DO 20 K = KMIN,KMAX IF (JA(K) .EQ. II) GO TO 40 20 CONTINUE 30 NUM = NUM + 1 40 CONTINUE 50 CONTINUE NZSUT = NUM RETURN C----------------------- End of Subroutine CNTNZU ---------------------- END *DECK SPRJS SUBROUTINE SPRJS (NEQ,Y,YH,NYH,EWT,FTEM,SAVF,WK,IWK,F,JAC) EXTERNAL F,JAC INTEGER NEQ, NYH, IWK REAL Y, YH, EWT, FTEM, SAVF, WK DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*), 1 WK(*), IWK(*) INTEGER IOWND, IOWNS, 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP, 1 IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA, 2 LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ, 3 NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU REAL ROWNS, 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND REAL CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH COMMON /SLS001/ ROWNS(209), 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 2 IOWND(6), IOWNS(6), 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU COMMON /SLSS01/ CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH, 1 IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP, 2 IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA, 3 LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ, 4 NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU INTEGER I, IMUL, J, JJ, JOK, JMAX, JMIN, K, KMAX, KMIN, NG REAL CON, DI, FAC, HL0, PIJ, R, R0, RCON, RCONT, 1 SRUR, SVNORM C----------------------------------------------------------------------- C SPRJS is called to compute and process the matrix C P = I - H*EL(1)*J , where J is an approximation to the Jacobian. C J is computed by columns, either by the user-supplied routine JAC C if MITER = 1, or by finite differencing if MITER = 2. C if MITER = 3, a diagonal approximation to J is used. C if MITER = 1 or 2, and if the existing value of the Jacobian C (as contained in P) is considered acceptable, then a new value of C P is reconstructed from the old value. In any case, when MITER C is 1 or 2, the P matrix is subjected to LU decomposition in SCDRV. C P and its LU decomposition are stored (separately) in WK. C C In addition to variables described previously, communication C with SPRJS uses the following: C Y = array containing predicted values on entry. C FTEM = work array of length N (ACOR in SSTODE). C SAVF = array containing f evaluated at predicted y. C WK = real work space for matrices. On output it contains the C inverse diagonal matrix if MITER = 3, and P and its sparse C LU decomposition if MITER is 1 or 2. C Storage of matrix elements starts at WK(3). C WK also contains the following matrix-related data: C WK(1) = SQRT(UROUND), used in numerical Jacobian increments. C WK(2) = H*EL0, saved for later use if MITER = 3. C IWK = integer work space for matrix-related data, assumed to C be equivalenced to WK. In addition, WK(IPRSP) and IWK(IPISP) C are assumed to have identical locations. C EL0 = EL(1) (input). C IERPJ = output error flag (in Common). C = 0 if no error. C = 1 if zero pivot found in SCDRV. C = 2 if a singular matrix arose with MITER = 3. C = -1 if insufficient storage for SCDRV (should not occur here). C = -2 if other error found in SCDRV (should not occur here). C JCUR = output flag showing status of (approximate) Jacobian matrix: C = 1 to indicate that the Jacobian is now current, or C = 0 to indicate that a saved value was used. C This routine also uses other variables in Common. C----------------------------------------------------------------------- HL0 = H*EL0 CON = -HL0 IF (MITER .EQ. 3) GO TO 300 C See whether J should be reevaluated (JOK = 0) or not (JOK = 1). ------ JOK = 1 IF (NST .EQ. 0 .OR. NST .GE. NSLJ+MSBJ) JOK = 0 IF (ICF .EQ. 1 .AND. ABS(RC - 1.0E0) .LT. CCMXJ) JOK = 0 IF (ICF .EQ. 2) JOK = 0 IF (JOK .EQ. 1) GO TO 250 C C MITER = 1 or 2, and the Jacobian is to be reevaluated. --------------- 20 JCUR = 1 NJE = NJE + 1 NSLJ = NST IPLOST = 0 CONMIN = ABS(CON) GO TO (100, 200), MITER C C If MITER = 1, call JAC, multiply by scalar, and add identity. -------- 100 CONTINUE KMIN = IWK(IPIAN) DO 130 J = 1, N KMAX = IWK(IPIAN+J) - 1 DO 110 I = 1,N 110 FTEM(I) = 0.0E0 CALL JAC (NEQ, TN, Y, J, IWK(IPIAN), IWK(IPJAN), FTEM) DO 120 K = KMIN, KMAX I = IWK(IBJAN+K) WK(IBA+K) = FTEM(I)*CON IF (I .EQ. J) WK(IBA+K) = WK(IBA+K) + 1.0E0 120 CONTINUE KMIN = KMAX + 1 130 CONTINUE GO TO 290 C C If MITER = 2, make NGP calls to F to approximate J and P. ------------ 200 CONTINUE FAC = SVNORM(N, SAVF, EWT) R0 = 1000.0E0 * ABS(H) * UROUND * N * FAC IF (R0 .EQ. 0.0E0) R0 = 1.0E0 SRUR = WK(1) JMIN = IWK(IPIGP) DO 240 NG = 1,NGP JMAX = IWK(IPIGP+NG) - 1 DO 210 J = JMIN,JMAX JJ = IWK(IBJGP+J) R = MAX(SRUR*ABS(Y(JJ)),R0/EWT(JJ)) 210 Y(JJ) = Y(JJ) + R CALL F (NEQ, TN, Y, FTEM) DO 230 J = JMIN,JMAX JJ = IWK(IBJGP+J) Y(JJ) = YH(JJ,1) R = MAX(SRUR*ABS(Y(JJ)),R0/EWT(JJ)) FAC = -HL0/R KMIN =IWK(IBIAN+JJ) KMAX =IWK(IBIAN+JJ+1) - 1 DO 220 K = KMIN,KMAX I = IWK(IBJAN+K) WK(IBA+K) = (FTEM(I) - SAVF(I))*FAC IF (I .EQ. JJ) WK(IBA+K) = WK(IBA+K) + 1.0E0 220 CONTINUE 230 CONTINUE JMIN = JMAX + 1 240 CONTINUE NFE = NFE + NGP GO TO 290 C C If JOK = 1, reconstruct new P from old P. ---------------------------- 250 JCUR = 0 RCON = CON/CON0 RCONT = ABS(CON)/CONMIN IF (RCONT .GT. RBIG .AND. IPLOST .EQ. 1) GO TO 20 KMIN = IWK(IPIAN) DO 275 J = 1,N KMAX = IWK(IPIAN+J) - 1 DO 270 K = KMIN,KMAX I = IWK(IBJAN+K) PIJ = WK(IBA+K) IF (I .NE. J) GO TO 260 PIJ = PIJ - 1.0E0 IF (ABS(PIJ) .GE. PSMALL) GO TO 260 IPLOST = 1 CONMIN = MIN(ABS(CON0),CONMIN) 260 PIJ = PIJ*RCON IF (I .EQ. J) PIJ = PIJ + 1.0E0 WK(IBA+K) = PIJ 270 CONTINUE KMIN = KMAX + 1 275 CONTINUE C C Do numerical factorization of P matrix. ------------------------------ 290 NLU = NLU + 1 CON0 = CON IERPJ = 0 DO 295 I = 1,N 295 FTEM(I) = 0.0E0 CALL SCDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN), 1 WK(IPA),FTEM,FTEM,NSP,IWK(IPISP),WK(IPRSP),IESP,2,IYS) IF (IYS .EQ. 0) RETURN IMUL = (IYS - 1)/N IERPJ = -2 IF (IMUL .EQ. 8) IERPJ = 1 IF (IMUL .EQ. 10) IERPJ = -1 RETURN C C If MITER = 3, construct a diagonal approximation to J and P. --------- 300 CONTINUE JCUR = 1 NJE = NJE + 1 WK(2) = HL0 IERPJ = 0 R = EL0*0.1E0 DO 310 I = 1,N 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) CALL F (NEQ, TN, Y, WK(3)) NFE = NFE + 1 DO 320 I = 1,N R0 = H*SAVF(I) - YH(I,2) DI = 0.1E0*R0 - H*(WK(I+2) - SAVF(I)) WK(I+2) = 1.0E0 IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320 IF (ABS(DI) .EQ. 0.0E0) GO TO 330 WK(I+2) = 0.1E0*R0/DI 320 CONTINUE RETURN 330 IERPJ = 2 RETURN C----------------------- End of Subroutine SPRJS ----------------------- END *DECK SSOLSS SUBROUTINE SSOLSS (WK, IWK, X, TEM) INTEGER IWK REAL WK, X, TEM DIMENSION WK(*), IWK(*), X(*), TEM(*) INTEGER IOWND, IOWNS, 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP, 1 IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA, 2 LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ, 3 NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU REAL ROWNS, 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND REAL RLSS COMMON /SLS001/ ROWNS(209), 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 2 IOWND(6), IOWNS(6), 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU COMMON /SLSS01/ RLSS(6), 1 IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP, 2 IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA, 3 LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ, 4 NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU INTEGER I REAL DI, HL0, PHL0, R C----------------------------------------------------------------------- C This routine manages the solution of the linear system arising from C a chord iteration. It is called if MITER .ne. 0. C If MITER is 1 or 2, it calls SCDRV to accomplish this. C If MITER = 3 it updates the coefficient H*EL0 in the diagonal C matrix, and then computes the solution. C communication with SSOLSS uses the following variables: C WK = real work space containing the inverse diagonal matrix if C MITER = 3 and the LU decomposition of the matrix otherwise. C Storage of matrix elements starts at WK(3). C WK also contains the following matrix-related data: C WK(1) = SQRT(UROUND) (not used here), C WK(2) = HL0, the previous value of H*EL0, used if MITER = 3. C IWK = integer work space for matrix-related data, assumed to C be equivalenced to WK. In addition, WK(IPRSP) and IWK(IPISP) C are assumed to have identical locations. C X = the right-hand side vector on input, and the solution vector C on output, of length N. C TEM = vector of work space of length N, not used in this version. C IERSL = output flag (in Common). C IERSL = 0 if no trouble occurred. C IERSL = -1 if SCDRV returned an error flag (MITER = 1 or 2). C This should never occur and is considered fatal. C IERSL = 1 if a singular matrix arose with MITER = 3. C This routine also uses other variables in Common. C----------------------------------------------------------------------- IERSL = 0 GO TO (100, 100, 300), MITER 100 CALL SCDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN), 1 WK(IPA),X,X,NSP,IWK(IPISP),WK(IPRSP),IESP,4,IERSL) IF (IERSL .NE. 0) IERSL = -1 RETURN C 300 PHL0 = WK(2) HL0 = H*EL0 WK(2) = HL0 IF (HL0 .EQ. PHL0) GO TO 330 R = HL0/PHL0 DO 320 I = 1,N DI = 1.0E0 - R*(1.0E0 - 1.0E0/WK(I+2)) IF (ABS(DI) .EQ. 0.0E0) GO TO 390 320 WK(I+2) = 1.0E0/DI 330 DO 340 I = 1,N 340 X(I) = WK(I+2)*X(I) RETURN 390 IERSL = 1 RETURN C C----------------------- End of Subroutine SSOLSS ---------------------- END *DECK SSRCMS SUBROUTINE SSRCMS (RSAV, ISAV, JOB) C----------------------------------------------------------------------- C This routine saves or restores (depending on JOB) the contents of C the Common blocks SLS001, SLSS01, which are used C internally by one or more ODEPACK solvers. C C RSAV = real array of length 224 or more. C ISAV = integer array of length 71 or more. C JOB = flag indicating to save or restore the Common blocks: C JOB = 1 if Common is to be saved (written to RSAV/ISAV) C JOB = 2 if Common is to be restored (read from RSAV/ISAV) C A call with JOB = 2 presumes a prior call with JOB = 1. C----------------------------------------------------------------------- INTEGER ISAV, JOB INTEGER ILS, ILSS INTEGER I, LENILS, LENISS, LENRLS, LENRSS REAL RSAV, RLS, RLSS DIMENSION RSAV(*), ISAV(*) SAVE LENRLS, LENILS, LENRSS, LENISS COMMON /SLS001/ RLS(218), ILS(37) COMMON /SLSS01/ RLSS(6), ILSS(34) DATA LENRLS/218/, LENILS/37/, LENRSS/6/, LENISS/34/ C IF (JOB .EQ. 2) GO TO 100 DO 10 I = 1,LENRLS 10 RSAV(I) = RLS(I) DO 15 I = 1,LENRSS 15 RSAV(LENRLS+I) = RLSS(I) C DO 20 I = 1,LENILS 20 ISAV(I) = ILS(I) DO 25 I = 1,LENISS 25 ISAV(LENILS+I) = ILSS(I) C RETURN C 100 CONTINUE DO 110 I = 1,LENRLS 110 RLS(I) = RSAV(I) DO 115 I = 1,LENRSS 115 RLSS(I) = RSAV(LENRLS+I) C DO 120 I = 1,LENILS 120 ILS(I) = ISAV(I) DO 125 I = 1,LENISS 125 ILSS(I) = ISAV(LENILS+I) C RETURN C----------------------- End of Subroutine SSRCMS ---------------------- END *DECK SODRV subroutine sodrv * (n, ia,ja,a, p,ip, nsp,isp, path, flag) c 5/2/83 c*********************************************************************** c sodrv -- driver for sparse matrix reordering routines c*********************************************************************** c c description c c sodrv finds a minimum degree ordering of the rows and columns c of a matrix m stored in (ia,ja,a) format (see below). for the c reordered matrix, the work and storage required to perform c gaussian elimination is (usually) significantly less. c c note.. sodrv and its subordinate routines have been modified to c compute orderings for general matrices, not necessarily having any c symmetry. the miminum degree ordering is computed for the c structure of the symmetric matrix m + m-transpose. c modifications to the original sodrv module have been made in c the coding in subroutine mdi, and in the initial comments in c subroutines sodrv and md. c c if only the nonzero entries in the upper triangle of m are being c stored, then sodrv symmetrically reorders (ia,ja,a), (optionally) c with the diagonal entries placed first in each row. this is to c ensure that if m(i,j) will be in the upper triangle of m with c respect to the new ordering, then m(i,j) is stored in row i (and c thus m(j,i) is not stored), whereas if m(i,j) will be in the c strict lower triangle of m, then m(j,i) is stored in row j (and c thus m(i,j) is not stored). c c c storage of sparse matrices c c the nonzero entries of the matrix m are stored row-by-row in the c array a. to identify the individual nonzero entries in each row, c we need to know in which column each entry lies. these column c indices are stored in the array ja. i.e., if a(k) = m(i,j), then c ja(k) = j. to identify the individual rows, we need to know where c each row starts. these row pointers are stored in the array ia. c i.e., if m(i,j) is the first nonzero entry (stored) in the i-th row c and a(k) = m(i,j), then ia(i) = k. moreover, ia(n+1) points to c the first location following the last element in the last row. c thus, the number of entries in the i-th row is ia(i+1) - ia(i), c the nonzero entries in the i-th row are stored consecutively in c c a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1), c c and the corresponding column indices are stored consecutively in c c ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1). c c when the coefficient matrix is symmetric, only the nonzero entries c in the upper triangle need be stored. for example, the matrix c c ( 1 0 2 3 0 ) c ( 0 4 0 0 0 ) c m = ( 2 0 5 6 0 ) c ( 3 0 6 7 8 ) c ( 0 0 0 8 9 ) c c could be stored as c c - 1 2 3 4 5 6 7 8 9 10 11 12 13 c ---+-------------------------------------- c ia - 1 4 5 8 12 14 c ja - 1 3 4 2 1 3 4 1 3 4 5 4 5 c a - 1 2 3 4 2 5 6 3 6 7 8 8 9 c c or (symmetrically) as c c - 1 2 3 4 5 6 7 8 9 c ---+-------------------------- c ia - 1 4 5 7 9 10 c ja - 1 3 4 2 3 4 4 5 5 c a - 1 2 3 4 5 6 7 8 9 . c c c parameters c c n - order of the matrix c c ia - integer one-dimensional array containing pointers to delimit c rows in ja and a. dimension = n+1 c c ja - integer one-dimensional array containing the column indices c corresponding to the elements of a. dimension = number of c nonzero entries in (the upper triangle of) m c c a - real one-dimensional array containing the nonzero entries in c (the upper triangle of) m, stored by rows. dimension = c number of nonzero entries in (the upper triangle of) m c c p - integer one-dimensional array used to return the permutation c of the rows and columns of m corresponding to the minimum c degree ordering. dimension = n c c ip - integer one-dimensional array used to return the inverse of c the permutation returned in p. dimension = n c c nsp - declared dimension of the one-dimensional array isp. nsp c must be at least 3n+4k, where k is the number of nonzeroes c in the strict upper triangle of m c c isp - integer one-dimensional array used for working storage. c dimension = nsp c c path - integer path specification. values and their meanings are - c 1 find minimum degree ordering only c 2 find minimum degree ordering and reorder symmetrically c stored matrix (used when only the nonzero entries in c the upper triangle of m are being stored) c 3 reorder symmetrically stored matrix as specified by c input permutation (used when an ordering has already c been determined and only the nonzero entries in the c upper triangle of m are being stored) c 4 same as 2 but put diagonal entries at start of each row c 5 same as 3 but put diagonal entries at start of each row c c flag - integer error flag. values and their meanings are - c 0 no errors detected c 9n+k insufficient storage in md c 10n+1 insufficient storage in sodrv c 11n+1 illegal path specification c c c conversion from real to double precision c c change the real declarations in sodrv and ssro to double precision c declarations. c c----------------------------------------------------------------------- c integer ia(*), ja(*), p(*), ip(*), isp(*), path, flag, * v, l, head, tmp, q real a(*) c... double precision a(*) logical dflag c c----initialize error flag and validate path specification flag = 0 if (path.lt.1 .or. 5.lt.path) go to 111 c c----allocate storage and find minimum degree ordering if ((path-1) * (path-2) * (path-4) .ne. 0) go to 1 max = (nsp-n)/2 v = 1 l = v + max head = l + max next = head + n if (max.lt.n) go to 110 c call md * (n, ia,ja, max,isp(v),isp(l), isp(head),p,ip, isp(v), flag) if (flag.ne.0) go to 100 c c----allocate storage and symmetrically reorder matrix 1 if ((path-2) * (path-3) * (path-4) * (path-5) .ne. 0) go to 2 tmp = (nsp+1) - n q = tmp - (ia(n+1)-1) if (q.lt.1) go to 110 c dflag = path.eq.4 .or. path.eq.5 call ssro * (n, ip, ia, ja, a, isp(tmp), isp(q), dflag) c 2 return c c ** error -- error detected in md 100 return c ** error -- insufficient storage 110 flag = 10*n + 1 return c ** error -- illegal path specified 111 flag = 11*n + 1 return end subroutine md * (n, ia,ja, max, v,l, head,last,next, mark, flag) c*********************************************************************** c md -- minimum degree algorithm (based on element model) c*********************************************************************** c c description c c md finds a minimum degree ordering of the rows and columns of a c general sparse matrix m stored in (ia,ja,a) format. c when the structure of m is nonsymmetric, the ordering is that c obtained for the symmetric matrix m + m-transpose. c c c additional parameters c c max - declared dimension of the one-dimensional arrays v and l. c max must be at least n+2k, where k is the number of c nonzeroes in the strict upper triangle of m + m-transpose c c v - integer one-dimensional work array. dimension = max c c l - integer one-dimensional work array. dimension = max c c head - integer one-dimensional work array. dimension = n c c last - integer one-dimensional array used to return the permutation c of the rows and columns of m corresponding to the minimum c degree ordering. dimension = n c c next - integer one-dimensional array used to return the inverse of c the permutation returned in last. dimension = n c c mark - integer one-dimensional work array (may be the same as v). c dimension = n c c flag - integer error flag. values and their meanings are - c 0 no errors detected c 9n+k insufficient storage in md c c c definitions of internal parameters c c ---------+--------------------------------------------------------- c v(s) - value field of list entry c ---------+--------------------------------------------------------- c l(s) - link field of list entry (0 =) end of list) c ---------+--------------------------------------------------------- c l(vi) - pointer to element list of uneliminated vertex vi c ---------+--------------------------------------------------------- c l(ej) - pointer to boundary list of active element ej c ---------+--------------------------------------------------------- c head(d) - vj =) vj head of d-list d c - 0 =) no vertex in d-list d c c c - vi uneliminated vertex c - vi in ek - vi not in ek c ---------+-----------------------------+--------------------------- c next(vi) - undefined but nonnegative - vj =) vj next in d-list c - - 0 =) vi tail of d-list c ---------+-----------------------------+--------------------------- c last(vi) - (not set until mdp) - -d =) vi head of d-list d c --vk =) compute degree - vj =) vj last in d-list c - ej =) vi prototype of ej - 0 =) vi not in any d-list c - 0 =) do not compute degree - c ---------+-----------------------------+--------------------------- c mark(vi) - mark(vk) - nonneg. tag .lt. mark(vk) c c c - vi eliminated vertex c - ei active element - otherwise c ---------+-----------------------------+--------------------------- c next(vi) - -j =) vi was j-th vertex - -j =) vi was j-th vertex c - to be eliminated - to be eliminated c ---------+-----------------------------+--------------------------- c last(vi) - m =) size of ei = m - undefined c ---------+-----------------------------+--------------------------- c mark(vi) - -m =) overlap count of ei - undefined c - with ek = m - c - otherwise nonnegative tag - c - .lt. mark(vk) - c c----------------------------------------------------------------------- c integer ia(*), ja(*), v(*), l(*), head(*), last(*), next(*), * mark(*), flag, tag, dmin, vk,ek, tail equivalence (vk,ek) c c----initialization tag = 0 call mdi * (n, ia,ja, max,v,l, head,last,next, mark,tag, flag) if (flag.ne.0) return c k = 0 dmin = 1 c c----while k .lt. n do 1 if (k.ge.n) go to 4 c c------search for vertex of minimum degree 2 if (head(dmin).gt.0) go to 3 dmin = dmin + 1 go to 2 c c------remove vertex vk of minimum degree from degree list 3 vk = head(dmin) head(dmin) = next(vk) if (head(dmin).gt.0) last(head(dmin)) = -dmin c c------number vertex vk, adjust tag, and tag vk k = k+1 next(vk) = -k last(ek) = dmin - 1 tag = tag + last(ek) mark(vk) = tag c c------form element ek from uneliminated neighbors of vk call mdm * (vk,tail, v,l, last,next, mark) c c------purge inactive elements and do mass elimination call mdp * (k,ek,tail, v,l, head,last,next, mark) c c------update degrees of uneliminated vertices in ek call mdu * (ek,dmin, v,l, head,last,next, mark) c go to 1 c c----generate inverse permutation from permutation 4 do 5 k=1,n next(k) = -next(k) 5 last(next(k)) = k c return end subroutine mdi * (n, ia,ja, max,v,l, head,last,next, mark,tag, flag) c*********************************************************************** c mdi -- initialization c*********************************************************************** integer ia(*), ja(*), v(*), l(*), head(*), last(*), next(*), * mark(*), tag, flag, sfs, vi,dvi, vj c c----initialize degrees, element lists, and degree lists do 1 vi=1,n mark(vi) = 1 l(vi) = 0 1 head(vi) = 0 sfs = n+1 c c----create nonzero structure c----for each nonzero entry a(vi,vj) do 6 vi=1,n jmin = ia(vi) jmax = ia(vi+1) - 1 if (jmin.gt.jmax) go to 6 do 5 j=jmin,jmax vj = ja(j) if (vj-vi) 2, 5, 4 c c------if a(vi,vj) is in strict lower triangle c------check for previous occurrence of a(vj,vi) 2 lvk = vi kmax = mark(vi) - 1 if (kmax .eq. 0) go to 4 do 3 k=1,kmax lvk = l(lvk) if (v(lvk).eq.vj) go to 5 3 continue c----for unentered entries a(vi,vj) 4 if (sfs.ge.max) go to 101 c c------enter vj in element list for vi mark(vi) = mark(vi) + 1 v(sfs) = vj l(sfs) = l(vi) l(vi) = sfs sfs = sfs+1 c c------enter vi in element list for vj mark(vj) = mark(vj) + 1 v(sfs) = vi l(sfs) = l(vj) l(vj) = sfs sfs = sfs+1 5 continue 6 continue c c----create degree lists and initialize mark vector do 7 vi=1,n dvi = mark(vi) next(vi) = head(dvi) head(dvi) = vi last(vi) = -dvi nextvi = next(vi) if (nextvi.gt.0) last(nextvi) = vi 7 mark(vi) = tag c return c c ** error- insufficient storage 101 flag = 9*n + vi return end subroutine mdm * (vk,tail, v,l, last,next, mark) c*********************************************************************** c mdm -- form element from uneliminated neighbors of vk c*********************************************************************** integer vk, tail, v(*), l(*), last(*), next(*), mark(*), * tag, s,ls,vs,es, b,lb,vb, blp,blpmax equivalence (vs, es) c c----initialize tag and list of uneliminated neighbors tag = mark(vk) tail = vk c c----for each vertex/element vs/es in element list of vk ls = l(vk) 1 s = ls if (s.eq.0) go to 5 ls = l(s) vs = v(s) if (next(vs).lt.0) go to 2 c c------if vs is uneliminated vertex, then tag and append to list of c------uneliminated neighbors mark(vs) = tag l(tail) = s tail = s go to 4 c c------if es is active element, then ... c--------for each vertex vb in boundary list of element es 2 lb = l(es) blpmax = last(es) do 3 blp=1,blpmax b = lb lb = l(b) vb = v(b) c c----------if vb is untagged vertex, then tag and append to list of c----------uneliminated neighbors if (mark(vb).ge.tag) go to 3 mark(vb) = tag l(tail) = b tail = b 3 continue c c--------mark es inactive mark(es) = tag c 4 go to 1 c c----terminate list of uneliminated neighbors 5 l(tail) = 0 c return end subroutine mdp * (k,ek,tail, v,l, head,last,next, mark) c*********************************************************************** c mdp -- purge inactive elements and do mass elimination c*********************************************************************** integer ek, tail, v(*), l(*), head(*), last(*), next(*), * mark(*), tag, free, li,vi,lvi,evi, s,ls,es, ilp,ilpmax c c----initialize tag tag = mark(ek) c c----for each vertex vi in ek li = ek ilpmax = last(ek) if (ilpmax.le.0) go to 12 do 11 ilp=1,ilpmax i = li li = l(i) vi = v(li) c c------remove vi from degree list if (last(vi).eq.0) go to 3 if (last(vi).gt.0) go to 1 head(-last(vi)) = next(vi) go to 2 1 next(last(vi)) = next(vi) 2 if (next(vi).gt.0) last(next(vi)) = last(vi) c c------remove inactive items from element list of vi 3 ls = vi 4 s = ls ls = l(s) if (ls.eq.0) go to 6 es = v(ls) if (mark(es).lt.tag) go to 5 free = ls l(s) = l(ls) ls = s 5 go to 4 c c------if vi is interior vertex, then remove from list and eliminate 6 lvi = l(vi) if (lvi.ne.0) go to 7 l(i) = l(li) li = i c k = k+1 next(vi) = -k last(ek) = last(ek) - 1 go to 11 c c------else ... c--------classify vertex vi 7 if (l(lvi).ne.0) go to 9 evi = v(lvi) if (next(evi).ge.0) go to 9 if (mark(evi).lt.0) go to 8 c c----------if vi is prototype vertex, then mark as such, initialize c----------overlap count for corresponding element, and move vi to end c----------of boundary list last(vi) = evi mark(evi) = -1 l(tail) = li tail = li l(i) = l(li) li = i go to 10 c c----------else if vi is duplicate vertex, then mark as such and adjust c----------overlap count for corresponding element 8 last(vi) = 0 mark(evi) = mark(evi) - 1 go to 10 c c----------else mark vi to compute degree 9 last(vi) = -ek c c--------insert ek in element list of vi 10 v(free) = ek l(free) = l(vi) l(vi) = free 11 continue c c----terminate boundary list 12 l(tail) = 0 c return end subroutine mdu * (ek,dmin, v,l, head,last,next, mark) c*********************************************************************** c mdu -- update degrees of uneliminated vertices in ek c*********************************************************************** integer ek, dmin, v(*), l(*), head(*), last(*), next(*), * mark(*), tag, vi,evi,dvi, s,vs,es, b,vb, ilp,ilpmax, * blp,blpmax equivalence (vs, es) c c----initialize tag tag = mark(ek) - last(ek) c c----for each vertex vi in ek i = ek ilpmax = last(ek) if (ilpmax.le.0) go to 11 do 10 ilp=1,ilpmax i = l(i) vi = v(i) if (last(vi)) 1, 10, 8 c c------if vi neither prototype nor duplicate vertex, then merge elements c------to compute degree 1 tag = tag + 1 dvi = last(ek) c c--------for each vertex/element vs/es in element list of vi s = l(vi) 2 s = l(s) if (s.eq.0) go to 9 vs = v(s) if (next(vs).lt.0) go to 3 c c----------if vs is uneliminated vertex, then tag and adjust degree mark(vs) = tag dvi = dvi + 1 go to 5 c c----------if es is active element, then expand c------------check for outmatched vertex 3 if (mark(es).lt.0) go to 6 c c------------for each vertex vb in es b = es blpmax = last(es) do 4 blp=1,blpmax b = l(b) vb = v(b) c c--------------if vb is untagged, then tag and adjust degree if (mark(vb).ge.tag) go to 4 mark(vb) = tag dvi = dvi + 1 4 continue c 5 go to 2 c c------else if vi is outmatched vertex, then adjust overlaps but do not c------compute degree 6 last(vi) = 0 mark(es) = mark(es) - 1 7 s = l(s) if (s.eq.0) go to 10 es = v(s) if (mark(es).lt.0) mark(es) = mark(es) - 1 go to 7 c c------else if vi is prototype vertex, then calculate degree by c------inclusion/exclusion and reset overlap count 8 evi = last(vi) dvi = last(ek) + last(evi) + mark(evi) mark(evi) = 0 c c------insert vi in appropriate degree list 9 next(vi) = head(dvi) head(dvi) = vi last(vi) = -dvi if (next(vi).gt.0) last(next(vi)) = vi if (dvi.lt.dmin) dmin = dvi c 10 continue c 11 return end subroutine ssro * (n, ip, ia,ja,a, q, r, dflag) c*********************************************************************** c ssro -- symmetric reordering of sparse symmetric matrix c*********************************************************************** c c description c c the nonzero entries of the matrix m are assumed to be stored c symmetrically in (ia,ja,a) format (i.e., not both m(i,j) and m(j,i) c are stored if i ne j). c c ssro does not rearrange the order of the rows, but does move c nonzeroes from one row to another to ensure that if m(i,j) will be c in the upper triangle of m with respect to the new ordering, then c m(i,j) is stored in row i (and thus m(j,i) is not stored), whereas c if m(i,j) will be in the strict lower triangle of m, then m(j,i) is c stored in row j (and thus m(i,j) is not stored). c c c additional parameters c c q - integer one-dimensional work array. dimension = n c c r - integer one-dimensional work array. dimension = number of c nonzero entries in the upper triangle of m c c dflag - logical variable. if dflag = .true., then store nonzero c diagonal elements at the beginning of the row c c----------------------------------------------------------------------- c integer ip(*), ia(*), ja(*), q(*), r(*) real a(*), ak c... double precision a(*), ak logical dflag c c c--phase 1 -- find row in which to store each nonzero c----initialize count of nonzeroes to be stored in each row do 1 i=1,n 1 q(i) = 0 c c----for each nonzero element a(j) do 3 i=1,n jmin = ia(i) jmax = ia(i+1) - 1 if (jmin.gt.jmax) go to 3 do 2 j=jmin,jmax c c--------find row (=r(j)) and column (=ja(j)) in which to store a(j) ... k = ja(j) if (ip(k).lt.ip(i)) ja(j) = i if (ip(k).ge.ip(i)) k = i r(j) = k c c--------... and increment count of nonzeroes (=q(r(j)) in that row 2 q(k) = q(k) + 1 3 continue c c c--phase 2 -- find new ia and permutation to apply to (ja,a) c----determine pointers to delimit rows in permuted (ja,a) do 4 i=1,n ia(i+1) = ia(i) + q(i) 4 q(i) = ia(i+1) c c----determine where each (ja(j),a(j)) is stored in permuted (ja,a) c----for each nonzero element (in reverse order) ilast = 0 jmin = ia(1) jmax = ia(n+1) - 1 j = jmax do 6 jdummy=jmin,jmax i = r(j) if (.not.dflag .or. ja(j).ne.i .or. i.eq.ilast) go to 5 c c------if dflag, then put diagonal nonzero at beginning of row r(j) = ia(i) ilast = i go to 6 c c------put (off-diagonal) nonzero in last unused location in row 5 q(i) = q(i) - 1 r(j) = q(i) c 6 j = j-1 c c c--phase 3 -- permute (ja,a) to upper triangular form (wrt new ordering) do 8 j=jmin,jmax 7 if (r(j).eq.j) go to 8 k = r(j) r(j) = r(k) r(k) = k jak = ja(k) ja(k) = ja(j) ja(j) = jak ak = a(k) a(k) = a(j) a(j) = ak go to 7 8 continue c return end *DECK SCDRV subroutine scdrv * (n, r,c,ic, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag) c*** subroutine scdrv c*** driver for subroutines for solving sparse nonsymmetric systems of c linear equations (compressed pointer storage) c c c parameters c class abbreviations are-- c n - integer variable c f - real variable c v - supplies a value to the driver c r - returns a result from the driver c i - used internally by the driver c a - array c c class - parameter c ------+---------- c - c the nonzero entries of the coefficient matrix m are stored c row-by-row in the array a. to identify the individual nonzero c entries in each row, we need to know in which column each entry c lies. the column indices which correspond to the nonzero entries c of m are stored in the array ja. i.e., if a(k) = m(i,j), then c ja(k) = j. in addition, we need to know where each row starts and c how long it is. the index positions in ja and a where the rows of c m begin are stored in the array ia. i.e., if m(i,j) is the first c nonzero entry (stored) in the i-th row and a(k) = m(i,j), then c ia(i) = k. moreover, the index in ja and a of the first location c following the last element in the last row is stored in ia(n+1). c thus, the number of entries in the i-th row is given by c ia(i+1) - ia(i), the nonzero entries of the i-th row are stored c consecutively in c a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1), c and the corresponding column indices are stored consecutively in c ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1). c for example, the 5 by 5 matrix c ( 1. 0. 2. 0. 0.) c ( 0. 3. 0. 0. 0.) c m = ( 0. 4. 5. 6. 0.) c ( 0. 0. 0. 7. 0.) c ( 0. 0. 0. 8. 9.) c would be stored as c - 1 2 3 4 5 6 7 8 9 c ---+-------------------------- c ia - 1 3 4 7 8 10 c ja - 1 3 2 2 3 4 4 4 5 c a - 1. 2. 3. 4. 5. 6. 7. 8. 9. . c c nv - n - number of variables/equations. c fva - a - nonzero entries of the coefficient matrix m, stored c - by rows. c - size = number of nonzero entries in m. c nva - ia - pointers to delimit the rows in a. c - size = n+1. c nva - ja - column numbers corresponding to the elements of a. c - size = size of a. c fva - b - right-hand side b. b and z can the same array. c - size = n. c fra - z - solution x. b and z can be the same array. c - size = n. c c the rows and columns of the original matrix m can be c reordered (e.g., to reduce fillin or ensure numerical stability) c before calling the driver. if no reordering is done, then set c r(i) = c(i) = ic(i) = i for i=1,...,n. the solution z is returned c in the original order. c if the columns have been reordered (i.e., c(i).ne.i for some c i), then the driver will call a subroutine (snroc) which rearranges c each row of ja and a, leaving the rows in the original order, but c placing the elements of each row in increasing order with respect c to the new ordering. if path.ne.1, then snroc is assumed to have c been called already. c c nva - r - ordering of the rows of m. c - size = n. c nva - c - ordering of the columns of m. c - size = n. c nva - ic - inverse of the ordering of the columns of m. i.e., c - ic(c(i)) = i for i=1,...,n. c - size = n. c c the solution of the system of linear equations is divided into c three stages -- c snsfc -- the matrix m is processed symbolically to determine where c fillin will occur during the numeric factorization. c snnfc -- the matrix m is factored numerically into the product ldu c of a unit lower triangular matrix l, a diagonal matrix c d, and a unit upper triangular matrix u, and the system c mx = b is solved. c snnsc -- the linear system mx = b is solved using the ldu c or factorization from snnfc. c snntc -- the transposed linear system mt x = b is solved using c the ldu factorization from nnf. c for several systems whose coefficient matrices have the same c nonzero structure, snsfc need be done only once (for the first c system). then snnfc is done once for each additional system. for c several systems with the same coefficient matrix, snsfc and snnfc c need be done only once (for the first system). then snnsc or snntc c is done once for each additional right-hand side. c c nv - path - path specification. values and their meanings are -- c - 1 perform snroc, snsfc, and snnfc. c - 2 perform snnfc only (snsfc is assumed to have been c - done in a manner compatible with the storage c - allocation used in the driver). c - 3 perform snnsc only (snsfc and snnfc are assumed c - to have been done in a manner compatible with c - the storage allocation used in the driver). c - 4 perform snntc only (snsfc and snnfc are assumed c - to have been done in a manner compatible with c - the storage allocation used in the driver). c - 5 perform snroc and snsfc. c c various errors are detected by the driver and the individual c subroutines. c c nr - flag - error flag. values and their meanings are -- c - 0 no errors detected c - n+k null row in a -- row = k c - 2n+k duplicate entry in a -- row = k c - 3n+k insufficient storage in snsfc -- row = k c - 4n+1 insufficient storage in snnfc c - 5n+k null pivot -- row = k c - 6n+k insufficient storage in snsfc -- row = k c - 7n+1 insufficient storage in snnfc c - 8n+k zero pivot -- row = k c - 10n+1 insufficient storage in scdrv c - 11n+1 illegal path specification c c working storage is needed for the factored form of the matrix c m plus various temporary vectors. the arrays isp and rsp should be c equivalenced. integer storage is allocated from the beginning of c isp and real storage from the end of rsp. c c nv - nsp - declared dimension of rsp. nsp generally must c - be larger than 8n+2 + 2k (where k = (number of c - nonzero entries in m)). c nvira - isp - integer working storage divided up into various arrays c - needed by the subroutines. isp and rsp should be c - equivalenced. c - size = lratio*nsp. c fvira - rsp - real working storage divided up into various arrays c - needed by the subroutines. isp and rsp should be c - equivalenced. c - size = nsp. c nr - esp - if sufficient storage was available to perform the c - symbolic factorization (snsfc), then esp is set to c - the amount of excess storage provided (negative if c - insufficient storage was available to perform the c - numeric factorization (snnfc)). c c c conversion to double precision c c to convert these routines for double precision arrays.. c (1) use the double precision declarations in place of the real c declarations in each subprogram, as given in comment cards. c (2) change the data-loaded value of the integer lratio c in subroutine scdrv, as indicated below. c (3) change e0 to d0 in the constants in statement number 10 c in subroutine snnfc and the line following that. c integer r(*), c(*), ic(*), ia(*), ja(*), isp(*), esp, path, * flag, d, u, q, row, tmp, ar, umax real a(*), b(*), z(*), rsp(*) c double precision a(*), b(*), z(*), rsp(*) c c set lratio equal to the ratio between the length of floating point c and integer array data. e. g., lratio = 1 for (real, integer), c lratio = 2 for (double precision, integer) c data lratio/1/ c if (path.lt.1 .or. 5.lt.path) go to 111 c******initialize and divide up temporary storage ******************* il = 1 ijl = il + (n+1) iu = ijl + n iju = iu + (n+1) irl = iju + n jrl = irl + n jl = jrl + n c c ****** reorder a if necessary, call snsfc if flag is set ********** if ((path-1) * (path-5) .ne. 0) go to 5 max = (lratio*nsp + 1 - jl) - (n+1) - 5*n jlmax = max/2 q = jl + jlmax ira = q + (n+1) jra = ira + n irac = jra + n iru = irac + n jru = iru + n jutmp = jru + n jumax = lratio*nsp + 1 - jutmp esp = max/lratio if (jlmax.le.0 .or. jumax.le.0) go to 110 c do 1 i=1,n if (c(i).ne.i) go to 2 1 continue go to 3 2 ar = nsp + 1 - n call snroc * (n, ic, ia,ja,a, isp(il), rsp(ar), isp(iu), flag) if (flag.ne.0) go to 100 c 3 call snsfc * (n, r, ic, ia,ja, * jlmax, isp(il), isp(jl), isp(ijl), * jumax, isp(iu), isp(jutmp), isp(iju), * isp(q), isp(ira), isp(jra), isp(irac), * isp(irl), isp(jrl), isp(iru), isp(jru), flag) if(flag .ne. 0) go to 100 c ****** move ju next to jl ***************************************** jlmax = isp(ijl+n-1) ju = jl + jlmax jumax = isp(iju+n-1) if (jumax.le.0) go to 5 do 4 j=1,jumax 4 isp(ju+j-1) = isp(jutmp+j-1) c c ****** call remaining subroutines ********************************* 5 jlmax = isp(ijl+n-1) ju = jl + jlmax jumax = isp(iju+n-1) l = (ju + jumax - 2 + lratio) / lratio + 1 lmax = isp(il+n) - 1 d = l + lmax u = d + n row = nsp + 1 - n tmp = row - n umax = tmp - u esp = umax - (isp(iu+n) - 1) c if ((path-1) * (path-2) .ne. 0) go to 6 if (umax.lt.0) go to 110 call snnfc * (n, r, c, ic, ia, ja, a, z, b, * lmax, isp(il), isp(jl), isp(ijl), rsp(l), rsp(d), * umax, isp(iu), isp(ju), isp(iju), rsp(u), * rsp(row), rsp(tmp), isp(irl), isp(jrl), flag) if(flag .ne. 0) go to 100 c 6 if ((path-3) .ne. 0) go to 7 call snnsc * (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l), * rsp(d), isp(iu), isp(ju), isp(iju), rsp(u), * z, b, rsp(tmp)) c 7 if ((path-4) .ne. 0) go to 8 call snntc * (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l), * rsp(d), isp(iu), isp(ju), isp(iju), rsp(u), * z, b, rsp(tmp)) 8 return c c ** error.. error detected in snroc, snsfc, snnfc, or snnsc 100 return c ** error.. insufficient storage 110 flag = 10*n + 1 return c ** error.. illegal path specification 111 flag = 11*n + 1 return end subroutine snroc (n, ic, ia, ja, a, jar, ar, p, flag) c c ---------------------------------------------------------------- c c yale sparse matrix package - nonsymmetric codes c solving the system of equations mx = b c c i. calling sequences c the coefficient matrix can be processed by an ordering routine c (e.g., to reduce fillin or ensure numerical stability) before using c the remaining subroutines. if no reordering is done, then set c r(i) = c(i) = ic(i) = i for i=1,...,n. if an ordering subroutine c is used, then snroc should be used to reorder the coefficient c matrix. the calling sequence is -- c ( (matrix ordering)) c (snroc (matrix reordering)) c snsfc (symbolic factorization to determine where fillin will c occur during numeric factorization) c snnfc (numeric factorization into product ldu of unit lower c triangular matrix l, diagonal matrix d, and unit c upper triangular matrix u, and solution of linear c system) c snnsc (solution of linear system for additional right-hand c side using ldu factorization from snnfc) c (if only one system of equations is to be solved, then the c subroutine trk should be used.) c c ii. storage of sparse matrices c the nonzero entries of the coefficient matrix m are stored c row-by-row in the array a. to identify the individual nonzero c entries in each row, we need to know in which column each entry c lies. the column indices which correspond to the nonzero entries c of m are stored in the array ja. i.e., if a(k) = m(i,j), then c ja(k) = j. in addition, we need to know where each row starts and c how long it is. the index positions in ja and a where the rows of c m begin are stored in the array ia. i.e., if m(i,j) is the first c (leftmost) entry in the i-th row and a(k) = m(i,j), then c ia(i) = k. moreover, the index in ja and a of the first location c following the last element in the last row is stored in ia(n+1). c thus, the number of entries in the i-th row is given by c ia(i+1) - ia(i), the nonzero entries of the i-th row are stored c consecutively in c a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1), c and the corresponding column indices are stored consecutively in c ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1). c for example, the 5 by 5 matrix c ( 1. 0. 2. 0. 0.) c ( 0. 3. 0. 0. 0.) c m = ( 0. 4. 5. 6. 0.) c ( 0. 0. 0. 7. 0.) c ( 0. 0. 0. 8. 9.) c would be stored as c - 1 2 3 4 5 6 7 8 9 c ---+-------------------------- c ia - 1 3 4 7 8 10 c ja - 1 3 2 2 3 4 4 4 5 c a - 1. 2. 3. 4. 5. 6. 7. 8. 9. . c c the strict upper (lower) triangular portion of the matrix c u (l) is stored in a similar fashion using the arrays iu, ju, u c (il, jl, l) except that an additional array iju (ijl) is used to c compress storage of ju (jl) by allowing some sequences of column c (row) indices to used for more than one row (column) (n.b., l is c stored by columns). iju(k) (ijl(k)) points to the starting c location in ju (jl) of entries for the kth row (column). c compression in ju (jl) occurs in two ways. first, if a row c (column) i was merged into the current row (column) k, and the c number of elements merged in from (the tail portion of) row c (column) i is the same as the final length of row (column) k, then c the kth row (column) and the tail of row (column) i are identical c and iju(k) (ijl(k)) points to the start of the tail. second, if c some tail portion of the (k-1)st row (column) is identical to the c head of the kth row (column), then iju(k) (ijl(k)) points to the c start of that tail portion. for example, the nonzero structure of c the strict upper triangular part of the matrix c d 0 x x x c 0 d 0 x x c 0 0 d x 0 c 0 0 0 d x c 0 0 0 0 d c would be represented as c - 1 2 3 4 5 6 c ----+------------ c iu - 1 4 6 7 8 8 c ju - 3 4 5 4 c iju - 1 2 4 3 . c the diagonal entries of l and u are assumed to be equal to one and c are not stored. the array d contains the reciprocals of the c diagonal entries of the matrix d. c c iii. additional storage savings c in snsfc, r and ic can be the same array in the calling c sequence if no reordering of the coefficient matrix has been done. c in snnfc, r, c, and ic can all be the same array if no c reordering has been done. if only the rows have been reordered, c then c and ic can be the same array. if the row and column c orderings are the same, then r and c can be the same array. z and c row can be the same array. c in snnsc or snntc, r and c can be the same array if no c reordering has been done or if the row and column orderings are the c same. z and b can be the same array. however, then b will be c destroyed. c c iv. parameters c following is a list of parameters to the programs. names are c uniform among the various subroutines. class abbreviations are -- c n - integer variable c f - real variable c v - supplies a value to a subroutine c r - returns a result from a subroutine c i - used internally by a subroutine c a - array c c class - parameter c ------+---------- c fva - a - nonzero entries of the coefficient matrix m, stored c - by rows. c - size = number of nonzero entries in m. c fva - b - right-hand side b. c - size = n. c nva - c - ordering of the columns of m. c - size = n. c fvra - d - reciprocals of the diagonal entries of the matrix d. c - size = n. c nr - flag - error flag. values and their meanings are -- c - 0 no errors detected c - n+k null row in a -- row = k c - 2n+k duplicate entry in a -- row = k c - 3n+k insufficient storage for jl -- row = k c - 4n+1 insufficient storage for l c - 5n+k null pivot -- row = k c - 6n+k insufficient storage for ju -- row = k c - 7n+1 insufficient storage for u c - 8n+k zero pivot -- row = k c nva - ia - pointers to delimit the rows of a. c - size = n+1. c nvra - ijl - pointers to the first element in each column in jl, c - used to compress storage in jl. c - size = n. c nvra - iju - pointers to the first element in each row in ju, used c - to compress storage in ju. c - size = n. c nvra - il - pointers to delimit the columns of l. c - size = n+1. c nvra - iu - pointers to delimit the rows of u. c - size = n+1. c nva - ja - column numbers corresponding to the elements of a. c - size = size of a. c nvra - jl - row numbers corresponding to the elements of l. c - size = jlmax. c nv - jlmax - declared dimension of jl. jlmax must be larger than c - the number of nonzeros in the strict lower triangle c - of m plus fillin minus compression. c nvra - ju - column numbers corresponding to the elements of u. c - size = jumax. c nv - jumax - declared dimension of ju. jumax must be larger than c - the number of nonzeros in the strict upper triangle c - of m plus fillin minus compression. c fvra - l - nonzero entries in the strict lower triangular portion c - of the matrix l, stored by columns. c - size = lmax. c nv - lmax - declared dimension of l. lmax must be larger than c - the number of nonzeros in the strict lower triangle c - of m plus fillin (il(n+1)-1 after snsfc). c nv - n - number of variables/equations. c nva - r - ordering of the rows of m. c - size = n. c fvra - u - nonzero entries in the strict upper triangular portion c - of the matrix u, stored by rows. c - size = umax. c nv - umax - declared dimension of u. umax must be larger than c - the number of nonzeros in the strict upper triangle c - of m plus fillin (iu(n+1)-1 after snsfc). c fra - z - solution x. c - size = n. c c ---------------------------------------------------------------- c c*** subroutine snroc c*** reorders rows of a, leaving row order unchanged c c c input parameters.. n, ic, ia, ja, a c output parameters.. ja, a, flag c c parameters used internally.. c nia - p - at the kth step, p is a linked list of the reordered c - column indices of the kth row of a. p(n+1) points c - to the first entry in the list. c - size = n+1. c nia - jar - at the kth step,jar contains the elements of the c - reordered column indices of a. c - size = n. c fia - ar - at the kth step, ar contains the elements of the c - reordered row of a. c - size = n. c integer ic(*), ia(*), ja(*), jar(*), p(*), flag real a(*), ar(*) c double precision a(*), ar(*) c c ****** for each nonempty row ******************************* do 5 k=1,n jmin = ia(k) jmax = ia(k+1) - 1 if(jmin .gt. jmax) go to 5 p(n+1) = n + 1 c ****** insert each element in the list ********************* do 3 j=jmin,jmax newj = ic(ja(j)) i = n + 1 1 if(p(i) .ge. newj) go to 2 i = p(i) go to 1 2 if(p(i) .eq. newj) go to 102 p(newj) = p(i) p(i) = newj jar(newj) = ja(j) ar(newj) = a(j) 3 continue c ****** replace old row in ja and a ************************* i = n + 1 do 4 j=jmin,jmax i = p(i) ja(j) = jar(i) 4 a(j) = ar(i) 5 continue flag = 0 return c c ** error.. duplicate entry in a 102 flag = n + k return end subroutine snsfc * (n, r, ic, ia,ja, jlmax,il,jl,ijl, jumax,iu,ju,iju, * q, ira,jra, irac, irl,jrl, iru,jru, flag) c*** subroutine snsfc c*** symbolic ldu-factorization of nonsymmetric sparse matrix c (compressed pointer storage) c c c input variables.. n, r, ic, ia, ja, jlmax, jumax. c output variables.. il, jl, ijl, iu, ju, iju, flag. c c parameters used internally.. c nia - q - suppose m* is the result of reordering m. if c - processing of the ith row of m* (hence the ith c - row of u) is being done, q(j) is initially c - nonzero if m*(i,j) is nonzero (j.ge.i). since c - values need not be stored, each entry points to the c - next nonzero and q(n+1) points to the first. n+1 c - indicates the end of the list. for example, if n=9 c - and the 5th row of m* is c - 0 x x 0 x 0 0 x 0 c - then q will initially be c - a a a a 8 a a 10 5 (a - arbitrary). c - as the algorithm proceeds, other elements of q c - are inserted in the list because of fillin. c - q is used in an analogous manner to compute the c - ith column of l. c - size = n+1. c nia - ira, - vectors used to find the columns of m. at the kth c nia - jra, step of the factorization, irac(k) points to the c nia - irac head of a linked list in jra of row indices i c - such that i .ge. k and m(i,k) is nonzero. zero c - indicates the end of the list. ira(i) (i.ge.k) c - points to the smallest j such that j .ge. k and c - m(i,j) is nonzero. c - size of each = n. c nia - irl, - vectors used to find the rows of l. at the kth step c nia - jrl of the factorization, jrl(k) points to the head c - of a linked list in jrl of column indices j c - such j .lt. k and l(k,j) is nonzero. zero c - indicates the end of the list. irl(j) (j.lt.k) c - points to the smallest i such that i .ge. k and c - l(i,j) is nonzero. c - size of each = n. c nia - iru, - vectors used in a manner analogous to irl and jrl c nia - jru to find the columns of u. c - size of each = n. c c internal variables.. c jlptr - points to the last position used in jl. c juptr - points to the last position used in ju. c jmin,jmax - are the indices in a or u of the first and last c elements to be examined in a given row. c for example, jmin=ia(k), jmax=ia(k+1)-1. c integer cend, qm, rend, rk, vj integer ia(*), ja(*), ira(*), jra(*), il(*), jl(*), ijl(*) integer iu(*), ju(*), iju(*), irl(*), jrl(*), iru(*), jru(*) integer r(*), ic(*), q(*), irac(*), flag c c ****** initialize pointers **************************************** np1 = n + 1 jlmin = 1 jlptr = 0 il(1) = 1 jumin = 1 juptr = 0 iu(1) = 1 do 1 k=1,n irac(k) = 0 jra(k) = 0 jrl(k) = 0 1 jru(k) = 0 c ****** initialize column pointers for a *************************** do 2 k=1,n rk = r(k) iak = ia(rk) if (iak .ge. ia(rk+1)) go to 101 jaiak = ic(ja(iak)) if (jaiak .gt. k) go to 105 jra(k) = irac(jaiak) irac(jaiak) = k 2 ira(k) = iak c c ****** for each column of l and row of u ************************** do 41 k=1,n c c ****** initialize q for computing kth column of l ***************** q(np1) = np1 luk = -1 c ****** by filling in kth column of a ****************************** vj = irac(k) if (vj .eq. 0) go to 5 3 qm = np1 4 m = qm qm = q(m) if (qm .lt. vj) go to 4 if (qm .eq. vj) go to 102 luk = luk + 1 q(m) = vj q(vj) = qm vj = jra(vj) if (vj .ne. 0) go to 3 c ****** link through jru ******************************************* 5 lastid = 0 lasti = 0 ijl(k) = jlptr i = k 6 i = jru(i) if (i .eq. 0) go to 10 qm = np1 jmin = irl(i) jmax = ijl(i) + il(i+1) - il(i) - 1 long = jmax - jmin if (long .lt. 0) go to 6 jtmp = jl(jmin) if (jtmp .ne. k) long = long + 1 if (jtmp .eq. k) r(i) = -r(i) if (lastid .ge. long) go to 7 lasti = i lastid = long c ****** and merge the corresponding columns into the kth column **** 7 do 9 j=jmin,jmax vj = jl(j) 8 m = qm qm = q(m) if (qm .lt. vj) go to 8 if (qm .eq. vj) go to 9 luk = luk + 1 q(m) = vj q(vj) = qm qm = vj 9 continue go to 6 c ****** lasti is the longest column merged into the kth ************ c ****** see if it equals the entire kth column ********************* 10 qm = q(np1) if (qm .ne. k) go to 105 if (luk .eq. 0) go to 17 if (lastid .ne. luk) go to 11 c ****** if so, jl can be compressed ******************************** irll = irl(lasti) ijl(k) = irll + 1 if (jl(irll) .ne. k) ijl(k) = ijl(k) - 1 go to 17 c ****** if not, see if kth column can overlap the previous one ***** 11 if (jlmin .gt. jlptr) go to 15 qm = q(qm) do 12 j=jlmin,jlptr if (jl(j) - qm) 12, 13, 15 12 continue go to 15 13 ijl(k) = j do 14 i=j,jlptr if (jl(i) .ne. qm) go to 15 qm = q(qm) if (qm .gt. n) go to 17 14 continue jlptr = j - 1 c ****** move column indices from q to jl, update vectors *********** 15 jlmin = jlptr + 1 ijl(k) = jlmin if (luk .eq. 0) go to 17 jlptr = jlptr + luk if (jlptr .gt. jlmax) go to 103 qm = q(np1) do 16 j=jlmin,jlptr qm = q(qm) 16 jl(j) = qm 17 irl(k) = ijl(k) il(k+1) = il(k) + luk c c ****** initialize q for computing kth row of u ******************** q(np1) = np1 luk = -1 c ****** by filling in kth row of reordered a *********************** rk = r(k) jmin = ira(k) jmax = ia(rk+1) - 1 if (jmin .gt. jmax) go to 20 do 19 j=jmin,jmax vj = ic(ja(j)) qm = np1 18 m = qm qm = q(m) if (qm .lt. vj) go to 18 if (qm .eq. vj) go to 102 luk = luk + 1 q(m) = vj q(vj) = qm 19 continue c ****** link through jrl, ****************************************** 20 lastid = 0 lasti = 0 iju(k) = juptr i = k i1 = jrl(k) 21 i = i1 if (i .eq. 0) go to 26 i1 = jrl(i) qm = np1 jmin = iru(i) jmax = iju(i) + iu(i+1) - iu(i) - 1 long = jmax - jmin if (long .lt. 0) go to 21 jtmp = ju(jmin) if (jtmp .eq. k) go to 22 c ****** update irl and jrl, ***************************************** long = long + 1 cend = ijl(i) + il(i+1) - il(i) irl(i) = irl(i) + 1 if (irl(i) .ge. cend) go to 22 j = jl(irl(i)) jrl(i) = jrl(j) jrl(j) = i 22 if (lastid .ge. long) go to 23 lasti = i lastid = long c ****** and merge the corresponding rows into the kth row ********** 23 do 25 j=jmin,jmax vj = ju(j) 24 m = qm qm = q(m) if (qm .lt. vj) go to 24 if (qm .eq. vj) go to 25 luk = luk + 1 q(m) = vj q(vj) = qm qm = vj 25 continue go to 21 c ****** update jrl(k) and irl(k) *********************************** 26 if (il(k+1) .le. il(k)) go to 27 j = jl(irl(k)) jrl(k) = jrl(j) jrl(j) = k c ****** lasti is the longest row merged into the kth *************** c ****** see if it equals the entire kth row ************************ 27 qm = q(np1) if (qm .ne. k) go to 105 if (luk .eq. 0) go to 34 if (lastid .ne. luk) go to 28 c ****** if so, ju can be compressed ******************************** irul = iru(lasti) iju(k) = irul + 1 if (ju(irul) .ne. k) iju(k) = iju(k) - 1 go to 34 c ****** if not, see if kth row can overlap the previous one ******** 28 if (jumin .gt. juptr) go to 32 qm = q(qm) do 29 j=jumin,juptr if (ju(j) - qm) 29, 30, 32 29 continue go to 32 30 iju(k) = j do 31 i=j,juptr if (ju(i) .ne. qm) go to 32 qm = q(qm) if (qm .gt. n) go to 34 31 continue juptr = j - 1 c ****** move row indices from q to ju, update vectors ************** 32 jumin = juptr + 1 iju(k) = jumin if (luk .eq. 0) go to 34 juptr = juptr + luk if (juptr .gt. jumax) go to 106 qm = q(np1)