LoginSignup
5
11

More than 3 years have passed since last update.

@Baruim

BLAS/LAPACKのすすめ (量子物理学応用の観点から)

Last updated at Posted at 2020-05-25

$\phantom{}$$
\newcommand\nr{\notag\\}%タグなし改行用
\newcommand\ret{\notag\\&\qquad}%数式改行用
\newcommand\I{\mathrm{i}}
\newcommand\D{\mathrm{d}}
\newcommand\E{\mathrm{e}}
\newcommand\hc{\mathrm{h.c.}}
\newcommand\cc{\mathrm{c.c.}}
\newcommand\O[1]{\mathscr{O}\left(#1\right)}
\newcommand\abs[1]{{\left\rvert #1 \right\lvert}}
\newcommand\Res{\mathop{\mathrm{Res}}}
\newcommand\bra[1]{\mathinner{\langle{#1}|}}
\newcommand\ket[1]{\mathinner{|{#1}\rangle}}
\newcommand\braket[1]{\mathinner{\langle{#1}\rangle}}
\newcommand\Bra[1]{\left\langle#1\right|}
\newcommand\Ket[1]{\left|#1\right\rangle}
\newcommand\Braket[1]{\left\langle #1 \right\rangle}
\newcommand|{\middle|}%\Braket用
\DeclareMathOperator{\Log}{Log}
\DeclareMathOperator{\Arg}{Arg}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\Tr}{Tr}
\newcommand\dashint{\mathchoice
{\int\kern-10pt-}
{\int\kern-8.5pt-}
{\int\kern-6.1pt-}
{\int\kern-4.58pt-}}
\newcommand\for[1]{\quad\mathrm{for}\quad #1}
\newcommand\set[1]{\left\{#1\right\}}
$　量子力学や量子統計力学に基づいた理論計算は線型代数計算に帰着される場合が多い．
このような計算をFortranを利用して実行する場合は，BLASやLAPACKが大変便利である．
77の時代とは違い，最近のFortranでは呼び出しが非常に簡単で使い勝手もよい．
ここでは，やや具体性に欠くが，物理量を簡単に記述できることを示すための応用例を挙げる．
(というよりも一般論であり，Hilbert空間が有限の問題の全てがこの方法で解ける)

なお，Intel Math Kernel Library (MKL) と Intel Fortran の組み合わせを想定したが，
もちろんその他の実装やコンパイラでも同様な記述が可能である．

BLAS/LAPACK

BLAS(Basic Linear Algebra Subprograms): その名の通り，線型代数の基本的な計算コードが用意されている．
Level 1, Level 2, および Level 3 にカテゴライズされている．
リファレンスは重いため，先にこのページを見てからあたるとよい (Intel MKLのリファレンスマニュアル)．

LAPACK(Linear Algebra PACKage): 線型代数における種々の問題の数値解法が用意されている．

コンパイラ＆リンカオプションは，このページに自身の計算機環境の情報を入れると
以下のようにオプションが出てくるので，そのままコピペして使用する．

上記のように表示された場合，以下のようにしてビルド＆実行する．

$ ifort main.f90 -c  -I${MKLROOT}/include/intel64/lp64 -I${MKLROOT}/include
$ ifort main2.o  ${MKLROOT}/lib/intel64/libmkl_blas95_lp64.a ${MKLROOT}/lib/intel64/libmkl_lapack95_lp64.a -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_lp64.a ${MKLROOT}/lib/intel64/libmkl_sequential.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl
$ ./a.out

長いのでmakefileを書いた方が便利である．
以下はmain.f90をビルドする場合の最小構造である(タブに注意)．

makefile
COMPILER=ifort
LINK_OPTIONS=${MKLROOT}/lib/intel64/libmkl_blas95_lp64.a ${MKLROOT}/lib/intel64/libmkl_lapack95_lp64.a -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_lp64.a ${MKLROOT}/lib/intel64/libmkl_sequential.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl
COMPILE_OPTIONS=-I${MKLROOT}/include/intel64/lp64 -I${MKLROOT}/include

main.out: main.o
    $(COMPILER) -o $@ $^ $(LINK_OPTIONS)

%.o: %.f90
    $(COMPILER) -o $@ -c $< $(COMPILE_OPTIONS)

線型代数演算

内積と行列同士の積をとりあげる．

ベクトル同士の積

\begin{align}
z=\sum_i x_i y_i
\end{align}

の形の計算はdot,dotu,dotcで行える．例えば，

\begin{align}
\mathrm{dot}(\vec x,\vec y) &= \sum_i x_i y_i,\quad (実ベクトル同士)\\
\mathrm{dotu}(\vec x,\vec y) &= \sum_i x_i y_i,\quad (複素ベクトル同士)\\
\mathrm{dotc}(\vec x,\vec y) &= \sum_i x_i{}^* y_i,\quad (複素ベクトル同士)
\end{align}

のように使える．

行列同士の積

\begin{align}
C_{ij}=\sum_{k} A_{ik} B_{kj}
\end{align}

の形の計算はgemmで行える．例えば，

\begin{align}
C&=AB,\\
C&=A^\dagger B,\\
C&=AB^\dagger,
\end{align}

のような計算に使う．
使い方についてはコードを直接見たほうがわかりやすいであろう：

  use BLAS95
  integer,parameter :: wp = selected_real_kind(p=15) !倍精度
  integer,parameter :: dim = 100 ! 次元

  real(wp),dimension(dim)::u,v,w
  complex(wp),dimension(dim)::x,y,z
  complex(wp),dimension(dim,dim)::A, B, C

  w = dot(u,v)  ! w = u(1)v(1) + u(2)v(2) + … + u(dim)v(dim)
  z = dotu(x,y) ! z = x(1) y(1) + z(2) y(2) + … + x(dim) y(dim)
  z = dotc(x,y) ! z = conjg(x(1)) y(1) + conjg(x(2)) y(2) + … + conjg(x(dim)) y(dim)


  call gemm(A,B,C) ! C = A B

  call gemm(A,B,C,transa='T') ! C = A^T B, T:転置
  call gemm(A,B,C,transb='T') ! C = A B^T, T:転置

  call gemm(A,B,C,transa='C') ! C = A^† B, †:Hermite共役 (転置＋複素共役)
  call gemm(A,B,C,transb='C') ! C = A B^†, †:Hermite共役 (転置＋複素共役)

固有値問題

　問題によって様々な特殊ルーチンに分けられているが，
ここでは次の2つを取り上げる．

geev：一般複素行列用
heev：Hermite行列用（量子論で重要）

引数の仕様が微妙に違うので注意．

MKLのリファレンスマニュアルの一部を抜粋しよう：

jobはjobzの誤植？

こんな感じで書かれているので，確認してから使う．
上の画像では削ったが，入力パラメータと出力パラメータの説明も必読である．
Fortran77の記述は過去の遺産に触れる際に必要になるであろう．
我々はもちろんFortran95のインターフェイスを用いる．

一般複素行列の固有値と固有ベクトル

例えば3x3の複素行列が配列complex::M(3,3)に格納されているとき，

temp = M
call geev(temp,A) ! tempは破壊される

とすれば，complex::A(3)に固有値が3個格納される．
固有ベクトルも欲しいときは

temp = M
call geev(temp,A,vr=V) ! tempは破壊される

と書く．そうすると

\sum_{j=1}^3 M(i,j) V(j,k) = A(k) V(i,k),

を満たした複素数値がcomplex::V(3,3)とcomplex::A(3)に格納される．
左固有ベクトルが欲しいときは

temp = M
call geev(temp,A,vl=V) ! tempは破壊される

と書く．そうすると

\sum_{j=1}^3 V(j,k)^* M(j,i) = A(k) V(i,k)^*,

を満たした複素数値がcomplex::V(3,3)とcomplex::A(3)に格納される．
両方欲しいときは

temp = M
call geev(temp,A,vl=V1, vr=V2) ! tempは破壊される

と書く．

Hermite行列の固有値と固有ベクトル

例えば3x3のHermite行列が配列complex::M(3,3)に格納されているとき，

temp = M
call heev(temp,A) ! tempは破壊される

とすれば，real::A(3)に固有値が3個格納される．
固有ベクトルも欲しいときは

V = M 
call heev(V,A,jobz='V') ! VにMの固有ベクトルが上書きされる

と書く．そうすると，

\sum_{j=1}^3 M(i,j) V(j,k) = A(k) V(i,k),

を満たした実数値がreal::A(3)に，複素数値がcomplex::V(3,3)に格納される．
Hermite行列の固有値は実数であるから，vl=やvr=といったオプションはない．
すなわち，左固有ベクトルは上式両辺の複素共役を取ることにより，

\sum_{j=1}^3 V(j,k)^* M(j,i)  = A(k) V(i,k)^*,

となって求められるからである (ここで，$M(i,j)^* =M(j,i)$ と $A(k)^*=A(k)$ の関係に注意)．

逆行列

数値計算界隈では嫌われ者の逆行列であるが，Green関数法では主役である．
例えば，リゾルベントは複素数$z$とハミルトニアン$\mathcal{H}$から作られる行列 $z-\mathcal{H}$ の逆行列である．
逆行列もLAPACKで簡単に求めることができる．

以下の (純粋) 関数inv(mat)は，実行列matの逆行列を返す関数である．
LAPACKによるLU分解を経由して求めている．
なお，realを全てcomplexに書き換えれば複素行列版となる．
(moduleとinterfaceを使って総称版を作ると使い勝手がよくなる)．

    pure function inv(mat) ! pureをつけるとdo concurrentブロック中にも書ける．
        integer,parameter :: wp=selected_real_kind(p=15) !倍精度
        real(wp),intent(in)::mat(:,:)
        real(wp)::inv(size(mat,1),size(mat,1))

        integer::ipiv(size(mat,1))

        inv = mat
        call getrf(inv,ipiv) ! LAPACK95
        call getri(inv,ipiv) ! LAPACK95
    end function

完全なテストプログラムは以下のようである：

program example
    use LAPACK95
    implicit none
    integer,parameter :: wp=selected_real_kind(p=15) !倍精度

    real(wp),dimension(2,2) :: A =[1.0_wp,2.0_wp,3.0_wp,4.0_wp], B

    B = inv(A)
    write(6,'(100F10.5)') B(1,:)
    write(6,'(100F10.5)') B(2,:)

    contains
        pure function inv(mat) ! pureをつけるとdo concurrentブロック中にも書ける．
            real(wp),intent(in)::mat(:,:)
            real(wp)::inv(size(mat,1),size(mat,1))

            integer::ipiv(size(mat,1))

            inv = mat
            call getrf(inv,ipiv) ! LAPACK95
            call getri(inv,ipiv) ! LAPACK95
        end function
end program

これを実行すると次の結果が返ってくる

$ ./main.out
  -2.00000   1.50000
   1.00000  -0.50000

つまり，

\begin{align}
B:=
\begin{pmatrix}
1 & 3 \\
2 & 4
\end{pmatrix}^{-1}
=
\begin{pmatrix}
-2 & 3/2 \\
1 & -1/2
\end{pmatrix}
\end{align}

が正しく計算されている．

応用例

　最初述べた通り，量子力学や量子統計力学での計算は線型代数計算に帰着される場合が多い．
以下の3つのステップで有限温度における物理量の統計平均を計算することを考えてみよう．

1. Schödinger方程式を解く

　ハミルトニアンを$\mathcal{H}$とするとき，
次のSchödinger方程式：

\mathcal{H}\ket{\phi_k}=E_k\ket{\phi_k},

をあるHilbert空間内部 (次元を$N$とする) で解き，
固有ケット $\ket{\phi_k}$ とそれが属する固有エネルギ $E_k$ をLAPACKによって数値的に求めることを考えよう．
量子数は $k=1,2,\ldots,N$ のように通し番号とする．
$\mathcal{H}$は抽象的な演算子であり，$\ket{\phi_k}$は抽象的なベクトルであるから，このままでは計算機に入力できない．
まずは以下のように，具体的な数値として表現する必要がある．

　同じHilbert空間を張る別の正規直交完全基底を$\set{\ket{i}}$とする．
正規直交完全基底とは次式を満たすようなものである：

\begin{align}
\braket{i|j}&=\delta_{ij}, (正規直交性)\\
\sum_{i=1}^N\ket{i}\bra{i}&=1.　(完全性)
\end{align}

$\set{\ket{i}}$は我々がよく知っている基底とする1．
求める固有ケットを完備関係式 (上式の第2式) を用いて展開する：

\ket{\phi_k}=\sum_{j=1}^N \ket{j}\braket{j|\phi_k}.

ブラケットが完成している部分$\braket{j|\phi_k}$は複素数である2．
上式をSchödinger方程式に代入し，更にSchödinger方程式の両辺に左から$\bra{i}$を作用させると

\sum_{j=1}^N \braket{i|\mathcal{H}|j}\braket{j|\phi_k}=E_k\braket{i|\phi_k},

となる．これはちょうどLAPACKによって解かれる行列の固有値問題の形になっている．
つまり，配列として

M(i,j) \Leftarrow \braket{i|\mathcal{H}|j},

と与え，geev(M,a,vr=V)をcallすれば

\begin{align}
E_k&=a(k),\\
\ket{\phi_k}&=\sum_{i=1}^N\ket{i}V(i,k),
\end{align}

を得ることができる．
このように，簡単にSchödinger方程式の完全な解を得ることができる．
なお，$\set{\braket{i|\mathcal{H}|j}}$はHermite行列なので，geevでなく
(V=Mの後に) heev(V,a,jobz='V')をcallしても良い．

2. Schödinger方程式の固有ケットに関する行列を作る

　次にSchödinger方程式の固有ケットを使って，ある物理量$\mathcal{A}$の行列を求めてみよう．
すなわち，その行列要素：

A_\mathrm{eigen}(k,l):=\braket{\phi_k|\mathcal{A}|\phi_l},

を知りたいというわけである．
上式の対角成分は固有状態における$\mathcal{A}$の期待値となっており，それ自体が有用な物理量となっている．

さて，我々が予め具体的に数値で用意できる行列要素は $\set{\ket{i}}$ 基底での行列要素：

A(i,j) \Leftarrow \braket{i|\mathcal{A}|j},

であるから，基底を変換しなければならない．
基底の変換はブラケット記法ではまったく簡単な作業である．
すなわち，完備関係式を2つ挿入するだけである：

\braket{\phi_k|\mathcal{A}|\phi_l}
=
\sum_{i=1}^N\sum_{j=1}^N
\braket{\phi_k|i}
\braket{i|\mathcal{A}|j}
\braket{j|\phi_l}.

また，$\braket{\phi_k|i}=\braket{i|\phi_k}^*=V(i,k)^*$であるから，結果として

A_\mathrm{eigen}(k,l)
=
\sum_{i=1}^N\sum_{j=1}^N
V(i,k)^*
A(i,j)
V(j,l),

で計算できる．これはあからさまに行列の積の計算であるから

A_\mathrm{eigen}=V^\dagger A V,

によって全ての行列要素が一度に求められる (ユニタリ変換)．
つまり，次の2行で済む：

call gemm(V,A,temp,transa='C') ! A:given, temp:中間変数
call gemm(temp,V,Aeigen)

3. 有限温度での統計平均を計算する

　更に，逆温度$\beta$の熱浴に$\mathcal{H}$で記述される系が接触し，熱平衡状態を達成しているものとする．
このときの物理量の演算子 (Hermite演算子である) $\mathcal{A}$の統計平均：

\braket{\mathcal{A}}_\beta:=\Tr\rho(\beta) \mathcal{A},

を計算しよう．ここで統計演算子は次式で定義されている：

\rho(\beta):=\frac{\E^{-\beta \mathcal{H}}}{Z(\beta)},\quad Z(\beta):=\Tr \E^{-\beta\mathcal{H}}.

トレースは好きな基底でとることができるが，
heevによって$\mathcal{H}$の固有値問題は解かれているので，
$\ket{\phi_k}$で作られた行列のトレースをとろう：

\begin{align}
\braket{\mathcal{A}}_\beta
&=
\sum_{k=1}^N
\braket{\phi_k|\rho(\beta) \mathcal{A}|\phi_k}
=
\frac{1}{
\sum_{k'=1}^N
\E^{-\beta E_{k'}}
}
\sum_{k=1}^N
\E^{-\beta E_k}\braket{\phi_k|\mathcal{A}|\phi_k}.
\end{align}

$\braket{\phi_k|\mathcal{A}|\phi_k}$はすでに求めてあるのでdotを使って和を取るだけでよい．
ここで一般的な注意であるが，$E_k$は最小値が0以上である保証はないので，
次のようにエネルギをシフトさせておくのが常套手段である：

\begin{align}
\braket{\mathcal{A}}_\beta
=
\frac{1}{
\sum_{k'=1}^N
\E^{-\beta (E_{k'}-E_\mathrm{min})}
}
\sum_{k=1}^N
\E^{-\beta (E_k-E_\mathrm{min})}\braket{\phi_k|\mathcal{A}|\phi_k},
\end{align}

ここで$E_\mathrm{min}$は$\set{E_n}$の最小値 (基底エネルギ) である．
数学的には分母分子に$\E^{\beta E_\mathrm{min}}$を掛けただけだが，
これによって指数関数が低温領域の計算 ($\beta\to\infty$) でオーバフローすることを回避できる．

まとめると以下のようになる3．
このように，量子力学での計算のほとんどが行列計算であるため，BLAS/LAPACKの恩恵を多く享受できる．

program example
    use BLAS95
    use LAPACK95
    implicit none
    integer,parameter :: wp = selected_real_kind(p=15) !倍精度
    integer,parameter :: dim = 100 ! Hilbert空間の次元

    complex(wp),dimension(dim,dim)::Hamiltonian, Eigenvectors, A, Aeigen
    real(wp),dimension(dim)::Eigenvalues

    real(wp),dimension(dim)::BoltzmannFactors, StatisticalOperator
    real(wp)::beta, PartitionFunction, Aav

    ! ここに Hamiltonian と A を与えるコードを書く
    ! ここでは適当に乱数を入れている
    block
        real(wp),dimension(dim,dim)::x,y
        call random_number(x)
        call random_number(y)
        Hamiltonian = x + (0.0_wp,1.0_wp)*y
        Hamiltonian = Hamiltonian + transpose(conjg(Hamiltonian)) ! Hermite行列にする

        call random_number(x)
        call random_number(y)
        A = x + (0.0_wp,1.0_wp)*y
        A = A + transpose(conjg(A)) ! Hermite行列にする

        beta = 0.01_wp ! 計算する逆温度を設定
    end block

    block
        complex(wp)::temp(dim,dim)
        real(wp)::Aexp(dim)

        Eigenvectors = Hamiltonian
        call heev(Eigenvectors,Eigenvalues,jobz='V') ! Schödinger方程式を解く LAPACK95

        call gemm(Eigenvectors,A,temp,transa='C') ! temp = v^† A を計算 BLAS(Lv.3)
        call gemm(temp,Eigenvectors,Aeigen) ! Aeigen = temp v = v^† A v を計算 BLAS(Lv.3)

        Eigenvalues = Eigenvalues - minval(Eigenvalues) ! 基底エネルギを0にシフト
        BoltzmannFactors = exp(-beta * Eigenvalues) ! Boltzmann因子を計算
        PartitionFunction = sum( BoltzmannFactors ) ! 分配関数 Z(β)を計算 BLAS (Lv.1) の asum で置き換えても良い
        StatisticalOperator = BoltzmannFactors/PartitionFunction ! 統計演算子の行列ρ(β)を計算
        forall (integer::i=1:dim)  Aexp(i) = real(Aeigen(i,i)) ! 対角成分 (Aの期待値) から配列作成　…　diagみたいな関数はない？
        Aav = dot(StatisticalOperator, Aexp) ! 統計平均を計算 BLAS(Lv.1)
    end block

  end program

  1. いろいろな演算子の行列要素を具体的に数値で表現できるもの．数値計算の出発点となる基底．例えば，位置固有ケット基底や，そのFourier変換である平面波基底など． ↩

  2. この式は見慣れた$N$次元のベクトルの展開形：$\vec A = c_1 \vec a_1 + c_2 \vec a_2+\cdots c_N\vec a_N$と全く同様である．$\vec A$や$\vec a_i$は抽象ベクトルであるが，展開係数$c_i$は数値である (つまり$c_i$は計算機で計算できる)． ↩

  3. 関数型プログラミングを心がけたがforallのところはうまく書けなかった． ↩

5
11
0

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
Sign upLogin
5
11