[数值问题] 网上找了个求解矩阵的程序

      program main
      implicit none
      real(8)::a3(18,9),A3inv(9,18)
      integer::log
      a3=1.
      call BGINV(18,9,A3,A3inv, Log, 1.e-15)
      end

!!-----------------------------------------
!!         general inverse
!!-----------------------------------------
       SUBROUTINE BGINV(M,N,A,AA,L,EPS)
       double precision :: A(M,N),U(M,M),V(N,N),AA(N,M),EPS
       double precision,ALLOCATABLE :: S(:), E(:), WORK(:)
       KA=MAX(M,N)+1
       ALLOCATE(S(KA))
       ALLOCATE(E(KA))
       ALLOCATE(WORK(KA))
       CALL BMUAV(A,M,N,U,V,L,EPS,KA,S,E,WORK)
       IF (L.EQ.0) THEN
       K=1
10     IF (A(K,K).NE.0.0) THEN
         K=K+1
       IF (K.LE.MIN(M,N)) GOTO 10
       END IF
       K=K-1
       IF (K.NE.0) THEN
       DO 40 I=1,N
       DO 40 J=1,M
        AA(I,J)=0.0
        DO 30 II=1,K
30        AA(I,J)=AA(I,J)+V(II,I)*U(J,II)/A(II,II)
40      CONTINUE
       END IF
       END IF
       DEALLOCATE(S)
       DEALLOCATE(E)
       DEALLOCATE(WORK)
       RETURN
       END

!!-----------------------------------------
!!         singular decomposite
!!-----------------------------------------
       SUBROUTINE BMUAV(A,M,N,U,V,L,EPS,KA,S,E,WORK)
!略去
       end