how to use eigen library for lu decomposition c++

Question

I am using the standard Eigen library to compute the LU decomposition of a given matrix.

Yet I get confused about some of the functions. The problem is based on changing A to LU (A=LU).

I thought I know the logic behind LU decomposition. First step check whether u can do it directly, if not change the original matrix A to PA and do PA=LU; And I thought this step is being done by calling the function matrixLU(), which is exactly what is done in example 2. However this function gives me something that I don't understand in example 1 and 3.

Example 1 works, question is what's the use of the function matrixLU()?

Example 2 works, and I do know the use of matrixLU() here. (change to PA=LU)

Example 3 does not work, the result LU gives a matrix which changes the row of the original matrix. Why this one gives a wrong result?

Appreciate any explanations or math logic behind the function matrixLU().

#include <iostream>
#include <Eigen/Dense>
//#include <Dense>
#include <Eigen/LU>
using Eigen::MatrixXd;
using std::cout;
using std::endl;
using Eigen::VectorXd;
using std::cin;
int main()
{

    //example 1
    typedef Eigen::Matrix<double,4,4> M4x4;
    M4x4 p;
    p<< 7,3,-1,2,3,8,1,-1,-1,1,4,-1,2,-4,-1,6;
    cout<<p<<endl<<endl;
    // Create LU Decomposition template object for p
    Eigen::PartialPivLU<M4x4> LU(p);
    cout<<"LU MATRIX:\n"<<LU.matrixLU()<<endl<<endl;
    // Output L, the lower triangular matrix
    M4x4 l = MatrixXd::Identity(4,4);//默认 单位对角矩阵
    //开始填充
    l.block<4,4>(0,0).triangularView<Eigen::StrictlyLower>()=LU.matrixLU();
    cout<<"L MATRIX:\n"<<l<<endl <<endl;
    M4x4 u = LU.matrixLU().triangularView<Eigen::Upper>();
    cout<<"R MATRIX:\n"<<u<<endl<<endl;
    MatrixXd m0(4,4);
    m0 = l*u;
    cout<<"calculate the original matrix:\n"<<m0<<endl<<endl;
    //证明完毕
    //example 2
    typedef Eigen::Matrix<double,2,2> M2X2;
    M2X2 p0;
    p0<< 0,2,1,3;
    cout<<p0<<endl<<endl;
    Eigen::PartialPivLU<M2X2> LU0(p0);
    cout<<"LU MATRIX:\n"<<LU0.matrixLU()<<endl<<endl;//原来是在做PA的过程
    //一切结果从PA开始
    M2X2 l0 = MatrixXd::Identity(2,2);
    l0.block<2,2>(0,0).triangularView<Eigen::StrictlyLower>()=LU0.matrixLU();
    cout<<"L MATRIX:\n"<<l0<<endl <<endl;
    //以下省略N行
    //example 3
    typedef Eigen::Matrix<double,3,3> M3X3;
    M3X3 p1;
    p1<<3,-1,2,6,-1,5,-9,7,3;
    cout<<p1<<endl<<endl;
    Eigen::PartialPivLU<M3X3> LU1(p1);
    cout<<"LU MATRIX:\n"<<LU1.matrixLU()<<endl<<endl;//暂时没明白这步做的啥
    M3X3 l1 = MatrixXd::Identity(3,3);
    l1.block<3,3>(0,0).triangularView<Eigen::StrictlyLower>()=LU1.matrixLU();
    cout<<"L MATRIX:\n"<<l1<<endl <<endl;
    //直接up
    M3X3 u1 = LU1.matrixLU().triangularView<Eigen::Upper>();
    cout<<"R MATRIX:\n"<<u1<<endl<<endl;
    cout<<l1*u1<<endl;
    cin.get();
}

Answer 1

The factorization is PA=LU meaning that the product LU represents the matrix A after some row interchanges. This is consistent with all your three examples. See the respective documentation for more details.