4−４．線形数値演算パッケージLAPACK

この章ではコンピューターを用いて行列演算を行う方法について学習した。物理で扱う行列は要素の数が莫大になることが多く、その演算は精度と速度の両方が要求される。どちらも単に安定で演算量 の少ないアルゴリズムを用いればよいというものではなく、コンピューターのハードウエアやオペレーティングシステム、コンパイラの仕様に大きく依存した（本来の物理には全く本質的ではない）プログラミング技術が要求される。したがって行列演算に限っては、自分で苦労してプログラムを開発するより、すでに公開されている汎用のパッケージを用いることを強く勧める。そのほうが「高速で信頼性の高い」結果を得ることができるからである。そのようなパッケージとして

などがよく知られており、この講義ではその中でも最も普及していると思われるLAPACKを取り上げる。

　

LAPACK（レイパック）を使った線形数値演算

１．LAPACKとは

LAPACK(Linear Algebra PACKage、レイパックと読む)は、netlib [original, mirror]で公開されている高速で信頼性の高い（しかもフリーの）線形 演算ライブラリであり、線形方程式の直接解法に関しては最も優れた選択肢の一つである。京都大学をはじめ、世界中の主な計算機センターのマシンにははじめからインストールされていることが多く、その場合にはコンパイル時に ライブラリをリンクするだけですべての機能を利用することができる。

LAPACKは

連立1次方程式／線形最小二乗問題／固有値問題／特異値分解

を扱うことができる（残念ながらFFTは含まれない）。ソースがFORTRAN 77で記述されたサブルーチン集であるが、コンパイル済みのライブラリをリンクすることにより、各々のサブルーチンをC言語から呼び出して使うこともできる。またCコンパイラのみの開発環境でもライブラリを作成できるように、f2cを用いてソースをC言語化したCLAPACKというパッケージ、あるいはC++への対応を目指したLAPACK++というパッケージ（機能は限定される）も存在する。

LAPACKは，内部でBLAS (Basic Linear Algebra Subprograms) ライブラリを呼び出す。OSやCPUに依存する最適化をこのBLASによって吸収するしくみである。自分の使用するCPUに対して最適化されたBLASを用いれば、計算のパフォーマンスが飛躍的に向上する。BLAS には様々な種類があり、各CPUベンダーから最適化されたBLASが提供されている。 またそれとは別にインテルPentiumプロセッサー上のLinux用のBLASはここにある。たとえ自分のシステムに適合したBLASが見つからなくても、各システムに最適化したBLASを生成するATLAS というソフトが存在する。

２．LAPACKのインストール

ここではメディアセンターのLinux環境において、自分のディスク領域にBLASとLAPACKのコンパイル済みのライブラリを導入して使用する方法を説明する。メディアセンターでは個人に５０MBのディスクを割り当てているが、このライブラリだけで約１０MBの空きスペースを要するので注意すること。本来ならば/usr/libなどにインストールして共有するのが望ましいが、私が管理者権限を持たないので今回はこのような方法をとらざるを得ない。

  1. Linux上のg77、gccで動作するPentium用のコンパイル済みライブラリを用意した。http://stat.scphys.kyoto-u.ac.jp/~ryoichi/CP/lapack/lapack.tar.gz にアクセスし、lapack.tar.gzを任意のディレクトリ（例えば ~/lapack）に保存する。

    > mkdir lapack
    > cd lapack
    > netscape &　　（netscapeから上のファイルにアクセスし、保存）
    　
  2. 保存したディレクトリで以下を実行してライブラリを解凍し、不要なファイルを削除する。

    > tar xvfz lapack.tar.gz
    > ls
    lapack_linux.a   libblas.a    lsblaspii1.2f_03.00.a
    lapack.tar.gz    liblapack.a
    > rm lapack.tar.gz
    　
  3. 同じディレクトリにサンプルプログラムを保存し、以下を実行してちゃんと動作するか確認。

    [Fortran]
    > f77 sample.f -L. -llapack -lblas
    > ./a.out
    19.00 3.00 1.00 12.00 1.00 16.00 1.00 3.00 11.00 0.00
    -19.00 3.00 1.00 12.00 1.00 16.00 1.00 3.00 11.00 0.00
    -19.00 -3.00 1.00 12.00 1.00 16.00 1.00 3.00 11.00 1.00
    -19.00 -3.00 -1.00 12.00 1.00 16.00 1.00 3.00 11.00 0.00
    -19.00 -3.00 -1.00 -12.00 1.00 16.00 1.00 3.00 11.00 0.00
    -19.00 -3.00 -1.00 -12.00 -1.00 16.00 1.00 3.00 11.00 0.00
    -19.00 -3.00 -1.00 -12.00 -1.00 -16.00 1.00 3.00 11.00 0.00
    -19.00 -3.00 -1.00 -12.00 -1.00 -16.00 -1.00 3.00 11.00 1.00
    -19.00 -3.00 -1.00 -12.00 -1.00 -16.00 -1.00 -3.00 11.00 0.00
    DGESV works fine
    X 1 = 0.
    X 2 = -0.166666667
    X 3 = 0.5
    X 4 = 0.
    X 5 = 0.
    X 6 = 0.
    X 7 = -0.5
    X 8 = 0.166666667
    X 9 = 0.

    [C]
    > gcc sample.c -L. -L/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66 -llapack -lblas -lm -lg2c
    > ./a.out
    N = 3
    1.000000
    1.000000
    1.000000
    ただしコンパイルオプションの -L. は、カレントディレクトリ( . )をリンクされるべきライブラリのディレクトリリストに追加することを表し、-llapack -lblas はそれぞれ liblapack.a libblas.a をリンクすることを表す。

http://www.intel.com/software/products/
（ここのIntel Math Kernel LibraryにはPentium用に最適化されたBLASとLAPACKが含まれる。Linux用、Windows用があり、どちらも研究用途でのダウンロードが認められている。要登録。）

http://phase.hpcc.jp/phase/lapack-j/LAPACK3.0/lapack3.0_howto.html

http://www.view.human.nagoya-u.ac.jp/~hiroyuki/lab/clapack/

http://soil.en.a.u-tokyo.ac.jp/~tokida/Lapack.htm

http://www.a.mei.titech.ac.jp/~yonishi/Program/CLapack/

http://www.miv.t.u-tokyo.ac.jp/~yabuki/tip/lapack/lapack.htm

　

３．LAPACKの使い方（公式マニュアル）

LAPACKはDriver（一般な問題を解くための機能を提供する）、Computational（個々の問題を解くための機能を提供する）、およびAuxiliary（補助的な計算や共通に使用される手続きを提供する）というルーチンで構成されているが、ほとんどの場合Driverルーチンを使えば事が足りる。

各Driverルーチンの機能の説明、呼び出し方の一覧が京大の計算機センターのページにある[html, pdf]。扱う行列の型（一般、対称、３重対角…）や、変数の型（実数、複素数）、精度（単精度、倍精度）によって実に多くのサブルーチンが用意されているが、まずはこの中から自分の目的に合ったものを見つけ出すことが必要である。上のリストを参考にして

  1. 単度実数(S)の一般行列(GE)で与えられた連立方程式(SV)を解くためには SGESV
  2. 倍精度複素数(Z)の一般行列(GE)で与えられた対称固有値問題(EV)を解くには ZGEEV

とった風にすればよい。

Fortran、C言語でそれぞれサブルーチンの呼び方と行列を指定する方法が異なるので注意が必要である。次節で具体例をあげて解説する。