Fortran 95 source code to calculate an inverse matrix by cofactor matrix.

inverse_matrix.f95

!****************************************************
! 逆行列の計算（余因子行列）
!
!   DATE          AUTHOR          VERSION
!   2019.12.24    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2013 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const

module comp
  use const
  implicit none
  private
  public :: det, cof_mtx, inv_mtx

contains
  ! 行列式計算
  ! * 再帰的
  !
  ! :param(in) real(8) m(:, :): 行列
  ! :return    real(8)       d: 行列式
  recursive real(DP) function det(m)
    implicit none
    real(DP), intent(in) :: m(:, :)
    integer(SP) :: j, n

    det = 0.0_DP
    n = size(m(1, :))
    select case(n)
    case (1)
      det = m(1, 1)
    case (2)
      det = m(1, 1) * m(2, 2) &
        & - m(1, 2) * m(2, 1)
    case (3)
      det = m(1, 1) * m(2, 2) * m(3, 3) &
        & + m(1, 2) * m(2, 3) * m(3, 1) &
        & + m(1, 3) * m(2, 1) * m(3, 2) &
        & - m(1, 1) * m(2, 3) * m(3, 2) &
        & - m(1, 2) * m(2, 1) * m(3, 3) &
        & - m(1, 3) * m(2, 2) * m(3, 1)
    case (4:)
      det = 0
      do j = 1, n
        det = det + m(1, j) * cof(m, 1, j)
      end do
    end select
  end function det

  ! 余因子行列
  !
  ! :param(in)   m: 行列配列
  ! :param(out) cm: 余因子行列配列
  subroutine cof_mtx(m, cm)
    implicit none
    real(DP), intent(in)  ::  m(:, :)
    real(DP), intent(out) :: cm(:, :)
    integer(SP) :: i, j, n

    n = size(m(1, :))
    do i = 1, n
      do j = 1, n
        cm(i, j) = cof(m, j, i)
      end do
    end do
  end subroutine cof_mtx

  ! 逆行列
  !
  ! :param(in)   m: 行列配列
  ! :param(out) im: 逆行列配列
  subroutine inv_mtx(m, im)
    implicit none
    real(DP), intent(in)  ::  m(:, :)
    real(DP), intent(out) :: im(:, :)
    integer(SP) :: n
    real(DP), allocatable :: cm(:, :)

    n = size(m(1, :))
    if (n == 1) then
      im = 1.0_DP / m
      return
    end if
    allocate(cm(n, n))
    call cof_mtx(m, cm)
    im = cm / det(m)
    deallocate(cm)
  end subroutine inv_mtx

  ! 余因子（正方行列 A の (i, j) に対する）
  !
  ! :param(in)   m: 行列配列
  ! :param(in)   i: 行インデックス
  ! :param(in)   j: 列インデックス
  ! :param(out)  c: 余因子
  real(DP) function cof(m, i, j)
    implicit none
    real(DP),    intent(in)  :: m(:, :)
    integer(SP), intent(in)  :: i, j
    integer(SP) :: n
    real(DP), allocatable :: m2(:, :)

    n = size(m(1, :))
    allocate(m2(n - 1, n - 1))
    m2(1:i-1, 1:j-1) = m(1:i-1, 1:j-1)
    m2(1:i-1, j:n  ) = m(1:i-1, j+1:n)
    m2(i:n  , 1:j-1) = m(i+1:n, 1:j-1)
    m2(i:n  , j:n  ) = m(i+1:n, j+1:n)
    cof = (-1)**(i+j) * det(m2)
    deallocate(m2)
  end function cof
end module comp

program inverse_matrix
  use const
  use comp
  implicit none
  integer(SP) :: n, i
  real(DP)    :: d
  real(DP), allocatable :: m(:, :), cm(:, :), im(:, :)

  ! データ数読み込み
  read (*, *) n

  ! 配列用メモリ確保
  allocate(m(n, n))
  allocate(cm(n, n))
  allocate(im(n, n))

  ! データ読み込み
  do i = 1, n
    read (*, *) m(i, :)
  end do
  print *, "Matrix A ="
  do i = 1, n
    print *, m(i, :)
  end do

  ! 行列式計算
  d = det(m)
  print *, "det(A) =", d

  ! 余因子行列計算
  call cof_mtx(m, cm)
  print *, "Cofactor Matrix of A ="
  do i = 1, n
    print *, cm(i, :)
  end do

  ! 逆行列計算
  call inv_mtx(m, im)
  print *, "Inverse Matrix of A ="
  do i = 1, n
    print *, im(i, :)
  end do

  ! 配列用メモリ解放
  deallocate(m)
  deallocate(cm)
  deallocate(im)
end program inverse_matrix