FortranでLapackを使った行列の対角化

Fortranで数値計算をしなければならない場合があるかもしれません。
例えば、研究室で代々FortranやFORTRANを使っているとかですね。
あとは、スパコンを使って計算する場合でしょうか。
最近のFortranではオブジェクト指向もできますし、バグの少ない書きやすいコードも書けるようになってきました。
"Fortranからみたオブジェクト指向：オブジェクト指向でFortranコードを書く"
https://qiita.com/cometscome_phys/items/f87080286c6fdf72e49a

今回の記事では、一番基本である行列の対角化について書きます。
行列の対角化は、通常Lapackを使います。最近はMKLが無料で使えるようになりましたので、
https://software.intel.com/en-us/mkl/choose-download
FortranならMKLで行列を対角化するのが良いでしょう。

Lapackでのサブルーチンには引数がたくさんありますので、どのように指定すればいいか初心者には難しいかもしれません。
そこで、eigen(H,e)と呼べば固有値eが,eigen(H,e,V)と呼べば固有ベクトルが呼ばれるサブルーチンを作ってみましょう。
Fortran90のInterfaceを使えば、Hが実数か複素数かをうまいこと判別してくれるコードを書くことができます。

#バージョン
gfortran gcc version 8.2.0 (Homebrew GCC 8.2.0)

#コンパイルコマンド
MKLのコンパイルコマンドは状況によって異なりますが、

gfortran test.f90 -lmkl_intel_lp64  -lmkl_sequential -lmkl_core -lpthread -lm -ldl 

のような感じでコンパイルします。

#コード

コードは以下の通りです。

module linearalgebra
    implicit none
    interface eigen
        module procedure eigen_dble_sym_eigenvalues
        module procedure eigen_dble_sym_both
        module procedure eigen_complex_her_eigenvalues
        module procedure eigen_complex_her_both
    end interface eigen

    contains

    subroutine eigen_dble_sym_eigenvalues(H,e)
        implicit none
        real(8),intent(in)::H(:,:)
        real(8),intent(out)::e(:)

        real(8),allocatable::V(:,:),work(:)
        integer::n,lwork,info
        
        n = ubound(H,1)
        lwork = 3*n-1
        write(*,*) n
        allocate(V(1:n,1:n))
        allocate(work(lwork))
        V = H
        call dsyev('N', 'U', n, V, n, e, work, lwork, info)
        if (info .ne. 0) then
            write(*,*) "error! in eigen. info = ",info
        end if
        deallocate(V,work)

        return

    end subroutine eigen_dble_sym_eigenvalues


    subroutine eigen_dble_sym_both(H,e,V)
        implicit none
        real(8),intent(in)::H(:,:)
        real(8),intent(out)::e(:)
        real(8),intent(out)::V(:,:)

        real(8),allocatable::work(:)
        integer::n,lwork,info
        
    
        n = ubound(H,1)
        lwork = 3*n-1
        allocate(work(lwork))
        V = H
        call dsyev('V', 'U', n, V, n, e, work, lwork, info)
        if (info .ne. 0) then
            write(*,*) "error! in eigen. info = ",info
        end if
        deallocate(work)

        return

    end subroutine eigen_dble_sym_both  

    subroutine eigen_complex_her_eigenvalues(H,e)
        implicit none
        complex(8),intent(in)::H(:,:)
        real(8),intent(out)::e(:)

        complex(8),allocatable::V(:,:),work(:)
        integer::n,lwork,info
        integer,allocatable::rwork(:)
        
        n = ubound(H,1)
        lwork = 3*n-1
        allocate(V(1:n,1:n))
        allocate(work(lwork))
        allocate(rwork(2*n-2))
        V = H
        call zheev('N', 'U', n, V, n, e, work, lwork, rwork, info)
        if (info .ne. 0) then
            write(*,*) "error! in eigen. info = ",info
        end if

        deallocate(work,rwork,V)

        return

    end subroutine eigen_complex_her_eigenvalues  
    
    subroutine eigen_complex_her_both(H,e,V)
        implicit none
        complex(8),intent(in)::H(:,:)
        real(8),intent(out)::e(:)

        complex(8),intent(out)::V(:,:)
        complex(8),allocatable::work(:)
        integer::n,lwork,info
        integer,allocatable::rwork(:)

        n = ubound(H,1)
        lwork = 3*n-1
        allocate(work(lwork))
        allocate(rwork(2*n-2))
        V = H
        call zheev('V', 'U', n, V, n, e, work, lwork, rwork, info)
        if (info .ne. 0) then
            write(*,*) "error! in eigen. info = ",info
        end if

        deallocate(work,rwork)

        return

    end subroutine eigen_complex_her_both  

end module linearalgebra


subroutine test()
    use linearalgebra
    implicit none
    integer::n
    real(8),allocatable::Hs(:,:)
    real(8),allocatable::evec(:)
    real(8),allocatable::Vs(:,:)
    complex(8),allocatable::V(:,:)
    complex(8),allocatable::H(:,:)
    real(8),allocatable::Hr(:,:)
    real(8),allocatable::Hi(:,:)
    complex,parameter::ci = (0d0,1d0)
    n = 2
    allocate(Hs(1:n,1:n))
    
    allocate(V(1:n,1:n))
    allocate(Vs(1:n,1:n))
    allocate(evec(1:n))

    call random_number(Hs)
    
    Hs = (Hs + transpose(Hs))/2
    !write(*,*) Hs

    call eigen(Hs,evec)
    write(*,*) evec

    call eigen(Hs,evec,Vs)
    write(*,*) evec

    allocate(H(1:n,1:n))
    allocate(Hi(1:n,1:n))
    allocate(Hr(1:n,1:n))
    call random_number(Hr)
    call random_number(Hi)
    H = Hr + ci*Hi

    H = (H + conjg(transpose(H)))/2

    call eigen(H,evec)
    write(*,*) evec

    call eigen(H,evec,V)
    write(*,*) evec



    return
end subroutine


program main
    use linearalgebra
    implicit none
    call test()
end program

moduleのlinearalgebraをuseすると、call eigenで対角化ができるようになります。
ここで、

    interface eigen
        module procedure eigen_dble_sym_eigenvalues
        module procedure eigen_dble_sym_both
        module procedure eigen_complex_her_eigenvalues
        module procedure eigen_complex_her_both
    end interface eigen

としているのは、引数の数や型によって呼ぶサブルーチンを変えるためです。

追記

上記の固有値問題は、対称行列あるいはエルミート行列を対象としています。