macOSにlapackを導入し線型方程式を解く #C++

macOSにlapackを導入し使えるようにするまでの手順をまとめておく。具体例として連立線形方程式を解くところまで試してみた。
インストール方法を示した後、ビルドできるサンプルコードを示し、その後サンプルコードの内容を解説する。

lapackのインストール

homebrewで簡単にインストールできる。

$ brew install lapack

ただしlapackはKeg-only（/usr/local以下にシンボリックリンクを作らない）なので、lapackへのインクルードパスやリンクパスは自前で指定する必要がある。
brew infoするとどのようにパスを指定するべきか情報がちゃんと表示される。

$ brew info lapack

...
For compilers to find lapack you may need to set:
  export LDFLAGS="-L/usr/local/opt/lapack/lib"
  export CPPFLAGS="-I/usr/local/opt/lapack/include"

サンプルコードのビルド

線型方程式を解くサンプルコードは以下の通り。

dgesv.cpp

#include <iostream>
#include <lapacke.h>

int main (int argc, const char * argv[])
{
   // solve Ax = B
   const lapack_int N = 2, NRHS = 1;
   double a[N*N] = {1,2,3,4};  // NxN matrix
   double b[N*NRHS] = {5,6};   // NxNRHS matrix
   // column major
   // A = (1 3), B = (5), ANS = (-1)
   //     (2 4)      (6)        ( 2)
   const lapack_int LDA = N, LDB = N;
   lapack_int IPIV[N];
   lapack_int info;

   info = LAPACKE_dgesv(LAPACK_COL_MAJOR, N, NRHS, a, LDA, IPIV, b, LDB);

   for (int i = 0; i < N; i++) {
      for (int j = 0; j < NRHS; j++) {
         std::cout << b[i + j*N] << ' ';
      }
      std::cout << "\n";
   }

   return(info);
}

以下の手順でビルドできる。
インクルードパス、リンクパスを指定し、liblapackeをリンクしているだけ。

#!/bin/bash -eux

export LDFLAGS="-L/usr/local/opt/lapack/lib"
export CPPFLAGS="-I/usr/local/opt/lapack/include"
g++ $LDFLAGS $CPPFLAGS dgesv.cpp -llapacke

実行して以下のように結果が表示されれば成功である。

$ ./a.out
-1
2

解説

上記のサンプルコードでは Ax = b という連立線型方程式を解いている。
ここで

A = \begin{pmatrix}
  1 & 3 \\
  2 & 4
\end{pmatrix},
b = \begin{pmatrix}
  5 \\
  6
\end{pmatrix}

であり、解$x$は

x = \begin{pmatrix}
  -1 \\
  2
\end{pmatrix}

である。

まずヘッダファイルをインクルードする。

#include <lapacke.h>

ここではlapackではなく、lapackeを利用する。lapackeとはlapackをC言語から利用しやすいようにラップしたもの。
C,C++から使う場合、特に理由がなければlapackeを使えばよい。brew install lapackでlapackeもインストールされる。

lapackはFortranで実装されているため、CやC++から使う場合に注意を要する点がいくつかある。
まず、lapackを使う時の注意点として行列を"column major"で用意する必要がある。
lapackでは行列要素を１次元配列に保存しているが、i行j列の要素は A[i + j*N]に保存される。
C言語のようなrow-majorな言語では通常A[i][j] == A[i*N + j]であるので注意が必要。

https://www.cc.u-tokyo.ac.jp/events/lectures/104/20180919-4.pdf より引用。"Colomn"は"Column"のタイポだと思われる。

また、lapackのインデックスはFortranと同様に１から始まる。
関数にデータを渡すときには配列の先頭のアドレスを渡せば良いので特に気にする必要はないが、関数からインデックスの値が返ってくる場合、そのインデックスの要素にアクセスするためにはa[idx-1]とする必要がある。
今回のサンプルコードでは特にインデックスを利用している箇所はないが、より複雑な処理をする場合には注意が必要である。

lapackには線形代数演算を処理するルーチン群が用意されているが、それらの名前はある規約に基づいて決められている。
今回使った関数はDGESVというlapackのルーチンを使っているが、最初の"D"がdoubleの型であることを表しており、"GE"が一般(general)の行列、"SV"が線型方程式を解く(SolvE)という処理を表している。この規約を知っていれば行いたい処理がどんな名前なのか検討がつくので、関数を探すときに役に立つだろう。
lapackのルーチンを探す際にはこちらのページが便利である。

さて、lapackの関数名がわかったら、lapackeから使うためにはLAPACKE_{lapackの関数名を小文字にしたもの}という名前の関数を呼べば良い。
今回の場合はLAPACKE_dgesvを呼ぶ事になる。DGESVの引数はこちらのページから参照できる。
https://www.hpc.nec/documents/sdk/SDK_NLC/UsersGuide/man/dgesv.html

lapackeから呼ぶときには以下のように引数を渡す。

lapack_int info = LAPACKE_dgesv(LAPACK_COL_MAJOR, N, NRHS, a, LDA, IPIV, b, LDB)

第一引数に配列が"column major"か"row major"かを指定するフラグがある。
- さきほどcolumn majorで行列要素を準備する必要があると書いたが、実はrow majorで構築してもこのフラグをLAPACK_ROW_MAJORにすることで正しく解くことができる。しかし、その場合には転置の処理が行われるので計算コストを要する。自分のコードでできるだけcolumn majorで書いた方がわかりやすいだろう。
第二引数は行列サイズ。AはNxNの行列になる。
第３引数はBの列サイズ。BはNxNRHSの行列になる。
第４引数は行列のポインタA。入力でもあり、LU分解した後の結果も格納される。
- A = P*L*UとなるようにAが分解され、その結果が格納される。（Lの対角要素は1なので、Aには格納されない。）
第５引数LDAはAのleading dimension。Aとしてある巨大な行列の部分行列を指定する場合に使う。
- 部分行列を使う必要がない場合は、Nと同じ値を指定しておけばよい。
第６引数IPIVはサイズNのintegerの配列の出力。何行目と何行目が入れ替えられたかの情報を持っている。
第７引数は行列Bのポインタ。inputかつoutputであり、実行後の解はこの行列に保存される。
第８引数のLDBは行列Bのleding dimension。通常はNを指定すればよい。
返り値のinfoは正常終了したかどうかのフラグ。