SUBROUTINE ADISET(IM,IN,IMU,LM,LX,RM,RX,KK,IFN,IFUT,INFO,PARM,ARR) C*********************************************************************** C SUBROUTINE SETS UP COEFFICIENT ARRAYS FOR 2 STATE ADI ALGORITHM C ARGUMENTS: IM NUMBER OF GRID INTERVALS IN L DIRECTION C IN NUMBER OF GRID INTERVALS IN R DIRECTION C IMU FIRST DIMENSION - 1 OF U SOLN. MATRIX IF C DIFFERENT FROM IM (DEFAULTS TO IM) C LM MIN. VALUE OF STATE VARIABLE L C LX MAX. VALUE OF STATE VARIABLE L C RM MIN. VALUE OF STATE VARIABLE R C RX MAX. VALUE OF STATE VARIABLE R C KK STEP SIZE IN TIME DIRECTION C IFN MODEL NUMBER (IFNSET) SELECTING COEFFICIENT C IFUT FLAG SET TO 1 FOR FUTURES PRICES C INFO() FLAGS FOR BOUNDARY CONDITIONS (IULMN,IULMX, C IURMN,IURMX,ILMNRN,ILMNRX,ILMXRN,ILMXRX) C PARM() VECTOR OF MODEL PARAMETERS FOR CALCULATING COEF C ARRAY() OUTPUT ARRAY OF COEFFICIENTS SET FOR ADI STEP C ARRAY() CORRESPONDS TO BK1001 AS FOLLOWS: C ARR(1,.,.) = D1MTX(.,.) ARR(4,.,.) = DP1MTX() C ARR(2,.,.) = D2MTX(.,.) ARR(5,.,.) = DP2MTX() C ARR(3,.,.) = D3MTX(.,.) ARR(6,.,.) = DP3MTX() C ARR(7,.,.) = C1MTX(.,.) ARR(8,.,.) = C2MTX () C C AUTHOR: ROBERT JONES 3 OCTOBER 1988 C C*********************************************************************** IMPLICIT INTEGER (I-J,M-N) IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z) DIMENSION INFO(8), PARM(15) DIMENSION ARR( 8, 0:IM, 0:IN ) COMMON /ADICOM/ A,B,C,D,E,F,K, X1,X2,X3,X4,X5, 1 L,LMIN,LMAX,LSCALE, R,RMIN,RMAX,RSCALE, 2 IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX, 3 ILMXRN,ILMXRX,IFNSET,IFUTUR, M,N,MU SAVE ADICOM M = IM N = IN K = KK MU = MAX( IM, IMU ) IFNSET = IFN IFUTUR = IFUT IULMN = INFO(1) IULMX = INFO(2) IURMN = INFO(3) IURMX = INFO(4) ILMNRN = INFO(5) ILMNRX = INFO(6) ILMXRN = INFO(7) ILMXRX = INFO(8) LMIN = LM LMAX = LX RMIN = RM RMAX = RX DO 100 I3 = 0, N DO 100 I2 = 0, M DO 100 I1 = 1, 8 100 ARR(I1,I2,I3) = 0.0 LSCALE = DBLE(M) / ( LMAX - LMIN ) RSCALE = DBLE(N) / ( RMAX - RMIN ) X1 = LSCALE * LSCALE * K X2 = RSCALE * RSCALE * K X3 = LSCALE * RSCALE * K X4 = LSCALE * K X5 = RSCALE * K C**** COMPUTE INTERIOR COEFFICIENTS DO 200 I = 1, M-1 DO 200 J = 1, N-1 CALL ADI1(I,J,PARM) ARR(1,I,J) = D / 4.0 - A / 2.0 ARR(2,I,J) = 1.0 + A ARR(3,I,J) = -D / 4.0 - A / 2.0 ARR(7,I,J) = C / 4.0 ARR(8,I,J) = F - 2.0 * B - A ARR(4,I,J) = E / 4.0 - B / 2.0 ARR(5,I,J) = 1.0 + B ARR(6,I,J) = -E / 4.0 - B / 2.0 200 CONTINUE C**** COMPUTE COEFFICIENTS AT LMIN IF ( IULMN .EQ. 1 ) THEN C KNOWN SOLUTION AT LMIN DO 220 J = 1, N-1 220 ARR(2,0,J) = 1.0 ELSE IF ( IULMN .EQ. 2 ) THEN C SOLN AT LMIN SATISFIES 1 STATE VARIABLE PDE DO 240 J = 1, N-1 240 ARR(2,0,J) = 1.0 CALL ADI2(PARM,ARR) ELSE C QUADRATIC EXTRAPOL. DO 260 J = 1, N-1 G = ( 3.0 + ARR(2,2,J) / ARR(3,2,J) ) / ARR(3,1,J) ARR(3,0,J) = -3.0 + ARR(1,2,J)/ARR(3,2,J)-G*ARR(2,1,J) 260 ARR(2,0,J) = 1.0 - G * ARR(1,1,J) CALL ADI2(PARM,ARR) END IF C**** COMPUTE COEFFICIENTS AT LMAX IF ( IULMX .EQ. 1 ) THEN C KNOWN SOLUTION AT L = MAX DO 310 J= 1,N-1 310 ARR(2,M,J)= 1.0 ELSE C FOR QUADRATIC EXTRAPOLATION DO 320 J = 1, N-1 G = (3.0+ARR(2,M-2,J)/ARR(1,M-2,J)) / ARR(1,M-1,J) ARR(1,M,J)= -3.0+ARR(3,M-2,J)/ARR(1,M-2,J)-G*ARR(2,M-1,J) 320 ARR(2,M,J)= 1.0 - G * ARR(3,M-1,J) DO 330 J= 1,N-1 CALL ADI1(M,J,PARM) ARR(4,M,J) = E/4.0 - B/2.0 ARR(5,M,J) = 1.0 + B ARR(6,M,J) = -E/4.0 - B/2.0 C USED RHS OF BOUNDARY PDE 330 ARR(8,M,J) = B - F IF ( ILMXRN.EQ.1) THEN C SOLN. KNOWN AT L=LMAX, R=0: ARR(5,M,0) = 1.0 ELSE G = ( 3.0 + ARR(5,M,2)/ARR(6,M,2) ) / ARR(6,M,1) ARR(6,M,0) = -3.0 + ARR(4,M,2)/ARR(6,M,2) -G * ARR(5,M,1) ARR(5,M,0) = 1.0 - G * ARR(4,M,1) END IF IF ( ILMXRX .EQ. 1 ) THEN C KNOWN SOLUTION AT L=LMAX, R=RMAX: ARR(5,M,N) = 1.0 ELSE C QUADRATIC EXTRAPOLATION: G = ( 3.0 + ARR(5,M,N-2)/ARR(4,M,N-2) ) / ARR(4,M,N-1) ARR(4,M,N) =-3.0+ARR(6,M,N-2)/ARR(4,M,N-2)-G*ARR(5,M,N-1) ARR(5,M,N) = 1.0 - G * ARR(6,M,N-1) END IF END IF C**** COMPUTE COEFFICIENTS AT RMIN IF ( IURMN .EQ. 1 ) THEN C KNOWN SOLUTION AT R=RMIN: DO 410 I = 1, M-1 410 ARR(5,I,0) = 1.0 ELSE IF ( IURMN .EQ. 2 ) THEN C FOR BOUNDARY PDE SOLUTION: DO 420 I = 1, M-1 420 ARR(5,I,0) = 1.0 DO 430 I = 1, M-1 CALL ADI1(I,0,PARM) ARR(1,I,0) = D/4.0 - A/2.0 ARR(2,I,0) = 1.0 + A ARR(3,I,0) = -D/4.0 - A/2.0 C USED RHS OF BOUNDARY PDE 430 ARR(8,I,0) = A - F IF (ILMNRN.EQ.1) THEN C KNOWN SOLUTION AT L=LMIN, R=RMIN: ARR(2,0,0)= 1.0 ELSE G= (3.0+ARR(2,2,0)/ARR(3,2,0))/ARR(3,1,0) ARR(3,0,0)= -3.0+ARR(1,2,0)/ARR(3,2,0)-G*ARR(2,1,0) ARR(2,0,0)= 1.0-G*ARR(1,1,0) END IF IF ( ILMXRN .EQ. 1 ) THEN C KNOWN SOLUTION AT L=LMAX, R=RMIN: ARR(2,M,0) = 1.0 ELSE G = (3.0+ARR(2,M-2,0)/ARR(1,M-2,0))/ARR(1,M-1,0) ARR(1,M,0) = -3.0+ARR(3,M-2,0)/ARR(1,M-2,0)-G / *ARR(2,M-1,0) ARR(2,M,0) = 1.0 - G * ARR(3,M-1,0) END IF ELSE C QUADRATIC EXTRAPOLATION: DO 440 I = 1, M-1 G = (3.0+ARR(5,I,2)/ARR(6,I,2))/ARR(6,I,1) ARR(6,I,0) = -3.0 + ARR(4,I,2)/ARR(6,I,2)- G*ARR(5,I,1) 440 ARR(5,I,0) = 1.0 - G * ARR(4,I,1) END IF C**** COMPUTE COEFFICIENTS AT RMAX IF ( IURMX .EQ. 1 ) THEN C KNOWN SOLUTION AT R=MAX: DO 500 I = 1, M-1 500 ARR(5,I,N) = 1.0 ELSE C QUADRATIC EXTRAPOLATION: DO 510 I = 1, M-1 G = ( 3.0 + ARR(5,I,N-2)/ARR(4,I,N-2) ) / ARR(4,I,N-1) ARR(4,I,N) = -3.0 + ARR(6,I,N-2) / / ARR(4,I,N-2) - G * ARR(5,I,N-1) 510 ARR(5,I,N) = 1.0 - G * ARR(6,I,N-1) END IF RETURN END C**** SUBROUTINES CALLED BY ADISET() ********************************** C**** SUBROUTINE TO CALCULATE & SCALE PDE COEFFICIENTS AT (I,J): SUBROUTINE ADI1(I,J,PARM) IMPLICIT INTEGER (I-J,M-N) IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z) DIMENSION PARM(15) COMMON /ADICOM/ A,B,C,D,E,F,K, X1,X2,X3,X4,X5, 1 L,LMIN,LMAX,LSCALE, R,RMIN,RMAX,RSCALE, 2 IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX, 3 ILMXRN,ILMXRX,IFNSET,IFUTUR, M,N,MU SAVE ADICOM L = LMIN + DBLE(I) / LSCALE R = RMIN + DBLE(J) / RSCALE IF ( IFUTUR .EQ. 1 ) THEN F = 0.0 ELSE F = FNAKF(IFNSET,L,R,PARM) END IF K1 = 1.0 - F * K / 2.0 A = X1 * FNAKA(IFNSET,L,R,PARM) / K1 B = X2 * FNAKB(IFNSET,L,R,PARM) / K1 C = X3 * FNAKC(IFNSET,L,R,PARM) / K1 D = X4 * FNAKD(IFNSET,L,R,PARM) / K1 E = X5 * FNAKE(IFNSET,L,R,PARM) / K1 F = (2.0 - K1) / K1 RETURN END C**** FORMERLY SUBROUTINE AK2014(IFUTUR,PARM) SUBROUTINE ADI2(PARM,ARR) IMPLICIT INTEGER (I-J,M-N) IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z) COMMON /ADICOM/ A,B,C,D,E,F,K, X1,X2,X3,X4,X5, 1 L,LMIN,LMAX,LSCALE, R,RMIN,RMAX,RSCALE, 2 IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX, 3 ILMXRN,ILMXRX,IFNSET,IFUTUR, M,N,MU SAVE ADICOM DIMENSION PARM(15), ARR( 8 , 0:M , 0:* ) DO 100 J = 1, N-1 CALL ADI1(0,J,PARM) ARR(4,0,J) = E/4.0 - B/2.0 ARR(5,0,J) = 1.0 + B ARR(6,0,J) = -E/4.0 - B/2.0 C USED RHS OF BOUNDARY 100 ARR(8,0,J) = B-F IF ( ILMNRN .EQ. 1 ) THEN C IF GIVEN VALUE AT L=0, R=0: ARR(5,0,0) = 1.0 ELSE C QUADRATIC EXTRAPLATION: G = ( 3.0 + ARR(5,0,2)/ARR(6,0,2) ) / ARR(6,0,1) ARR(6,0,0) = -3.0 + ARR(4,0,2)/ARR(6,0,2) - G * ARR(5,0,1) ARR(5,0,0) = 1.0 - G * ARR(4,0,1) END IF IF ( ILMNRX .EQ. 1 ) THEN C IF GIVEN VALUE AT L=0, R=RMAX: ARR(5,0,N) = 1.0 ELSE C QUADRATIC EXTRAPLATION: G = ( 3.0 + ARR(5,0,N-2)/ARR(4,0,N-2) ) / ARR(4,0,N-1) ARR(4,0,N) = -3.0 + ARR(6,0,N-2)/ARR(4,0,N-2)-G*ARR(5,0,N-1) ARR(5,0,N) = 1.0 - G * ARR(6,0,N-1) END IF RETURN END SUBROUTINE ADSTEP(U,V,ARR) C*********************************************************************** C SUBROUTINE TAKES ONE TIME STEP USING 2 STATE ADI ALGORITHM. C CALL ADISET() FIRST TO SET UP NECESSARY COEFFICIENT ARRAYS ARR(). C ARGUMENTS: U() CURRENT SOLUTION ARRAY -- OVERWRITTEN BY ADSTEP C V() INTERMEDIATE SOLUTION ARRAY C ARR() INPUT ARRAY OF COEFFICIENTS SET BY ADISET C ARRAY() CORRESPONDS TO BK1001 AS FOLLOWS: C ARR(1,.,.) = D1MTX(.,.) ARR(4,.,.) = DP1MTX() C ARR(2,.,.) = D2MTX(.,.) ARR(5,.,.) = DP2MTX() C ARR(3,.,.) = D3MTX(.,.) ARR(6,.,.) = DP3MTX() C ARR(7,.,.) = C1MTX(.,.) ARR(8,.,.) = C2MTX () C AUTHOR: ROBERT JONES 3 OCTOBER 1988 C*********************************************************************** IMPLICIT INTEGER (I-J,M-N) IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z) C NOTE: INFORMATION IN ADICOM IS SET BY SUBROUTINE ADISET() C VARIABLES IN FIRST ROW OF ADICOM NOT USED IN ADSTEP() COMMON /ADICOM/ ZA,ZB,ZC,ZD,ZE,ZF,ZK, Z1,Z2,Z3,Z4,Z5, 1 L,LMIN,LMAX,LSCALE, R,RMIN,RMAX,RSCALE, 2 IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX, 3 ILMXRN,ILMXRX,IFNSET,IFUTUR, M,N,MU SAVE ADICOM C NOTE: PARAMETER NT MUST BE .GE. MAX(M,N) FOR TRIDAG ALGORITHM PARAMETER ( NT = 200 ) COMMON /TRICOM/ A(0:NT),B(0:NT),C(0:NT),D(0:NT),GAM(0:NT),X(0:NT) C NOTE: ARR,U,V ARE THE ADJUSTABLE DIMENSION ARRAYS DIMENSION ARR(8, 0:M, 0:N), U(0:MU, 0:N), V(0:M, 0:N) C**** FIRST SOLVE IN R DIRECTION *************************************** IF ( IULMN .EQ. 1 ) THEN C PUT LMIN BOUNDARY VALUES IN V(0,J): DO 100 J = 0, N 100 V(0,J) = FNULMN( RMIN + DBLE(J)/RSCALE ) ENDIF IF ( IULMN .EQ. 2 ) THEN C SOLVE BOUNDARY PDE AT L=0. PUT TRANSFORMED SOLN. IN V(0,J) DO 110 J = 0, N A(J) = ARR(4,0,J) B(J) = ARR(5,0,J) 110 C(J) = ARR(6,0,J) DO 120 J=1,N-1 120 D(J) = -A(J)*U(0,J-1) - ARR(8,0,J)*U(0,J) - C(J)*U(0,J+1) IF ( ILMNRN .EQ. 1 ) THEN D(0) = FNLNRN ELSE D(0) = D(2) / C(2) - D(1) * ( 3.0 + B(2) / C(2) ) / C(1) ENDIF IF ( ILMNRX .EQ. 1 ) THEN D(N) = FNLNRX ELSE D(N) = D(N-2)/A(N-2)-D(N-1)*(3.0 + B(N-2)/A(N-2))/A(N-1) ENDIF CALL ADI4(A,B,C,D,GAM,X,N) C TRANSFORM RESULT BACK TO INTERMEDIATE SOLN VALUES DO 130 J = 1, N-1 130 V(0,J) = U(0,J) + ARR(4,0,J) * ( X(J-1) - U(0,J-1) ) / + ARR(5,0,J) * ( X(J) - U(0,J) ) / + ARR(6,0,J) * (X(J+1) - U(0,J+1)) V(0,0) = U(0,0) + ARR(5,0,0) * (X(0)-U(0,0)) + / ARR(6,0,0) * (X(1)-U(0,1)) V(0,N) = U(0,N) + ARR(4,0,N) * (X(N-1) - U(0,N-1)) + / ARR(5,0,N) * (X(N) - U(0,N)) ENDIF IF ( IULMX .EQ. 1 ) THEN C PUT UMAX BOUNDARY VALUES IN V(M,J) DO 140 J = 0, N 140 V(M,J) = FNULMX( RMIN + DBLE(J)/RSCALE ) ENDIF C SOLVE FOR INTERMEDIATE SOLUTION VALUES DO 180 J = 1, N-1 DO 150 I = 0, M A(I) = ARR(1,I,J) B(I) = ARR(2,I,J) C(I) = ARR(3,I,J) 150 CONTINUE DO 160 I = 1, M-1 D(I) = ( U(I+1,J+1) - U(I+1,J-1) - U(I-1,J+1) + U(I-1,J-1) ) / * ARR(7,I,J) / - U(I+1,J) * ARR(3,I,J) / - U(I-1,J) * ARR(1,I,J) / + U(I,J) * ARR(8,I,J) / - U(I,J+1) * ARR(6,I,J) * 2.0 / - U(I,J-1) * ARR(4,I,J) * 2.0 160 CONTINUE IF ( IULMN .EQ. 1 .OR. IULMN .EQ. 2 ) THEN D(0) = V(0,J) ELSE D(0) = D(2) / C(2) - D(1) * ( 3.0 + B(2)/C(2) ) / C(1) ENDIF IF ( IULMX .EQ. 1 ) THEN D(M) = V(M,J) ELSE D(M) = D(M-2)/A(M-2)-D(M-1)*(3.0 + B(M-2)/A(M-2))/A(M-1) ENDIF CALL ADI4(A,B,C,D,GAM,X,M) DO 170 I = 0, M 170 V(I,J) = X(I) 180 CONTINUE C**** ALTERNATE DIRECTION TO GET FINAL SOLUTION VALUES ***************** IF ( IURMN .EQ. 1 ) THEN C PUT RMIN BOUNDARY VALUES IN V(I,0) DO 200 I = 0, M 200 V(I,0) = FNURMN( LMIN + DBLE(I)/LSCALE ) ENDIF IF ( IURMN .EQ. 2 ) THEN C SOLVE BOUNDARY PDE AT R=0. STORE SOLUTION IN V(I,0) DO 210 I = 0, M A(I) = ARR(1,I,0) B(I) = ARR(2,I,0) 210 C(I) = ARR(3,I,0) DO 220 I=1,M-1 220 D(I)= - A(I)*U(I-1,0) - ARR(8,I,0)*U(I,0) - C(I)*U(I+1,0) IF ( ILMNRN .EQ. 1 ) THEN D(0) = FNLNRN ELSE D(0) = D(2)/C(2) - D(1) * ( 3.0+B(2) / C(2) ) / C(1) ENDIF IF ( ILMXRN .EQ. 1 ) THEN D(M) = FNLXRN ELSE D(M) = D(M-2)/A(M-2) - D(M-1)*(3.0+B(M-2)/A(M-2)) / A(M-1) ENDIF CALL ADI4(A,B,C,D,GAM,X,M) DO 230 I = 0, M 230 V(I,0) = X(I) ENDIF IF ( IURMX .EQ. 1 ) THEN C PUT RMAX BOUNDARY VALUES IN V(I,N) DO 240 I = 0, M 240 V(I,N) = FNURMX( LMIN + DBLE(I)/LSCALE ) ENDIF C SOLVE FOR FINAL U VALUES DO 250 I=1,M-1 250 CALL ADI3(I,U,V,ARR) IF ( IULMN .EQ. 1 ) THEN C COPY V TO U: DO 260 J = 0, N 260 U(0,J) = V(0,J) ELSEIF ( IULMN .EQ. 2 ) THEN CALL ADI3(0,U,V,ARR) ELSE C QUADRATIC EXTRAPOLATION AT LMIN: DO 270 J = 0, N 270 U(0,J) = 3.0 * U(1,J) - 3.0 * U(2,J) + U(3,J) ENDIF IF ( IULMX .EQ. 1 ) THEN C COPY V VALUES TO U FOR I=M: DO 280 J = 0, N 280 U(M,J) = V(M,J) ELSEIF ( IULMX .EQ. 2 ) THEN CALL ADI3(M,U,V,ARR) ELSE C QUADRATIC EXTRAPOLATION AT LMAX: DO 290 J = 0, N 290 U(M,J) = 3.0 * U(M-1,J) - 3.0 * U(M-2,J) + U(M-3,J) ENDIF RETURN END C**** SUBROUTINES CALLED BY ADSTEP (PREVIOUSLY BK4301 & BK4500 ********* SUBROUTINE ADI3(I,U,V,ARR) IMPLICIT INTEGER (I-J,M-N) IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z) COMMON /ADICOM/ ZA,ZB,ZC,ZD,ZE,ZF,ZK, Z1,Z2,Z3,Z4,Z5, 1 L,LMIN,LMAX,LSCALE, R,RMIN,RMAX,RSCALE, 2 IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX, 3 ILMXRN,ILMXRX,IFNSET,IFUTUR, M,N,MU SAVE ADICOM DIMENSION ARR(8, 0:M, 0:N), U(0:MU, 0:N), V(0:M, 0:N) C NOTE: PARAMETER NT MUST BE .GE. MAX(M,N) FOR TRIDAG ALGORITHM PARAMETER ( NT = 200 ) COMMON /TRICOM/ A(0:NT),B(0:NT),C(0:NT),D(0:NT),GAM(0:NT),X(0:NT) DO 100 J = 0, N A(J) = ARR(4,I,J) B(J) = ARR(5,I,J) 100 C(J) = ARR(6,I,J) DO 110 J = 1, N-1 D(J) = V(I,J) + ARR(4,I,J) * U(I,J-1) / + (ARR(5,I,J) - 1.0) * U(I,J) + ARR(6,I,J) * U(I,J+1) 110 CONTINUE IF ( IURMN .EQ. 1 .OR. IURMN .EQ. 2 ) THEN D(0) = V(I,0) ELSE D(0) = D(2)/C(2) - D(1) * ( 3.0 + B(2)/C(2) ) / C(1) ENDIF IF ( IURMX .EQ. 1 ) THEN D(N) = V(I,N) ELSE D(N) = D(N-2)/A(N-2) - D(N-1) * (3.0 + B(N-2)/A(N-2)) / A(N-1) ENDIF CALL ADI4(A,B,C,D,GAM,X,N) DO 120 J = 0, N 120 U(I,J) = X(J) RETURN END C**** TRIDIAGONAL SOLN. ALGORITHM FROM "NUMERICAL RECIPES", P. 40 ****** C SOLVES: (ABC)X = D FOR X. N=DIMENSION. A,B,C,D, NOT ALTERED C NOTE: SUBSCRIPTS RUN FROM 0 AND SCRATCH VECTOR GAM VARIABLE DIMEN. SUBROUTINE ADI4(A,B,C,D,GAM,X,N) IMPLICIT INTEGER (I-J,M-N) IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z) DIMENSION A(0:*),B(0:*),C(0:*),D(0:*),GAM(0:*),X(0:*) BET = B(0) X(0) = D(0) / BET DO 11 J = 1, N GAM(J) = C(J-1) / BET BET = B(J) - A(J) * GAM(J) X(J) = (D(J) - A(J) * X(J-1)) / BET 11 CONTINUE DO 12 J = N-1, 0, -1 12 X(J) = X(J) - GAM(J+1) * X(J+1) RETURN END