LoginSignup
3

More than 3 years have passed since last update.

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

Last updated at Posted at 2019-05-26

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

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

追記

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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
3