2019/06/29, 2022/02/20
疎行列用のソルバー SparseLU を使っているのは特に意味はありません。なお、#include "Eigen/SparseLU" に関して、このページの下のほうに追加情報があります。
// Ax=b を解く
#include <iostream>
#include "Eigen/SparseLU"
using namespace std;
int main(void){
typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;
// SuperLU (http://crd-legacy.lbl.gov/~xiaoye/SuperLU/) のドキュメントにある例題の行列
double arrayB[] = {1,1,1,1,1}; // 右辺値（普通のC++配列で与える例）
const double s=19, u=21, p=16, e=5, r=18, l=12;
tripletList.push_back(T(0,0, s));
tripletList.push_back(T(1,0, l));
tripletList.push_back(T(4,0, l));
tripletList.push_back(T(1,1, u));
tripletList.push_back(T(2,1, l));
tripletList.push_back(T(4,1, l));
tripletList.push_back(T(0,2, u));
tripletList.push_back(T(2,2, p));
tripletList.push_back(T(0,3, u));
tripletList.push_back(T(3,3, e));
tripletList.push_back(T(3,4, u));
tripletList.push_back(T(4,4, r));
// A行列の組成
Eigen::SparseMatrix<double> A(5,5); // サイズは5行5列
A.setFromTriplets(tripletList.begin(), tripletList.end());
std::cout << "matrix (given)" << std::endl << A << std::endl;
// 右辺値ベクトルの組成
Eigen::VectorXd b = Eigen::Map<Eigen::VectorXd>(&arrayB[0], 5);
std::cout << "right hand side vector (given)" << std::endl << b << std::endl;
// ソルバーオブジェクトの組成とA行列の分解
Eigen::SparseLU<Eigen::SparseMatrix<double> > solver;
solver.compute(A);
if(solver.info()!=Eigen::Success) {
// decomposition failed
std::cout << "decomposition failed" << std::endl;
return -1;
}
// 求解
Eigen::VectorXd x = solver.solve(b);
if(solver.info()!=Eigen::Success) {
// solving failed
std::cout << "solving failed" << std::endl;
return -1;
}
std::cout << "solution vector" << std::endl << x << std::endl;
return 0;
}
[2022/02/20] Gitlab の master branch 550af393 のもの（たぶん 3.4.1 になろうとしている最中のもの）を MS-VC++2017 (Express) でコンパイルする分には、後述の警告レベル 4（/W4）でのコンパイルにおける警告メッセージの問題がすべて解決されたようです。
一方で、 Shift_JIS メインの環境では warning C4819 が出るようになりました。これは当該ファイルが UTF-8 の BOM 無しで保存されており ASCII 外の文字を含むためです。SparseLU に関してだけ言えば該当箇所はコメント以外にないことが判りましたので、当座は無視で構いません。
プロジェクトが Shift_JIS で書かれている場合の当座しのぎとしては、警告が出る Eigen ソースファイルを BOM 付き UTF-8 で保存しなおすと、コンパイラ cl.exe の自動認識が機能してくれるようになって問題なくなります。（Shift_JIS で保存する手もあるが、将来を考えると BOM 付き UTF-8 が良いと思います。）
将来的には自身が書いているソースファイルを全面的に UTF-8 で保存しなおして cl /MD /EHsc /W4 /source-charset:utf-8 EigenTest1.cpp とすれば、Eigen が BOM 無しのままでも問題なくなります。ただし UTF-8 における BOM は Unicode 規格上では「容認されるが、必須でも勧められるものでもない」（from Wikipedia.jp）ということのようであり、何が望ましいのかは解釈が混在しているようです。
以下は Eigen-3.3.x および Eigen 3.4.0（リリース版）に関する古くなりつつある情報です。なお Eigen 3.4.0 から MS-VC++2013 ではコンパイルできなくなったようです。
上記のソースコード（EigenTest1.cpp というファイル名で保存されていると仮定、ただしこのリンク先のファイルは Shift_JIS ですので注意）を MS-VC++2013 (Express) で警告レベル 4（/W4）でコンパイルすると、多数の警告メッセージが出る。（/W3 またはそれ以下だと出ない）
cl /MD /EHsc /W4 EigenTest1.cpp
このうち「warning C4127: 条件式が定数です。」は、Eigen のソースコードを見れば、何ら問題ないことが自明なものばかりである。最初の１個だけを挙げるが、RK は enum の中での定義であるため、コンパイル時には定数になる（ただし、RK の値自体はコンパイラやターゲットマシンに依存して決まることに注意）。
| SparseLU_gemm_kernel.h | |
|---|---|
36 |
RK = NumberOfRegisters>=16 ? 4 : 2, // register blocking
|
77 |
if(RK==4) { b20 = pset1<Packet>(Bc0[2]); }
|
一方、「warning C4512: 代入演算子を生成できません。」は、「コンパイルされたクラスの テンプレート のインスタンス化：～～～の参照を確認してください」という膨大な量の警告メッセージに埋もれてしまっている。
「c:\_\eigentest\eigen\src/SparseLU/SparseLU_SupernodalMatrix.h(156) : 'Eigen::internal::MappedSuperNodalMatrix<double,int>::InnerIterator' の宣言を確認してください。」という指摘はされているものの、正直なところ、Eigen のエンドユーザーとしては、Eigen 自体の修正を待つ以外に為す術がない。警告レベル 3（/W3）にデグレードすればこのメッセージを見なくても済むとはいえ、C/C++ の初心者プログラマが陥りやすい様々な潜在バグの発見には、警告レベル 4（/W4）の設定が欠かせない。
警告はともかくとして Eigen の SparseLU の現状の実装が動作自体は正しいのであろうと信用するしかないが、以下のようにすると、Eigen のインクルードに関してだけ局部的に警告レベルのデグレードができる。これで、無為な警告メッセージを抑制し、他の（自分が記述した）部分での潜在バグ発見を重視できる。
| 仮対処 |
#pragma warning(push) #pragma warning(disable:4127) // 「warning C4127: 条件式が定数です。」 #pragma warning(disable:4512) // 「warning C4512: 代入演算子を生成できません。」 #include "Eigen/SparseLU" #pragma warning(pop) |
pragma による警告メッセージの制御については以下のリンクを参照。
supress が「非表示します。」に化けてしまう。
この文書が誰か・何かの役に立てば幸いです。さらに願わくば、この文書の公開（並びに Lis / SuperLU / MUMPS の各事例）が、日本における数値計算ライブラリ利用裾野を広げ、日本全体の底上げにつながることを、静かに期待します。
トップページ /
計算機関係の雑多なこと /
Eigenで線形方程式を解く使い方のごく簡単な例
このサイトについて・著者連絡先