编辑代码

       PROGRAM TestBRINV
        IMPLICIT NONE
        INTEGER, PARAMETER :: N=6
        DOUBLE PRECISION, DIMENSION(N, N) :: A
        DOUBLE PRECISION, DIMENSION(N, N) :: ExpectedB
        DOUBLE PRECISION, DIMENSION(N) :: IS, JS
        DOUBLE PRECISION :: PROPS(100)  ! 假设 PROPS 用于存储矩阵数据
        INTEGER :: i, j

        ! 初始化矩阵 A
        DATA A / &
     & 1.0D0, 2.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, &
     & 0.0D0, 1.0D0, 3.0D0, 0.0D0, 0.0D0, 0.0D0, &
     & 0.0D0, 0.0D0, 1.0D0, 4.0D0, 0.0D0, 0.0D0, &
     & 0.0D0, 0.0D0, 0.0D0, 1.0D0, 5.0D0, 0.0D0, &
     & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 1.0D0, 6.0D0, &
     & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 1.0D0 /

        ! 初始化预期的逆矩阵 B
        DATA ExpectedB / &
     & 1.0D0, -2.0D0, 6.0D0, -24.0D0, 120.0D0, -720.0D0, &
     & 0.0D0, 1.0D0, -3.0D0, 12.0D0, -60.0D0, 360.0D0, &
     & 0.0D0, 0.0D0, 1.0D0, -4.0D0, 20.0D0, -120.0D0, &
     & 0.0D0, 0.0D0, 0.0D0, 1.0D0, -5.0D0, 30.0D0, &
     & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 1.0D0, -6.0D0, &
     & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 1.0D0 /

        ! 将矩阵 A 的数据传递给 PROPS 数组
        ! 假设 PROPS 从索引 1 到 36 存储 A 的元素，按列优先顺序
        DO j = 1, N
          DO i = 1, N
            PROPS((j-1)*N + i) = A(i, j)
          END DO
        END DO

        ! 调用 BRINV 子程序对矩阵 A 进行求逆
        CALL BRINV(PROPS, N)

        ! 将求逆后的结果从 PROPS 数组恢复到矩阵 A
        DO j = 1, N
          DO i = 1, N
            A(i, j) = PROPS((j-1)*N + i)
          END DO
        END DO

        ! 比较求逆后的矩阵 A 与预期的逆矩阵 B
        PRINT *, '验证 BRINV 子程序的结果：'
        DO i = 1, N
          DO j = 1, N
            IF (ABS(A(i,j) - ExpectedB(i,j)) > 1.0E-6D0) THEN
              PRINT *, '不匹配位置 (', i, ',', j, '): A(i,j) = ', A(i,j), ' /= B(i,j) = ', ExpectedB(i,j)
            END IF
          END DO
        END DO

        PRINT *, '矩阵求逆验证完成。'

      END PROGRAM TestBRINV
*
      SUBROUTINE BRINV(A,N)
	DIMENSION A(N,N),IS(N),JS(N)

	B=1.0
	DO 100 K=1,N
	  D=0.0
	  DO 10 I=K,N
	  DO 10 J=K,N
	    IF (ABS(A(I,J)).GT.D) THEN
	      D=ABS(A(I,J))
	      IS(K)=I
	      JS(K)=J
	    END IF
10	  CONTINUE
	  IF (D+1.0.EQ.1.0) THEN
	    B=0.0
	    WRITE(*,20)
	    RETURN
	  END IF
20	  FORMAT(1X,'ERR**NOT INV')
	  DO 30 J=1,N
	    T=A(K,J)
	    A(K,J)=A(IS(K),J)
	    A(IS(K),J)=T
30	  CONTINUE
	  DO 40 I=1,N
	    T=A(I,K)
	    A(I,K)=A(I,JS(K))
	    A(I,JS(K))=T
40	  CONTINUE
	  A(K,K)=1/A(K,K)
	  DO 50 J=1,N
	    IF (J.NE.K) THEN
	      A(K,J)=A(K,J)*A(K,K)
	    END IF
50	  CONTINUE
	  DO 70 I=1,N
	    IF (I.NE.K) THEN
	      DO 60 J=1,N
	        IF (J.NE.K) THEN
	          A(I,J)=A(I,J)-A(I,K)*A(K,J)
	        END IF
60	      CONTINUE
	    END IF
70	  CONTINUE
	  DO 80 I=1,N
	    IF (I.NE.K) THEN
	      A(I,K)=-A(I,K)*A(K,K)
	    END IF
80	  CONTINUE
100	CONTINUE
	DO 130 K=N,1,-1
	  DO 110 J=1,N
	    T=A(K,J)
	    A(K,J)=A(JS(K),J)
	    A(JS(K),J)=T
110	  CONTINUE
	  DO 120 I=1,N
	    T=A(I,K)
	    A(I,K)=A(I,IS(K))
	    A(I,IS(K))=T
120	  CONTINUE
130	CONTINUE
	RETURN
      END      
end