LAPACKとは

「Liner Algebra PACKage」の略で、プロが作った線形代数のサブルーチン集（ライブラリ）。 線形代数のライブラリとしては事実上の標準。 様々なルーチンが含まれるが、「行列の固有値問題を解く（対角化）」、「線形方程式を解く（逆行列）」、 の２つの用途だけ知っておけばよいだろう。

古代の遺物、先人の遺産。

このLAPACK、かなり敷居が高い。 予備知識の無い初心者であれば、３日以内に使えるようになれば上出来、と考えるべきだろう。 正しく使うためにはそれなりの理解が必要だが、日本語で書かれたわかりやすい説明がほとんど存在せず、 また中身を直接理解しようとしても素人にはほぼ解読不可能である。 また古いFORTRAN77のコードで汎用性が高くなるよう書かれているため、引数が無駄に多くて使いづらい。
しかしながら、プロと呼ばれる人々が作っているため、速度と精度は折り紙付き。 また世界中の研究者が実際に使っているので、バグが残っている可能性はほぼ無いと考えて良い。 中身が見えないので最初は不安になるが、信頼性は非常に高い。

LAPACKを手に入れる

LAPACKに含まれるルーチンのソースコードはフリーで公開されている。
サブルーチン間に複雑な依存関係があるため、 生のソースコードを個別に持ってきてコンパイルする、という使い方は（普通は）せず、誰かがまとめて「ライブラリ」にしたものを利用する。 有料のFortran Compilerを購入すると、最初からライブラリとなったもの（例えば「mkl(Math Kernel Library)」）が付属する。

具体的な使い方

１．インターネットなり本なりで、「必要なルーチンの名前」を調べる。 命名規則：
必ず、頭には「単精度実数、倍精度実数、単精度複素数、倍精度複素数」を表す「S,D,C,Z」という接頭語が付く。
行列を扱うルーチンの場合、続く2文字は「一般、対称、エルミート」=「GE、SY、HE」など行列の性質を表す。
残る文字がルーチンの役割を示す。「固有値問題」=「EV (EigenValue)」

２．名前が分かったら、そのルーチンの使い方、特に引数を調べる。
分かりにくい引数：
　LDA：配列Aが実際に定義されているサイズ。
　 　F77では「配列を大きめに定義し、一部だけ使う」というのが普通だったため。
　 　f90できちんとallocateしているなら、そのまま配列のサイズ『n』だと思って良い。
　WORK：内部処理に使う作業領域。入力でも出力でもない。
　 　F77では動的にメモリを確保できなかったため。
　 　型と配列のサイズの定義を間違えなければ、値は気にしなくて良い。

３．（中身ではなく、使い方を）十分に理解したら、自分のプログラムに組み込む。 使い方は普通の（裸の）サブルーチンと一緒。callするだけ。 外部手続き扱いのため整合性チェックは行われないので、引数の型を間違えないように。

４．コンパイル。「ライブラリのパスを切る」「ライブラリを指定する」必要がある。 環境によって違うため一概には言えない。先生か先輩に教わろう。
Macの場合、

 $ ifort -framework Accelerate main.f90

とすれば良い。
ワークステーション、もしくはMacでMKLを使う場合は

 $ ifort main.f90 -lpthread -lm -mkl=sequential

とすれば動くかも知れない。
（＊環境による。詳しくはこちらを参照。）

対角化ルーチン

 ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO )

エルミート行列を対角化する。物理で使うのは、大抵これ。 使い方は下の例を参照の事。

