%PROCESS MARGINS(1,100); /* INVERSION OF POSITIVE DEFINITE MATRICES BY THE GAUSS-JORDAN METHOD */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 47. */ /* Variable IJ added to declarations in procedure GJDEF2, 22/3/2007 */ GJDEF1:PROC(N,A,FAIL) OPTIONS (REORDER); /*PAGE 47*/ /*00000010*/ DCL (N,I,J,K) FIXED BIN (31,0), /*00000020*/ (P,Q,HI) FLOAT BIN (21), /*00000030*/ FAIL LABEL, /*00000040*/ A(*,*) BIN FLOAT(21), /*00000050*/ H(LB:UB) BIN FLOAT(21) CTL, /*00000060*/ (LB,UB) BIN FIXED(31); /*00000070*/ LB=LBOUND(A,1); UB=LB+N-1; ALLOCATE H; /*00000080*/ DO K=UB BY -1 TO LB; /*00000090*/ P=A(LB,LB); /*00000100*/ IF P<=0 THEN GO TO FAIL; /*00000110*/ DO I=LB+1 TO UB; /*00000120*/ Q=A(I,LB); /*00000130*/ HI=Q/P; /*00000140*/ IF I<=K THEN HI=-HI; /*00000150*/ H(I)=HI; /*00000160*/ DO J=LB+1 TO I; /*00000170*/ A(I-1,J-1)=A(I,J)+Q*H(J); /*00000180*/ END; /*00000190*/ END; /*00000200*/ A(N,N)=1E0/P; /*00000210*/ DO I=LB+1 TO UB; /*00000220*/ A(N,I-1)=H(I); /*00000230*/ END; /*00000240*/ END; /*00000250*/ END GJDEF1; /*00000260*/ GJDEF2:PROC(N,A,FAIL) OPTIONS (REORDER); /*PAGE 47*/ /*00000010*/ DCL (N,I,II,IJ,K,M) FIXED BIN (31,0), /*00000020*/ (A(*), H(LB:UB) CTL, HI, /*00000030*/ (P,Q))FLOAT BIN (21), /*00000040*/ (LB,UB) BIN FIXED(31), /*00000050*/ FAIL LABEL; /*00000060*/ M=0; LB=LBOUND(A,1); UB=LB+N-1; ALLOCATE H; /*00000070*/ DO K=UB BY -1 TO LB; /*00000080*/ P=A(LB); /*00000090*/ IF P<=0E0 THEN GO TO FAIL; /*00000100*/ II=1; /*00000110*/ DO I=LB+1 TO UB; /*00000120*/ M=II;II=II+I;Q=A(M+1); /*00000130*/ HI=Q/P; /*00000140*/ IF I<=K THEN HI=-HI; /*00000150*/ H(I)=HI; /*00000160*/ DO IJ=M+2 TO II; /*00000170*/ A(IJ-I)=A(IJ)+Q*H(IJ-M); /*00000180*/ END; /*00000190*/ END; /*00000200*/ M=M-1;A(II)=1E0/P; /*00000210*/ DO I=LB+1 TO UB; /*00000220*/ A(M+I)=H(I); /*00000230*/ END; /*00000240*/ END; /*00000250*/ END GJDEF2; /*00000260*/