Created
February 5, 2020 02:29
-
-
Save komasaru/3fcb6d5ded8393ae6557229f6d3f0511 to your computer and use it in GitHub Desktop.
Fortran 95 source code to calculate a determinant by cofactor expansion.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
!**************************************************** | |
! 行列式の計算(余因子展開) | |
! | |
! DATE AUTHOR VERSION | |
! 2019.12.23 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 | |
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 | |
! 余因子(正方行列 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 determinant | |
use const | |
use comp | |
implicit none | |
integer(SP) :: n, i | |
real(DP) :: d | |
real(DP), allocatable :: m(:, :) | |
! データ数読み込み | |
read (*, *) n | |
! 配列用メモリ確保 | |
allocate(m(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 | |
! 配列用メモリ解放 | |
deallocate(m) | |
end program determinant |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment