CTITLEMATINV -- MATRIX INVERSION 00010002 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA 00020000 CA AUTHOR U.C.S.D. COMPUTER CENTER 00030000 CA LANGUAGE VS FORTRAN 00040001 CA WRITTEN 00/00/00 00050000 C REVISED 05/01/86 REP. CONVERT TO VS FORTRAN 00060001 CA 00070000 CA CALL MATINV (A, N, NN, B, M) 00080000 CA 00090000 CA INPUT A = SYMMETRIC REAL MATRIX R4 00100000 CA N = ORDER OF MATRIX I4 00110000 CA NN = ACTUAL DIMENSION OF MATRIX A I4 00120000 CA B = ARRAY OF COLUMN VECTORS R4 00130000 CA M = NUMBER OF COLUMN VECTORS IN B, MAY BE 0 I4 00140000 CA 00150000 CA OUTPUT A = INVERSE OF INPUT MATRIX R4 00160000 CA B = LINEAR EQUATION SOLUTION(S) R4 00170000 CA 00180000 CA MATINV PERFORMS A MATRIX INVERSION AND LINEAR EQUATION SOLUTION 00190000 CA OF THE FORM A * X = B, WHERE A COLUMN VECTOR SOLUTION X IS 00200000 CA FOUND FOR EACH INPUT COLUMN VECTOR B. 00210000 C 00220000 C=======================================================================00230000 C EJECT 00240000 C 00250000 SUBROUTINE MATINV (A, N, NN, B, M) 00260000 C 00270000 DIMENSION A(NN,NN),B(NN,M),IP(100),IX(100,2) 00280000 EQUIVALENCE (IR,JR),(IC,JC),(AX,T,SW) 00290000 C 00300000 DO 10 J=1,N 00310000 10 IP(J)=0 00320000 C 00330000 DO 130 I=1,N 00340000 AX=0 00350000 C 00360000 DO 40 J=1,I 00370000 IF (IP(J).EQ.1) GO TO 40 00380000 C 00390000 DO 30 K=1,N 00400000 IF (IP(K)-1) 20 , 30 , 160 00410000 20 IF (ABS(AX).GE.ABS(A(J,K)).AND.AX.NE.0) GO TO 30 00420000 IR=J 00430000 IC=K 00440000 AX=A(J,K) 00450000 30 CONTINUE 00460000 C 00470000 40 CONTINUE 00480000 C 00490000 IP(IC)=IP(IC)+1 00500000 IF (IR.EQ.IC) GO TO 70 00510000 C 00520000 DO 50 L=1,N 00530000 SW=A(IR,L) 00540000 A(IR,L)=A(IC,L) 00550000 50 A(IC,L)=SW 00560000 C 00570000 IF (M.LE.0) GO TO 70 00580000 C 00590000 DO 60 L=1,M 00600000 SW=B(IR,L) 00610000 B(IR,L)=B(IC,L) 00620000 60 B(IC,L)=SW 00630000 C 00640000 70 IX(I,1)=IR 00650000 IX(I,2)=IC 00660000 PV=A(IC,IC) 00670000 A(IC,IC)=1. 00680000 C 00690000 DO 80 L=1,N 00700000 80 A(IC,L)=A(IC,L)/PV 00710000 C 00720000 IF (M.LE.0) GO TO 100 00730000 C 00740000 DO 90 L=1,M 00750000 90 B(IC,L)=B(IC,L)/PV 00760000 C 00770000 100 DO 130 L1=1,N 00780000 IF (L1.EQ.IC) GO TO 130 00790000 T=A(L1,IC) 00800000 A(L1,IC)=0 00810000 C 00820000 DO 110 L=1,N 00830000 110 A(L1,L)=A(L1,L)-A(IC,L)*T 00840000 C 00850000 IF (M.LE.0) GO TO 130 00860000 C 00870000 DO 120 L=1,M 00880000 120 B(L1,L)=B(L1,L)-B(IC,L)*T 00890000 C 00900000 130 CONTINUE 00910000 C 00920000 DO 150 I=1,N 00930000 L=N+1-I 00940000 IF (IX(L,1).EQ.IX(L,2)) GO TO 150 00950000 JR=IX(L,1) 00960000 JC=IX(L,2) 00970000 C 00980000 DO 140 K=1,N 00990000 SW=A(K,JR) 01000000 A(K,JR)=A(K,JC) 01010000 140 A(K,JC)=SW 01020000 C 01030000 150 CONTINUE 01040000 C 01050000 160 RETURN 01060000 C 01070000 END 01080000