inmatrix.for

Fortran 的矩阵求逆程序

	!求矩阵的逆矩阵！
	PARAMETER (n=7)
	dimension xa(n,n)
	dimension xb(n-1,n-1)
	dimension xn(n,n)
	open(11,file="matrix0.ye")
	do 10 i=1,n
	read(11,*)(xa(i,j),j=1,n)
10 	continue
	close(11)
	call BSDET(xa,N,DET)
	xadet=det
	!求伴随矩阵开始，
	do 40 i=1,n
	do 40 j=1,n
	do 20 ii=1,n-1
	do 20 jj=1,n-1
	if (jj<j)then
	  if  (ii<i) then 
	    xb(ii,jj)=xa(ii,jj)
	  else
	    xb(ii,jj)=xa(ii+1,jj)
	  endif
	else
	  if  (ii<i) then 
	    xb(ii,jj)=xa(ii,jj+1)
	  else
	    xb(ii,jj)=xa(ii+1,jj+1)
	  endif
	endif
20	continue
	!求伴随矩阵结束。
	call BSDET(xb,N-1,DET)
	xbdet=det
	if (int((i+j)/2)*2.eq.(i+j)) then
		xn(j,i)=xbdet/xadet
	else
		xn(j,i)=xbdet/xadet*(-1)
	endif
40	continue
	open(12,file="matrixn0.ye")
	write(12,*)"|A|=",xadet
      write(12,*)"逆矩阵为："
	do 50 i=1,n
	write(12,*)(xn(i,j),j=1,n)
50	continue

	!另法求逆矩阵。
	do 70 i=1,n
	do 70 j=1,n
	xn(i,j)=xa(i,j)
70	continue
	call inv(xn,n)

	write(12,*)"另|A|=",xadet
      write(12,*)"逆矩阵为："

	do 80 i=1,n
	write(12,*)(xn(i,j),j=1,n)
80	continue

	write(*,100)
100	format("计算结束！")
	close(12)
	end


!	高斯列主元法求行列式值！
	SUBROUTINE BSDET(B,N,DET)
	DIMENSION B(N,N)
	DIMENSION A(N,N)
!	DOUBLE PRECISION A,DET,F,D,Q
	do 5 i=1,n
	do 5 j=1,n
	a(i,j)=b(i,j)
5	continue
	F=1.0
	DET=1.0
	DO 100 K=1,N-1
	  Q=0.0
	  DO 10 I=K,N
	  DO 10 J=K,N
	    IF (ABS(A(I,J)).GT.Q) THEN
	      Q=ABS(A(I,J))
	      IS=I
	      JS=J
	    END IF
10	  CONTINUE
	  IF (Q+1.0.EQ.1.0) THEN
	    DET=0.0
	    RETURN
	  END IF
	  IF (IS.NE.K) THEN
	    F=-F
	    DO 20 J=K,N
	      D=A(K,J)
	      A(K,J)=A(IS,J)
	      A(IS,J)=D
20	    CONTINUE
	  END IF
	  IF (JS.NE.K) THEN
	    F=-F
	    DO 30 I=K,N
	      D=A(I,JS)
	      A(I,JS)=A(I,K)
	      A(I,K)=D
30	    CONTINUE
	  END IF
	  DET=DET*A(K,K)
	  DO 50 I=K+1,N
	    D=A(I,K)/A(K,K)
	    DO 40 J=K+1,N
40	    A(I,J)=A(I,J)-D*A(K,J)
50	  CONTINUE
100	CONTINUE
	DET=F*DET*A(N,N)
	RETURN
	END

        subroutine inv(a,n)
        dimension a(n,n)
        common/syspar/iunniti,iunito,idebug
c       check the debug information.
	  maxn=550
	  ier=1
        if(idebug.ne.0)Then
        Write(iunito,10)
10      Format(/,15x,'Debug information',/,15x,30('-'),/,15x,
     1  'Subroutine ffw24 called')
        Write(iunito,20)n,maxn
20      Format(/,15x,'   n = ',t30,i6,/,15x,'maxn = ',t30,i6,
     1  /,15x,30('-'))
        Endif
c       Carry out some spot checks mainly on the dimensions.
        if(n.gt.maxn)then
        ier = 999
        write(iunito,30)ier
30      Format(/,20x,'Error detected in ffw24 - ier = ',i3)
        go to 120
	  endif
        do 50 k=1,n
        if(a(k,k).eq.0.0) then
        ier=998
        write(iunito,40)k,k
40      format(/,15x,'Fatal Error detected in ffw24',
     1  /,15x,'because a(',i3,1h,,i3,') is a zero element.',/)
        go to 120
        endif
50      continue
        do 110 k=1,n
        do 70 j=1,n
        if(j-k) 60,70,60
60     a(k,j)=a(k,j)/a(k,k)
70     continue
        a(k,k)=1.0/a(k,k)
        do 110 i=1,n
        if(i-k) 80,110,80
80     continue
        do 100 j=1,n
        if(j-k) 90,100,90
90     a(i,j)=a(i,j)-a(i,k)*a(k,j)
100     continue
        a(i,k)=-a(i,k)*a(k,k)
110     continue
        ier=0
120     return
        end