home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
ftp.barnyard.co.uk
/
2015.02.ftp.barnyard.co.uk.tar
/
ftp.barnyard.co.uk
/
cpm
/
walnut-creek-CDROM
/
SIMTEL
/
CPMUG
/
CPMUG049.ARK
/
MINV.RAT
< prev
next >
Wrap
Text File
|
1984-04-29
|
2KB
|
97 lines
#Test driver for MINV
PROGRAM MINVT
INTEGER*1 L,N
DIMENSION A(5,5),L(5)
include F77DEF
DATA N/5/
include CONSTS
DO I=1,N
$(DO J=1,N
A(J,I)=1./REAL(I+J-1)$)#Set up arbitrary matrix
write(kons,20) A
20 format(1x,5G12.6)
call MINV(A,N,L,DET)
write(kons,20)DET
write(kons,20)A
#Now try to get original back
call MINV(A,N,L,DET)
WRITE(KONS,20)DET
WRITE(KONS,20)A
STOP
END
#Gauss-Jordan inversion & determinant calculation
SUBROUTINE MINV(A,N,L,DET)
DIMENSION A(N,N),L(N)
#Use FORTRAN (col,row) subscript convention
INTEGER*1 L,J,K,N,IT,I
include CONSTS
ND=N#Convert to *2 for calls
#Tag original column
DO J=1,N
L(J)=J
DET=1.
DO K=1,N;$(
#Choose pivot row for maximum value and find if 0.
#Start at diagonal and look downwards
IF(FMXMGV(A(K,K),N,N-K+1,JR)==0.)
$(WRITE(KONS,900);BREAK$)
900 FORMAT(1x,'Singular Matrix')
IF(JR^=1)
$(J=K-1+JR#Convert from J counting from diagonal
#to J subscript
#Swap rows K & J and tags
IT=L(K)
L(K)=L(J)
L(J)=IT
CALL FVSWAP(A(1,K),1,A(1,J),1,ND)
#Swap changes determinant sign
DET=-DET$)#End of swapping block
#Determinant calculation
Q=A(K,K)
DET=Q*DET
A(K,K)=1.#Start changing this part of I matrix to inverse
# Divide row by diagonal, including both modified
# original and developing inverse
DO I=1,N
A(I,K)=A(I,K)/Q
#Eliminate off-diagonal elements of col K
FOR(J=1;J<=N;J=J+1)$(IF(J==K)NEXT
Q=-A(K,J)
#Bring in off-diagonal element of I matrix
A(K,J)=0.
DO I=1,N;A(I,J)=A(I,K)*Q+A(I,J)$)#End FOR
$)#End DO K=1,N
#Undo column swaps
N1=N-1
DO K=1,N1;$(
#Search for K
FOR(I=K;I<=N&L(I)^=K;I=I+1);#No body of FOR
CALL FVSWAP(A(K,1),ND,A(I,1),ND,ND)
L(I)=L(K)$)#Keep track of columns yet to go
RETURN
END
FUNCTION FMXMGV(A,IA,L,K)
#Start at A, at intervals IA find largest of L numbers
#and return position in K
DIMENSION A(1)#Actually (IA,L)
FMXMGV=0.
IS=1
DO I=1,L
$(ABSA=ABS(A(IS))
IF(ABSA>FMXMGV)
$(FMXMGV=ABSA;K=I$)
IS=IS+IA$)
RETURN
END
SUBROUTINE FVSWAP(A,IA,B,IB,N)
DIMENSION A(1),B(1)#Actually (IA,N),IB,N)
#Swap REAL vector starting at A interval IA
#with B interval IB, length N
ISA=1
ISB=1
DO I=1,N
$(T=A(ISA);A(ISA)=B(ISB);B(ISB)=T
ISA=ISA+IA;ISB=ISB+IB$)
RETURN
END