引数１（入） CHARACTER*1 :: JOBZ
　『'N'』：固有ベクトルを計算しない
　『'V'』：『A』に固有ベクトルが上書きされ、返される
引数２（入） CHARACTER*1 :: UPLO
　『'U'』：上三角行列を使う
　『'L'』：下三角行列を使う
引数３（入） INTEGER :: N
　nxn正方行列『A』のサイズ『n』
引数４（入出） COMPLEX*16 :: A(LDA,N)
　対角化したい（倍精度複素数型エルミートの）行列『A』
　JOBZ='V'のとき、固有ベクトル（を並べた配列）が上書きされる
　JOBZ='N'であっても、元のAは返って来ないので注意！
引数５（入） INTEGER :: LDA
　配列『A』の第一次元が実際に定義されているメモリ上の長さ。Fortran90ならば『N』と同じ。
引数６（出） REAL*8 :: W(N)
　固有値を並べたベクトルが出てくる
引数７（作業領域） COMPLEX*16 :: WORK(LWORK)
引数８（入） INTEGER :: LWORK
　内部的な作業領域のサイズ。最低でも2*N-1必要
引数９（作業領域） REAL*8 :: WORK(3*N-2)
引数１０（出） INTEGER :: INFO
　内部でエラーが発生した場合、『0』でない値が出てくる

固有値（と、対応する固有ベクトル）は、小さい方から順番にソートされる。

例１「倍精度複素数型エルミート行列の固有値問題」

固有値問題"***EV"系の例として、"ZHEEV"の例文。他のタイプもほぼ同様。

 use bacs
 implicit none
 integer :: i,j,n
 complex(kind=double), allocatable :: Amat(:,:), egVecs(:,:)
 real   (kind=double), allocatable :: egVals(:)
 !--- LAPACK private variable ---  !LAPACK用
 integer :: nzWork, info
 complex(kind=double), allocatable :: zWork(:)
 real   (kind=double), allocatable :: dWork(:)

    read *, n

    !--- allocation ---
    allocate (Amat(n,n),egVecs(n,n),egVals(n))
    !-- for 'ZHEEV' --  !LAPACK用
    nzWork = (64+1)*n
    allocate (zWork(nzWork),dWork(3*n-2))

    !--- define matrix A ---
    {Amat(i,j)=...} ! 行列を定義

    egVecs(:,:) = Amat(:,:)

    !--- use 'ZHEEV' ---
    call ZHEEV('V','U',n,egVecs,n,egVals,zWork,nzWork,dWork,info)
    if(info/=0) then ! エラーのシグナルが返ってきた場合
       print *, "### ERROR ###  process aborted "
       print *, " >>> ZHEEV : error code = ", info
       STOP   !実行を中断
    end if

    do i=1,n
      print *, egVals(i)
    end do

『Amat』を定義する部分を変えれば動く。それ以外はほぼパターンだと思って良く、自分なりの工夫は無駄なのでコピー＆ペーストをお勧めする。
この例では、行列『Amat』を上書きされないよう、『egVecs』にコピーしてから使っている。
『*work』の役割と値は我々の関知するところではない。が、型とサイズを間違えてはいけないので、変数の定義とallocationはしっかり写すこと。
内部でエラーが発生した場合のエラー処理も行っている（省略可）。

実行すると、

勉強：逆行列と連立方程式

式の上ではよく出てくる逆行列だが、逆行列そのものの値が必要になる場合はそれほど多くない。
例えば「y = A^-1 x」（Aの逆行列とxの積、yを求める）という式は、「A y = x」という連立方程式に帰着できる。
「逆行列を計算した後、行列積を計算する」よりも「連立方程式を解く」方が、アルゴリズム上、速度・精度の面で有利になる場合が多い。

逆行列ルーチン

 ZGETRF( M, N, A, LDA, IPIV, INFO )

行列Aを下処理(LU分解)する

 ZGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO )

下処理済の行列を使い、連立方程式 A X = B を解く

 ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )

下処理済の行列の逆行列を求める

まず『ZGETRF』で行列を下処理した後、『ZGETRS』もしくは『ZGETRI』を使う。 詳しい使い方はそれぞれのルーチンを参照の事。

例２「倍精度複素数型行列の逆行列」

