計算物理学２第 8 回：行列計算ライブラリver. 2018/6/13
連立一次方程式や固有値問題などの線形計算は物理学以外の分野でも現れる共通の問題であり、自ら計算コードを作らなくても使用できるライブラリがあります。これらは多くのユーザによって使われており、自分でコードを構築するよりも早く正確に計算することができます。対角化などの行列計算は行列の次元が大きくなると計算時間が増大するため、できる限り行列計算ライブラリを使用することをおすすめします。この資料では LAPACK ライブラリの使用方法について説明します。

1 LAPACKLAPACK(Linear Algebra PACKage, http://www.netlib.org/lapack/) は線形計算の Fortran パッケージです。C 言語から呼び出して使うことも可能です。LAPACK には様々なサブルーチンが含まれていますが、これらはドライバルーチン、計算ルーチン、補助ルーチンに分かれており、通常はドライバルーチンの中から目的にあったものを探して使います。ドライバルーチンには• 線形方程式問題• 線形最小二乗法問題・一般線形最小二乗法問題• 固有値問題、特異値問題・一般化固有値問題・一般化特異値問題などの問題に対するものが含まれています。

LAPACK のサブルーチンの名前は 6-7 文字程度で、用いる行列の種類、解きたい問題によってサブルーチンの名前が異なります。初めの一文字は変数型を表しており、• S: 単精度実数• D: 倍精度実数• C: 単精度複素数• Z: 倍精度複素数のいずれかになっています。次の二文字は行列の種類を示しており、• SY: 対称行列
• HE: エルミート行列• GT: 三重対角行列• ST: 対称三重対角行列• GE: 一般行列などがあります。
残りの文字は計算の種類を示しており、例えば線形方程式問題は SV で表され、倍精度実一般行列の線形方程式問題のルーチンは DGESV となります。倍精度複素エルミート行列の線形方程式問題を解く場合はZHESV などとなります。
その他、固有値問題は EV となり、倍精度実対称行列の固有値問題は DSYEV、倍精度一般行列の固有値問
題は DGEEV です。

1

LAPACK で提供されているルーチンのリストは例えば http://www.netlib.org/lapack/explore-html/にあります。これらの中からどのルーチンを使うべきかをまず決定します。

2 LAPACK のインストール
サテライト室のシステムには LAPACK ライブラリはインストールされており、gfortran で使うことができます。自分のコンピュータにインストールすることももちろんできますがここでは省略します。

3 LAPACK を使ったコーディング
LAPACK のルーチンはサブルーチン副プログラムとして呼び出すことができますので CALL 文を一文書くだけでライブラリを使うことができます。引数で計算に用いる行列の次元、行列、計算に一時的に必要となる配列などを形式通り渡す必要がありますので引数の形式を合わせる必要があります。使いたいルーチン名で検索をしても見つけることができますし、LAPACK のルーチンの引数 (およびソースコード) のリストはhttp://www.netlib.org/lapack/explore-html にもあります。
ここでは比較的よく使うと思われる倍精度実行列の線形方程式ルーチン DGESV、倍精度実対称行列の固有値問題ルーチン DSYEV、倍精度実一般行列の固有値問題ルーチン DGEEV について説明します。LAPACKのルーチンには正常終了したかどうかを返す変数がありますが、いずれの場合も正常終了しなかった場合は予定されている値が計算されていないため、サブルーチンを呼んだ後はエラー処理をするようにしてください。

3.1 DGESV
倍精度実一般行列 A による連立一次方程式

AX = Bを解くルーチンです。ここで X や B は一般に行列です。使うときは CALL 文で呼び出します。

CALL DGESV ( N,NRHS,A,LDA,IPIV,B,LDB,INFO )ここで引数はintegerN

integer

NRHS

double precision, dimension( LDA, N )integerA

LDA

integer, dimension( N ) IPIV
double precision, dimension( LDB, NRHS ) B
integer

LDB

integer

INFOの形式です。
このサブルーチンの入力変数は2(1)• N は配列 A の２つ目の次元• NRHS は配列 B の２つ目の次元 (解きたい N 元連立一次方程式の数)。
• A は LDA×N 次元の２次元配列 A(連立一次方程式の係数行列)。
• LDA は行列 A の１つ目の次元 (通常は N とします)。
• B は LDB×NRHS 次元の行列 B 、LDB は行列 B の１つ目の次元 (通常は N とします)。
このサブルーチンの出力変数は• ２次元配列 A には行列 A を代入した状態でサブルーチンに値を渡しますが、実行後には値が書き換えられ、行列 A を LU 分解した結果が代入されます。L 行列と U 行列はそれぞれ下三角行列、上三角行列で、対角要素を除いて一つの N × N 行列に値を入れることができますが、L の対角要素は 1 であるため、U の対角要素が A として返されます。

• ２次元配列 B には行列 B の値を入力しましたが、出力では方程式の解の X の値が返されます。
• IPIV は次元 N の整数型の一次元配列で、ピボットと呼ばれます。連立一次方程式を解く際に数値誤差を小さくするためにサブルーチンの中で行の順番を入れ替えることがありますが、その時の入れ替えた順番を保存するための配列です。

• INFO にはプログラムの実行結果が保存されます。正常終了した場合は 0 が入っています。負の数-i が返された場合は、i 行目で不正な値があったことを示します。また正の数 i が返された場合、LU 分解のU 行列の Uii がゼロとなり (つまり det A = 0 となり) 連立一次方程式の解が求まらなかったことを示します。
演習問題 30、31 で DGESV を使ったサンプルプログラムがありますので参照してください。

3.2 DSYEV
DSYEV は実対称行列 A の固有値問題
AU = U Dを解くサブルーチンです。呼び出すときはCALL DSYEV ( JOBZ,UPLO,N,A,LDA,W,WORK,LWORK,INFO )の引数を指定します。引数の型は

character

JOBZ

character

UPLO

integerN

double precision, dimension( LDA, N )integerA

LDA

double precision, dimension( N )Wdouble precision, dimension( LWORK )integer

LWORK

integer

INFO

WORK3(2)入力時に値を設定する引数は• JOBZ が ’N’ の場合は固有値のみを計算、’V’ の場合は固有値と固有ベクトルを計算します。
• UPLO: A は対称行列ですので’U’ の場合は A の上三角の部分に行列要素が入っているとして計算します。’L’ の場合は下三角の部分の行列要素を使って計算します。

• N: 配列 A の次元• A: 対角化したい行列 A で次元は LDA×N です。UPLO でどちらを指定するかによって計算時に参照される要素が変わります。

• LDA: 配列 A の１つ目の次元ですが通常 N とします。
• WORK: 作業用の倍精度実数型の LWORK 成分ある配列です。
• LWORK: WORK の要素数で、max(1, 3N − 1) よりも大きく取るようにします。より大きい方が計算の効率があがります。
出力は• A: 固有ベクトルを計算するオプションを指定していた場合はルーチン実行後の A には固有ベクトルが格納されます (U 行列が入ります)。固有ベクトルを計算しない場合ももともとの入力時の A の値は破壊されます。

• W: 倍精度実数の N 要素の一次元配列で、小さいものから順に固有値が代入されます。実対称行列の固有値は実数です。

• INFO: 正常終了したら 0 が入ります。ゼロでない値が返された場合はマニュアルを参照してください。
演習問題 32 に DSYEV を使ったサンプルプログラムがあります。

3.3 DGEEVDGEEV では実一般行列の固有値問題を解きます。DSYEV と類似していますが、対称行列に限らないため、左、右固有ベクトルを計算するオプションや、固有値は一般的に複素数となり、引数の構造が少し異なります。

CALL DGEEV ( JOBVL,JOBVR,N,A,LDA,WR,WI,VL,LDVL,VR,LDVR,WORK,LWORK,INFO )引数の形式はcharacter

JOBVL

character

JOBVR

integerN

double precision, dimension( LDA, N )integer

double precision, dimension( * )WRdouble precision, dimension( * )WIdouble precision, dimension( LDVL, * )integer

ALDAVLLDVL

4

double precision, dimension( LDVR, * )integer

VRLDVR

double precision, dimension( * )integer

LWORK

integer

INFO

WORK

入力時に値を設定する変数は• JOBVL 左固有ベクトルを計算する (’V’) かしない (’N’) かを与えます。
• JOBVR 右固有ベクトルを計算する (’V’) かしない (’N’) かを与えます。
• N A 行列の２つ目の次元です。
• A 対角化する LDA×N 行列を入力します。ルーチン実行後は値は破壊されます。
• LDA A 行列の１つ目の次元で通常 N に設定します・• WORK 作業用の要素数 LWORK の一次元配列です。
• LWORK WORK の次元。固有ベクトルを計算する場合は 4N よりも大きな値とします。
出力時に値が返ってくる変数は• WR 一次元倍精度実数配列 (要素数 N)。固有値の実部が入ります。
• WI 一次元倍精度実数配列 (要素数 N)。固有値の虚部が入ります。
• VL LDVL×N の二次元配列。左固有ベクトルを計算するオプションを選んでいた場合、左固有ベクトルが入ります。

• LDVL VL の次元。通常 N に設定します。
• VR LDVR×N の２次元配列。右固有ベクトルを計算するオプションを選んでいた場合、右固有ベクトルが入ります。

• LDVR VR の次元。通常 N に設定します。
• INFO 正常終了した場合はゼロが代入されます。そうでない場合はマニュアルを参照してください。
複素数固有値が出る場合は固有ベクトルへの値の格納の仕方が変わります。詳しくはマニュアルを参照してください。演習問題 33 に DGEEV を使ったサンプルプログラムがあります。

3.4 C 言語での LAPACK を用いたコーディング
LAPACK は Fortran で書かれていますので C 言語からは Fortran のサブルーチンを呼び出すという方法で用います。例えば DGESV を使う場合は C 言語では最後にアンダースコアをつけた関数 dgesv () を呼び出します。引数の形式は基本的に Fortran と同じですが、違いは多次元配列は一次元配列として定義し、すべてポインタ渡しをします。これは Fortran のサブルーチンは参照渡しをしているためです。dgesv の場合はdgesv_(&n, &nrhs, a, &lda, ipiv, b, &ldb, &info);
ここで a は要素数 lda*n の一次元配列、b は要素数 ldb*nhrs の一次元配列です。(配列の先頭要素のアドレ
スはそれぞれ a, b で与えられますが変数のアドレスは変数名に&をつけて与えられます)5

4 コンパイル
4.1 FortranFortran では LAPACK はサブルーチンとして呼び出されます。このままコンパイルすると、このサブルーチンがプログラムの中で定義されていないため、エラーがでます。例えば演習問題 30 のコードをコンパイル
すると$ gfortran exercise30.f90 -o exercise30.exe/tmp/ccsUdakv.o: 関数 ‘MAIN__’ 内:
exercise30.f90:(.text+0x410): ‘dgesv_’ に対する定義されていない参照です
collect2: error: ld returned 1 exit statusといった、中で呼び出されている”DGESV” が見つからないというエラーがでます。LAPACK を使うためにはコンパイル時にライブラリをリンクする必要があります。これはコンパイルオプション -llapack -lblas を使います。-l ライブラリ名 で指定されたライブラリを見に行きますが、LAPACK の中では BLAS が使われていますので BLAS ライブラリも同時にリンクする必要があります。-l とライブラリ名の間にスペースは要りません。

$ gfortran exercise30.f90 -o exercise30.exe -llapack -lblas

4.2 C 言語
C 言語でも通常のコンパイルを行うと$ gcc exercise30.c -o exercise30.exe/tmp/ccOqRgMv.o: 関数 ‘main’ 内:
exercise30.c:(.text+0x30e): ‘dgesv_’ に対する定義されていない参照です
collect2: error: ld returned 1 exit statusとエラーがでます。-llapack -lblas オプションを用いて$ gcc exercise30.c -o exercise30.exe -llapack -lblasのようにコンパイルします。

65 演習問題
10×10 行列
10 98
 9 10 9
 8 9 10
7 89
6 78A=
5 67
4 56
3 45
2 341 23789109876546789109876556789109876456789109873456789109823456789109


12

3

4

5

6

7

8

9
10(3)

を考える。

(30) 行列 A と列ベクトル

1
 2 


 3 


 4 


 5 


b=

 −6 
 −7 
 −8 
 −9 
−10(4)

Ax = b(5)からなる連立一次方程式

を LAPACK の実線形方程式のルーチン DGESV を用いて計算せよ。

(31) 行列 A の行列式 det A および逆行列 A−1 を DGESV を用いて求めよ。
(32) 行列 A を実対称行列の固有値問題のルーチン DSYEV を用いて対角化し、固有値と固有ベクトルを求めよ。

(33) 行列 A を実一般行列の固有値問題のルーチン DGEEV を用いて対角化し、固有値と固有ベクトルを求めよ。

7