

C++の行列演算ライブラリの A.lu().solve(b);と同じ処理をC言語で書きたいです。  
自分の理解のためにも、ライブラリを利用しないで実装をしたいのですが、この処理はLU分解を用いた連立方程式を解く処理と思っていたのですが、下記のサイトで紹介している実装で同じ処理をしても結果が一致しませんでした。  
<https://workspacememory.hatenablog.com/entry/2015/09/27/172005>

[code]

        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);
    
[/code]

の結果

[code]

    -0.3333333333
     0.6666666667
     0.0000000000
    
[/code]

になるのですが、

[code]

        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);
    
[/code]

の結果が

[code]

    6.34413e-17, 
    -5.55112e-17, 
    0.333333,
    
[/code]

と一致しません。  
他のサイトの実装を確認してみましたが、Eigenじゃない方の結果と一致します。  
Eigenのsolve()の関数コメントを読んでみると、LU分解して連立方程式の解を得るように読めるのですが、  
C言語で同じ結果を得るためにはどのようにしたらいいのでしょうか？


CC BY-SA 4.0
2021-04-17T14:39:40.930
75348
2021-04-17T14:39:40.930
44859
question
c++ c 行列
Eigenで　A.lu().solve(b);　とした処理をC言語で実現したいのですがどのように書いたらいいでしょうか？