"ZGETRI"の例文。他のタイプもほぼ同様。

 use bacs
 implicit none
 integer :: i,j,n
 complex(kind=double), allocatable :: Amat(:,:), Ainv(:,:)
 !--- LAPACK private variable ---  !-- LAPACK用
 integer                           :: nzWork, info
 integer,              allocatable :: iPiv(:)
 complex(kind=double), allocatable :: zWork(:)

    read *, n

    !--- allocation ---
    allocate (Amat(n,n),Ainv(n,n))
    !-- for 'ZGETRF,ZGETRI' --  !-- おまじない
    nzWork = 64*n
    allocate (iPiv(n),zWork(nzWork))

    !--- define matrix A ---
    {Amat(i,j)=...} !-- 行列を定義

    Ainv(:,:) = Amat(:,:) !-- 元のAを残しておく（ZGETRFで上書きされるため）

    !--- use 'ZGETRF' ---
    call ZGETRF(n,n,Ainv,n,iPiv,info)
    if(info/=0) then !-- エラーのシグナルが返ってきた場合
       print *, "### ERROR ###  process aborted "
       print *, " >>> ZGETRF : error code = ", info
       STOP   !-- 実行を中断
    end if
    !--- use 'ZGETRI' ---
    call ZGETRI(n,Ainv,n,iPiv,zWork,nzWork,info)
    if(info/=0) then !-- エラーのシグナルが返ってきた場合
       print *, "### ERROR ###  process aborted "
       print *, " >>> ZGETRI : error code = ", info
       STOP   !-- 実行を中断
    end if

    Amat = matmul(Amat,Ainv)
    do i=1,n
       print '(99f4.1)',  dble(Amat(i,1:n))
    end do

『Amat』を定義する部分を変えれば動くはず。
この例では、行列『Amat』を上書きされないよう、『Ainv』にコピーしてから使っている。最後は単位行列に戻るか確認しているだけ。
『iPiv』『nzWork』『zWork』の定義とallocationはしっかり写すこと。
エラー処理は省略しても良い。

実行すると、

能書き

LAPACKはとても奥が深い。 より簡単にプログラミングをするためには、ツールを自前で用意しておくと便利である。
以下に、私が普段使っているツールの一部を紹介する。参考にして自分で作るなり、私に教えを請うなりすると良い。

その１　固有値問題のwrapper: LAP_ev

LAPACKを直接呼ぶと、「おまじない」が多くて分かりづらい。
そこで、ややこしい部分はモジュールに押し込めてブラックボックス化してしまおう。
LAP_ev.f90 (v1.8.1, 限定解除)
これを使うと、次のようなプログラムになる。

program taikakuka
 use LAP_ev
 implicit none
 integer, parameter :: double=selected_real_kind(15)
 integer :: i,j,n,Lx
 real(kind=double), allocatable :: Hamiltonian(:,:),eigenStates(:,:)
 real(kind=double), allocatable :: eigenEnergies(:)
    print *, "length of system = ?"
    read *, Lx
    n = Lx  !--- size of matrix
    allocate (Hamiltonian(n,n),eigenStates(n,n),eigenEnergies(n))

    !--- define Hamiltonian in matrix form ---
    {Hamiltonian(i,j)=...} ! 行列を定義

    !--- solve eigenvalue problem ---
    call LAP__ev(n,Hamiltonian,eigenEnergies,eigenStates)

    do i=1,n
      print *, eigenEnergies(i)
    end do

    print *, "probability amplitude"
    j = n/2  !--- j番目の波動関数の確率振幅
    do i=1,n
      print *, i, abs(eigenStates(i,j))**2
    end do
end program

「LAP_ev」(_は１つ)というモジュールの
「LAP__ev」(__は２つ) というサブルーチンを使って
「n」次元の「行列」の「固有値」と「固有ベクトル」を求める、
という簡単な構造になる。皆さんはそれ以外の部分（ハミルトニアンの定義）とかに集中すれば良い。

コンパイルは次のようにする。

 $ ifort -framework Accelerate LAP_ev.f90 taikakuka.f90

その２　行列演算ライブラリ：matkob2

ちょっとした計算にいちいちLAPACKやBLASを呼び出して使うのは大げさである。
そこで、基本的な行列演算ルーチンは手作りしてまとめておくと便利である。
matkob2.f90 (v3.2.1, 手渡し限定)

 use matkob2

これを使うと例えば、「行列Aと、行列Bのエルミート共役の積の逆行列をCとする」が一行で書ける。

 C = matinv(matpro(A,matdag(B)))