Fortran 矩阵求逆：使用 LU 分解和高斯消元

下面是一个使用 Fortran 编写的程序，可以计算一个矩阵的逆矩阵：

program matrix_inverse
    implicit none
    integer, parameter :: n = 3 ! 矩阵的维度
    real :: matrix(n,n), inverse(n,n)
    integer :: ipiv(n), info
    integer :: i, j

    ! 输入矩阵
    print *, '输入矩阵：'
    do i = 1, n
        read *, matrix(i,:)
    end do

    ! 计算逆矩阵
    call sgetrf(n, n, matrix, n, ipiv, info)
    if (info /= 0) then
        print *, '无法计算逆矩阵'
        stop
    end if
    call sgetri(n, matrix, n, ipiv, inverse, n, info)
    if (info /= 0) then
        print *, '无法计算逆矩阵'
        stop
    end if

    ! 输出逆矩阵
    print *, '逆矩阵：'
    do i = 1, n
        print *, inverse(i,:)
    end do

contains

    ! LU分解
    subroutine sgetrf(m, n, a, lda, ipiv, info)
        implicit none
        integer, intent(in) :: m, n, lda
        real, intent(inout) :: a(lda,*)
        integer, intent(out) :: ipiv(n), info
        integer :: i, j, k, imax
        real :: tmp, maxval

        info = 0
        do j = 1, min(m,n)
            maxval = abs(a(j,j))
            imax = j
            do i = j+1, m
                if (abs(a(i,j)) > maxval) then
                    maxval = abs(a(i,j))
                    imax = i
                end if
            end do
            if (imax /= j) then
                do k = 1, n
                    tmp = a(j,k)
                    a(j,k) = a(imax,k)
                    a(imax,k) = tmp
                end do
            end if
            if (a(j,j) == 0.0) then
                info = j
                return
            end if
            ipiv(j) = j
            do i = j+1, m
                a(i,j) = a(i,j) / a(j,j)
                do k = j+1, n
                    a(i,k) = a(i,k) - a(i,j) * a(j,k)
                end do
            end do
        end do
    end subroutine sgetrf

    ! 计算逆矩阵
    subroutine sgetri(n, a, lda, ipiv, work, lwork, info)
        implicit none
        integer, intent(in) :: n, lda, lwork
        real, intent(inout) :: a(lda,*), work(lwork)
        integer, intent(in) :: ipiv(n)
        integer, intent(out) :: info
        integer :: i, j
        real :: wkopt

        call sgetri(n, a, lda, ipiv, & ! 输入参数
                    work, -1, info) ! 输出工作数组大小
        wkopt = work(1)
        lwork = max(1, int(wkopt))
        allocate(work(lwork))

        call sgetri(n, a, lda, ipiv, & ! 输入参数
                    work, lwork, info) ! 输出逆矩阵和信息

        deallocate(work)
    end subroutine sgetri

end program matrix_inverse

以上代码使用 LU 分解和高斯消元方法计算矩阵的逆矩阵。

LU 分解 是将矩阵分解为一个下三角矩阵 L 和一个上三角矩阵 U 的乘积，即 A = LU。

高斯消元 是一种常用的线性方程组求解方法，它通过对系数矩阵进行行变换将其转化为上三角矩阵，从而得到方程组的解。

在以上代码中，sgetrf 子程序实现 LU 分解，而 sgetri 子程序则使用 LU 分解和高斯消元法计算矩阵的逆矩阵。

代码说明：

• matrix: 输入矩阵
• inverse: 逆矩阵
• ipiv: 行交换指示符
• info: 错误信息
• sgetrf: LU 分解子程序
• sgetri: 逆矩阵计算子程序

本代码仅供参考，您可以根据自己的需求进行修改和完善。

注意：

• 该代码仅适用于方阵，即行数和列数相等的矩阵。
• 如果矩阵不可逆，则程序会输出错误信息并停止执行。