IOMATJ = J-1+1OMAT

lUMATJ = J-l+IUMAT

AD(J) = A(IOMATJ)

CD(J) = C(IUMATJ)

A(IOMATJ) =0.0

C(IUMATJ) = 0.0 20 CONTINUE

CALL MULPUT(AD,D(IIVEC),D(I0VEC),ORDER)

DO 22 J=l,ORDER

IOMATJ = IOMAT+(J-1)*ORDER

IIMATJ = I1MAT+(J-1)*0RDER

CALL MULPUT(AD,C(I1MATJ),B(IOMATJ),ORDER)

CALL MULPUT(AD,A(IIMATJ),A(IOMATJ)/ORDER) 22 CONTINUE

CALL LUDECO(B(I0MAT),ORDER)

CALL LUSOLV(B(IOMAT),D(IOVEC)jD(IOVEC),ORDER)

DO 24 J=l,ORDER

IOMATJ = IOMAT+(J-1)*ORDER

CALL LUSOLV(B(IOMAT),C(IOMATJ),C(IOMATJ),ORDER)

CALL LUSOLV(B(I0MAT),A(IOMATJ),A(IOMATJ),ORDER) 24 CONTINUE

CALL MULPUT(CD,D(I1VEC),D(IUVEC),ORDER)

DO 26 J=l,ORDER

lUMATJ = IUMAT+(J-1)*0RDER

IIMATJ = I1MAT+(J-1)*0RDER

CALL MULPUT(CD,A(I1MATJ),B(IUMATJ).ORDER)

CALL MULPUT(CD,С(IIMATJ),С(lUMATJ),ORDER) 26 CONTINUE 200 CONTINUE C...

DO 30 J=l,ORDSQ

lUMATJ = J-l+IUMAT

AD(J) = A(IUMATJ)+C(IUMATJ) 30 CONTINUE

CALL MULPUT(AD,D(IEVEC),D(IUVEC),ORDER)

DO 32 J=l,ORDER

lUMATJ = IUMAT+(J-1)*0RDER

lEMATJ = IEMAT+(J-1)*0RDER

CALL MULPUT(AD,C(IEMATJ),B(IUMATJ),ORDER)

CALL MULPUT(AD,A(IEMATJ),B(IUMATJ),ORDER) 32 CONTINUE

CALL LUDECO(B(lUMAT),ORDER)

CALL LUSOLV(B(IUMAT),D(IUVEC),D(IUVEC).ORDER)

C...

C...BACK SUBSTITUTION C...

DO 40 IBAC = IL,IE I = lE-IBAC+IL lOMAT = 1+(I-1)*0RDSQ lOVEC = 1+(I-1)*0RDER IIVEC = lOVEC+ORDER

CALL MULPUT(A(lOMAT),D(IUVEC),D(IOVEC).ORDER) CALL MULPUT(C(IOMAT),D(I1VEC),D(IOVEC).ORDER) 40 C0NTIN4JE C...

с...

с...SUBROUTINE ТО CALCULATE L-U DECOMPOSITION С...OF А GIVEN MATRIX А AND STORE RESULT IN A C...(NO PIVOTING STRATEGY IS EMPLOYED) C...

SUBROUTINE LUDECO(A,ORDER)

DIMENSION A(0RDER,1) INTEGER ORDER

DO"8 JC=2,ORDER 8 A(1,JC) = A(1,JC)/A(1,1)

JRJC = 1 10 CONTINUE

JRJC = JRJC+1

JRJCMl = JRJC-1

JRJCPr= JRJC+1

DO 14 JR=JRJC,ORDER

SUM = A(JR,JRJC)

DO 12 JM=1,JRJCMl 12 SUM = SUM-A(JR,JM)*A(JM,JRJC) 14 A(JR,JRJC) = SUM

IF (JRJC.EQ.ORDER) RETURN

DO 18 JC = JRJCPl,ORDER

SUM = A(JRJC,JC)

DO 16 JM=1,JRJCMl 16 SUM = SUM-A(JRJC,JM)*A(JM,JC) 18 A(JRJC,JC) = SUM/A(JRJC,JRJC)

GO TO 10

c...

C...SUBROUTINE TO MULTIPLY A VECTOR В BY A MATRIX A, C*. . . SUBTRACT RESULT FROM ANOTHER VECTOR С AND STORE C...RESULT IN C. THUS VECTOR С IS OVERWRITTEN. C...

SUBROUTINE MULPUT(A,B,C,ORDER)

DIMENSION A(1),B(1),C(1) INTEGER ORDER

DO 200 JR=1,ORDER

SUM = 0.0

DO 100 JC=1,ORDER

lA = JR+(JC-1)"ORDER 100 SUM = SUM+A(IA)VB(JC) 200 C(JR) = C(JR)-SUM

c...

RETURN END

с...

с...SUBROUTINE то SOLVE LINEAR ALGEBRAIC SYSTEM OF С...EQUATIONS A*C=B AND STORE RESULTS IN VECTOR C. C...MATRIX A IS INPUT IN L-U DECOMPOSITION FORM. C...(NO PIVOTING STRATEGY HAS BEEN EMPLOYED TO C...COMPUTE THE L-U DECOMPOSITION OF THE MATRIX A).

SUBROUTINE LUSOLV(A,B,C,ORDER)

DIMENSION A(0RDER,1),B(1),C(1) INTEGER ORDER

C...

C...FIRST L(INV)*B C. . .

C(l) » C(l)/A(l,l)

DO 14 JR=2,ORDER

JRMl s JR-1

SUM = B(JR)

DO 12 JM=1,JRM1 12 SUM a SUM-A(JR,JM)*C(JM) 14 C(JR) = SUM/A(JR,JR) C...

C...NEXT U(INV) OF L(INV)*B •

•.DO 18 JRJR=2,ORDER JR « ORDER-JRJR+1 JRPl « JR+1 SUM = C(JR)

DO 16 JMJM = JRPl,ORDER JM = ORDER-JMJM+JRPl

16 SUM = SUM-A(JR,JM)VC(JM)

18 C(JR) * SUM

RETURN END