C #INVTB V1A 17-APR-72 C LAST UPDATE 30-MAY-73. C SUBROUTINE INVTB(A,N,M,DET,WRK) INTEGER WRK DIMENSION A(M,M),WRK(N) C C *PURPOSE. C TO INVERT A REAL SQUARE MATRIX. C C *PARAMETERS. C C A - THE REAL ARRAY OF DIMENSION "M" X "M" CONTAINING THE C MATRIX OF SIZE "N" X "N" TO BE INVERTED. THE MATRIX C IS INVERTED IN PLACE. C N - THE LOGICAL SIZE OF THE MATRIX TO BE INVERTED. C M - THE PHYSICAL SIZE OF THE ARRAY CONTAINING THE MATRIX. C DET - A REAL VALUE RETURNING THE DETERMINANT OF THE ORIGINAL C MATRIX. C WRK - AN ARRAY CONTAINING AT LEAST "N" INTEGER ELEMENTS C USED AS WORKSPACE BY "INVTB". C C INPUT PARAMETERS: C A, N, M C OUTPUT PARAMETERS: C A, DET, WRK. C C *METHOD. C THE CROUT FORM OF GAUSSIAN ELIMINATION IS USED AS C DESIGNED TO INVERT THE MATRIX IN PLACE. PIVOTING IS DONE BY USING C THE LARGEST ELEMENT FOUND IN A SEARCH DOWN THE CURRENT COLUMN, AND C THE ROUTINE GUARDS AGAINST A ZERO PIVOT ELEMENT. C C *ACCURACY. C THE ACCURACY OF THIS METHOD DEPENDS VERY MUCH ON THE CONDITION C OF THE MATRIX AND GETS WORSE AS DET GETS SMALLER. HOWEVER THIS C ROUTINE WILL NORMALLY GIVE ACCURACY AT LEAST ONE ORDER OF MAGNITUDE C BETTER THAN INVTA AND IS LESS SENSITIVE TO SMALL VALUES OF DET THAN C THAT ROUTINE. C C *RESTRICTIONS. C C *ERROR CONDITIONS. C IF NO NON ZERO PIVOT CAN BE FOUND AT ANY STAGE THE C VALUE OF DET IS RETURNED AS MACHINE ZERO. C C *NON STANDARD ROUTINES CALLED. C C *TYPICAL TIMES. C THE TIME TAKEN IS ROUGHLY PROPORTIONAL TO N**3 AND THE EXTRA C TIME FOR PIVOTING IS GENERALY NEGLIGEABLE. C C *ORIGIN. M.R.MANNING. C C *COMMENTS. C AS ALREADY MENTIONED THIS ROUTINE SHOULD BE USED IN PREFERENCE C TO THE ROUTINE INVTA EXCEPT WHERE IT IS KNOWN THAT THAT ROUTINE WILL C PRODUCE ACCEPTABLE RESULTS. C C C #END. C C #START. DET=1. DO 1 K=1,N 1 WRK(K)=K DO 10 K=1,N FAC=0. DO 3 J=K,N IF(ABS(A(J,K)).LE.ABS(FAC)) GO TO 3 JST=J FAC=A(J,K) 3 CONTINUE IF(FAC) 5,4,5 4 DET=0. RETURN 5 DET=DET*FAC IF(JST.EQ.K) GO TO 7 J=WRK(JST) WRK(JST)=WRK(K) WRK(K)=J DO 6 I=1,N TEM=A(JST,I) A(JST,I)=A(K,I) 6 A(K,I)=TEM 7 A(K,K)=1./FAC DO 8 I=1,N IF(I.EQ.K) GO TO 8 A(K,I)=A(K,I)/FAC 8 CONTINUE DO 10 I=1,N IF(I.EQ.K) GO TO 10 AFAC=A(I,K) A(I,K)=-AFAC/FAC DO 9 J=1,N IF(J.EQ.K) GO TO 9 A(I,J)=A(I,J)-AFAC*A(K,J) 9 CONTINUE 10 CONTINUE C DO 20 I=1,N 15 K=WRK(I) IF(I.EQ.K) GO TO 20 DET=-DET WRK(I)=WRK(K) WRK(K)=K DO 18 J=1,N TEM= A(J, I) A(J, I)= A(J, K) 18 A(J, K)= TEM GO TO 15 C 20 CONTINUE RETURN END