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
```
となります。

# サンプルコード

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


```fortran: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
```

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

```fortran
call syev(Hs,evec,'N')
```
とすれば、対称行列の固有値が求まり、

```fortran
call syev(Hs,evec,'V')
```
とすれば固有ベクトルも求まります。簡単ですね。

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

```fortran
use lapack95
```
としています。

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