subroutine gauss(RO,AP,MT,M,N,*) c c CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION c RO IS A SQUARE MATRIX c THE CALCULATED PRODUCT IS STORED IN AP c ABNORMAL EXIT IF RO IS SINGULAR c integer MT, M, N, M1,I,J,IM, & I1,MI,MI1 double precision RO(MT,M),AP(MT,N),X,RM,PR, & Y if(M.ne.1)goto 10 X=RO(1,1) if(dabs(X).le.1.0D-13) return 1 X=1.0/X do 16 I=1,N 16 AP(1,I)=AP(1,I)*X return 10 continue M1=M-1 DO1 I=1,M1 IM=I RM=DABS(RO(I,I)) I1=I+1 do 2 J=I1,M if(DABS(RO(J,I)).LE.RM) goto 2 RM=DABS(RO(J,I)) IM=J 2 continue If(IM.eq.I)goto 17 do 3 J=1,N PR=AP(I,J) AP(I,J)=AP(IM,J) 3 AP(IM,J)=PR do 4 J=I,M PR=RO(I,J) RO(I,J)=RO(IM,J) 4 RO(IM,J)=PR 17 X=RO(I,I) if(dabs(X).le.1.0E-13) return 1 X=1.0/X do 5 J=1,N 5 AP(I,J)=X*AP(I,J) do 6 J=I1,M 6 RO(I,J)=X*RO(I,J) do 7 J=I1,M Y=RO(J,I) do 8 K=1,N 8 AP(J,K)=AP(J,K)-Y*AP(I,K) do 9 K=I1,M 9 RO(J,K)=RO(J,K)-Y*RO(I,K) 7 continue 1 continue X=RO(M,M) if(dabs(X).le.1.0E-13) return 1 X=1.0/X do 11 J=1,N 11 AP(M,J)=X*AP(M,J) do 12 I=1,M1 MI=M-I MI1=MI+1 do 14 J=1,N X=AP(MI,J) do 15 K=MI1,M 15 X=X-AP(K,J)*RO(MI,K) 14 AP(MI,J)=X 12 continue return end