MKLのLapack95インターフェースを使ってFortranで行列を対角化しよう

Fortran Advent Calendar 2019です。

はじめに

FortranでLapackを使う時によくお世話になるMKLですが、マニュアルを見るとよく「Fortran95インターフェース」というのを見かけると思います。
これを見ると引数がとても減っていて、「なんかよさそう。でもどうやって使うんだろう？」と思う方も多いかと思います。
今回は、このインターフェースをgfortranで使う方法について述べます。
これを使えば、たくさんの引数を使わずに、PythonやJuliaのように簡単に対角化ができます。

以前
「FortranでLapackを使った行列の対角化」
https://qiita.com/cometscome_phys/items/e22d329697e378f6763d
という記事を書きましたが、この時は自前でインターフェースを作りました。
今回はMKLに入っているインターフェースを使います。
gfortranで使う場合には、
https://www.xlsoft.com/jp/products/intel/perflib/mkl/11.2/mkl_userguide_lnx/GUID-C3CFB9D7-EC03-4ECC-8F37-7DAAE2C6C93A.htm
にあるように、まずビルドする必要があります。
ビルドして使えるようになるまで、やってみましょう。

OS

OS: Ubuntu 16.04.6 LTS (Xenial Xerus)
gfortran: gcc version 7.4.0 (Ubuntu 7.4.0-1ubuntu1~16.04~ppa1)

Lapack95インターフェースのインストール

まず、Lapack95インターフェースをgfortran用にビルドします。
MKLは入っているとします。
つまり、

echo $MKLROOT

とした時に、適切にMKLのインストール先が入っているとしましょう。
まず、ホームディレクトリに適当なlapack95のインストール先のディレクトリを作成します。
例えば、

cd ~/
mkdir mkl

などですね。

次に、

cd $MKLROOT/interfaces/lapack95/
make libintel64  INSTALL_DIR=~/mkl/ interface=lp64 FC=gfortran

とすると、いま作ったディレクトリにgfortran用のlapack95インターフェースのライブラリやモジュールがビルドされます。
次に、.bashrcなどにF95ROOTという環境変数を設定します。インストール先のディレクトリの情報です。
例えば、

export F95ROOT=$HOME/mkl

ですね。

これでgfortranでも使えるようになりました。

コンパイル

コンパイル方法は
https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor
でご自分の環境に合わせて調べてみましょう。
Ubuntuであれば、

Use this link line: に

 ${F95ROOT}/lib/intel64/libmkl_lapack95_lp64.a -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl

Compiler options:

 -I${F95ROOT}/include/intel64/lp64 -m64 -I${MKLROOT}/include

とあると思います。
これを使います。

lapacktest.f90というコードがあった時、コンパイルは、

gfortran lapacktest.f90  ${F95ROOT}/lib/intel64/libmkl_lapack95_lp64.a -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl  -I${F95ROOT}/include/intel64/lp64 -m64 -I${MKLROOT}/include

となります。

サンプルコード

以下のようなサンプルコードを用意しました。

lapacktest.f90

subroutine test()
    use lapack95
    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

    call syev(Hs,evec,'N')
    write(*,*) evec

    call syev(Hs,evec,'V')
    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 heev(H,evec)
    write(*,*) evec

    call heev(H,evec,'V')
    write(*,*) evec



    return
end subroutine

program main
    implicit none
    call test()
end program

以前の記事と同じように、対称行列の固有値と複素数エルミート行列の固有値を求めています。

call syev(Hs,evec,'N')

とすれば、対称行列の固有値が求まり、

call syev(Hs,evec,'V')

とすれば固有ベクトルも求まります。簡単ですね。

lapack95インターフェースを使う時は、

use lapack95

としています。

これで、mklのマニュアルにあるいろんなlapack95インターフェースが使えるようになりました。
もう面倒な引数からはおさらばできますね。