使用Fortran求一个矩阵的逆矩阵

下面是一个使用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分解和高斯消元方法计算矩阵的逆矩阵