--- /dev/null
+ 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