0

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言語で同じ結果を得るためにはどのようにしたらいいのでしょうか?

2
  • 2
    double A1, ..., 7, 7, 9 となっていますが、書き間違いでしょうか?
    – user39889
    2021年4月17日 15:50
  • これでしたね、書き間違いでした。恥ずかしい間違いですね。ありがとうございます。 2021年4月19日 11:49

0

この質問に回答するには、ログインする必要があります。

のタグが付いた他の質問を参照する。