C++の行列演算ライブラリの A.lu().solve(b);と同じ処理をC言語で書きたいです。
自分の理解のためにも、ライブラリを利用しないで実装をしたいのですが、この処理はLU分解を用いた連立方程式を解く処理と思っていたのですが、下記のサイトで紹介している実装で同じ処理をしても結果が一致しませんでした。
https://workspacememory.hatenablog.com/entry/2015/09/27/172005
MatrixXd A(DIM, DIM);
MatrixXd B(DIM, 1);
MatrixXd x(DIM, DIM);
A << 1, 2, 3, 4, 5, 6, 7, 8, 9;
B << 1, 2, 3;
x = A.lu().solve(B);
の結果
-0.3333333333
0.6666666667
0.0000000000
になるのですが、
double A[dim * dim] = {
1, 2, 3,
4, 5, 6,
7, 7, 9
};
double b[dim] = {
1,
2,
3,
};
int pivot[dim];
luDecomposition(dim, pivot, A);
luSolver(dim, pivot, A, b);
の結果が
6.34413e-17,
-5.55112e-17,
0.333333,
と一致しません。
他のサイトの実装を確認してみましたが、Eigenじゃない方の結果と一致します。
Eigenのsolve()の関数コメントを読んでみると、LU分解して連立方程式の解を得るように読めるのですが、
C言語で同じ結果を得るためにはどのようにしたらいいのでしょうか?
double A
が1, ..., 7, 7, 9
となっていますが、書き間違いでしょうか?