!------------------------------------------------------------------------ ! Copyright (c) 2005-2019 SPMODEL Development Group. All rights reserved. !------------------------------------------------------------------------ ! !表題 eigmatrix テストプログラム ! !履歴 2005/01/26 竹広真一 ! 2007/11/02 竹広真一 メッセージ, 正誤判定追加 ! 2008/12/05 佐々木洋平. 文字列 sort の長さを 2 に ! 2011/02/18 佐々木洋平 dc_test を使用するように修正 ! 2019/10/10 佐々木洋平 モジュール名変更 ! 2022/10/27 竹広真一 倍精度用テスト ! program eigmatrix_test use dc_message, only : MessageNotify use dc_test, only : AssertEqual use dc_types, only: DP use spml_utils, only: eigen implicit none real(DP), dimension(:,:), allocatable :: amatrix real(DP), dimension(:), allocatable :: eigval_r, eigval_i real(DP), dimension(:,:), allocatable :: eigvec_r, eigvec_i real(DP), dimension(:), allocatable :: eigval_r_sol, eigval_i_sol real(DP), dimension(:,:), allocatable :: eigvec_r_sol, eigvec_i_sol complex(kind(0d0)), dimension(:),allocatable :: cwork real(DP),dimension(:), allocatable :: rwork integer :: info real(DP) :: error integer :: i ! 判定誤差設定 integer, parameter :: check_digits = 10 integer, parameter :: ignore = -11 call MessageNotify('M','eigen_test(double)','Test of eigmatrix') !--------------- 3x3 行列 ------------------ ! VALUE = 1, 2, 3 ! VECTOR = (-15,12,4), (-16,13,4), (-4,3,1) allocate(amatrix(3,3)) allocate(eigval_r(3), eigval_i(3)) allocate(eigvec_r(3,3), eigvec_i(3,3)) allocate(eigval_r_sol(3), eigval_i_sol(3)) allocate(eigvec_r_sol(3,3), eigvec_i_sol(3,3)) !---- 行列 ---- amatrix(:,1) = (/ 33.0d0, -24.0d0, -8.0d0 /) amatrix(:,2) = (/ 16.0d0, -10.0d0, -4.0d0 /) amatrix(:,3) = (/ 72.0d0, -57.0d0, -17.0d0 /) !---- 正解 ---- eigval_r_sol = (/3.0D0,2.0D0,1.0D0/) eigval_i_sol = (/0.0D0,0.0D0,0.0D0/) eigvec_r_sol(:,1) = (/-4.0D0,3.0D0,1.0D0/) eigvec_r_sol(:,2) = (/-16.0D0,13.0D0,4.0D0/) eigvec_r_sol(:,3) = (/-15.0D0,12.0D0,4.0D0/) eigvec_i_sol = 0.0D0 !---- 固有値計算 ---- call eigen(amatrix,eigval_r,eigval_i,eigvec_r,eigvec_i,info,& sort=' R',reverse=.true.) do i=1,3 eigvec_r(:,i) = eigvec_r(:,i)/eigvec_r(1,i) eigvec_r_sol(:,i) = eigvec_r_sol(:,i)/eigvec_r_sol(1,i) enddo !---- 判定 ---- call AssertEqual( & message= '3x3 arrray eigenvalue(real part)', & answer = eigval_r_sol, & check = eigval_r, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual( & message= '3x3 arrray eigenvalue(imaginary part)', & answer = eigval_i_sol, & check = eigval_i, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual( & message= '3x3 arrray eigenvector(real part)', & answer = eigvec_r_sol, & check = eigvec_r, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual( & message= '3x3 arrray eigenvector(imaginary part)', & answer = eigvec_i_sol, & check = eigvec_i, & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','Test of eigmatrix','Test of 3x3 array succeeded!') deallocate(amatrix) deallocate(eigval_r, eigval_i) deallocate(eigvec_r, eigvec_i) deallocate(eigval_r_sol, eigval_i_sol) deallocate(eigvec_r_sol, eigvec_i_sol) !--------------- 4x4 行列 ------------------ ! VALUE = 12, 1+5I, 1-5i, 2 ! VECTOR = (1,-1,1,1), (1,i,i,-1), (1,-i,-i,-1), (1,1,-1,1) allocate(amatrix(4,4)) allocate(eigval_r(4), eigval_i(4)) allocate(eigvec_r(4,4), eigvec_i(4,4)) allocate(eigval_r_sol(4), eigval_i_sol(4)) allocate(eigvec_r_sol(4,4), eigvec_i_sol(4,4)) allocate(cwork(4),rwork(4)) !---- 行列 ---- amatrix(:,1) = (/ 4.0d0, 0.0d0, 5.0d0, 3.0d0 /) amatrix(:,2) = (/ -5.0d0, 4.0d0, -3.0d0, 0.0d0 /) amatrix(:,3) = (/ 0.0d0, -3.0d0, 4.0d0, 5.0d0 /) amatrix(:,4) = (/ 3.0d0, -5.0d0, 0.0d0, 4.0d0 /) !---- 正解 ---- eigval_r_sol = (/ 1.0D0,1.0D0,2.0D0,12.0D0/) eigval_i_sol = (/-5.0D0,5.0D0,0.0D0, 0.0D0/) eigvec_r_sol(:,1) = (/1.0D0, 0.0D0, 0.0D0,-1.0D0/) eigvec_i_sol(:,1) = (/0.0D0,-1.0D0,-1.0D0, 0.0D0/) eigvec_r_sol(:,2) = (/1.0D0, 0.0D0, 0.0D0,-1.0D0/) eigvec_i_sol(:,2) = (/0.0D0, 1.0D0, 1.0D0, 0.0D0/) eigvec_r_sol(:,3) = (/1.0D0, 1.0D0,-1.0D0, 1.0D0/) eigvec_i_sol(:,3) = (/0.0D0, 0.0D0, 0.0D0, 0.0D0/) eigvec_r_sol(:,4) = (/1.0D0,-1.0D0, 1.0D0, 1.0D0/) eigvec_i_sol(:,4) = (/0.0D0, 0.0D0, 0.0D0, 0.0D0/) !---- 固有値計算 ---- call eigen(amatrix,eigval_r,eigval_i,eigvec_r,eigvec_i,info,& sort='RA') do i=4,1,-1 cwork = (eigvec_r(i,:)+(0,1)*eigvec_i(i,:))& /(eigvec_r(1,:)+(0,1)*eigvec_i(1,:)) eigvec_r(i,:) = real(cwork) eigvec_i(i,:) = aimag(cwork) enddo !---- 判定 ---- if (eigval_i(1) > 0.0 ) then error = eigval_r_sol(1) eigval_r_sol(1) = eigval_r_sol(2) eigval_r_sol(2) = error error = eigval_i_sol(1) eigval_i_sol(1) = eigval_i_sol(2) eigval_i_sol(2) = error do i=1,4 rwork = eigvec_r_sol(:,1) eigvec_r_sol(:,1) = eigvec_r_sol(:,2) eigvec_r_sol(:,2) = rwork rwork = eigvec_i_sol(:,1) eigvec_i_sol(:,1) = eigvec_i_sol(:,2) eigvec_i_sol(:,2) = rwork end do end if !---- 判定 ---- call AssertEqual( & message= '4x4 arrray eigenvalue(real part)', & answer = eigval_r_sol, & check = eigval_r, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual( & message= '4x4 arrray eigenvalue(imaginary part)', & answer = eigval_i_sol, & check = eigval_i, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual( & message= '4x4 arrray eigenvector(real part)', & answer = eigvec_r_sol, & check = eigvec_r, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual( & message= '4x4 arrray eigenvector(imaginary part)', & answer = eigvec_i_sol, & check = eigvec_i, & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','Test of eigmatrix(double)','Test of 4x4 array succeeded!') end program eigmatrix_test