[Eigen] Day 10. LU decomposition with pivoting

 
//lu decomposition with pivoting
inline vxd::t lu_m(mxd::r _m){
 
    vxd::t pivot_( _m.rows() );
    pivot_.setZero();
 
    i::t idx = 0;
 
    while(idx < _m.rows()){
 
        mxd::b A00     = _m.block(0      , 0       , idx                     , idx                    );
        mxd::b a10     = _m.block(idx    , 0       , 1                       , idx                    );
        mxd::b A20     = _m.block(idx + 1, 0       , _m .rows() - (idx + 1)  , idx                    );
 
        mxd::b a01     = _m.block(0      , idx     , idx                     , 1                      );
        mxd::b alpha11 = _m.block(idx    , idx     , 1                       , 1                      );
        mxd::b a21     = _m.block(idx + 1, idx     , _m.rows() - (idx + 1)   , 1                      );
 
        mxd::b A02     = _m.block(0      , idx + 1 , idx                     , _m.cols() - (idx + 1)  );
        mxd::b a12     = _m.block(idx    , idx + 1 , 1                       , _m.cols() - (idx + 1)  );
        mxd::b A22     = _m.block(idx + 1, idx + 1 , _m.rows() - (idx + 1)   , _m.cols() - (idx + 1)  );
 
        mxd::b piv     = _m.block(idx, 0 , _m .rows() - (idx)            , _m.cols() );
 
        pivot_[idx]    = pivot_m(alpha11.value() , a21);
        piv = pivotM_ip(pivot_.row(idx).value(), pivot_.rows() - idx) * piv;
 
        a21 = a21 / alpha11.value();
        A22 = A22 - a21 * a12;
 
        idx++;
    }
    return pivot_;
}
 
 

main

 
    mxd::t m2(3,3);
 
    m2 <<   1, 2, 3,
            1, 2, 4,
            2, 1, 2;
 
 
    co(lu_m(m2));
    co(m2);
 
check validity:

    mxd::t l(3,3);
    mxd::t u(3,3);
    mxd::t p(3,3);
 
    l <<    1,0,0,
            2,1,0,
            1,0,1;
 
    u <<    1, 2, 3,
            0,-3,-4,
            0, 0, 1;
 
    p <<    1,0,0,
            0,0,1,
            0,1,0;
 
co( p * l * u,"p * l * u");
 

result : 

output = 
0
1
0
output = 
 1  2  3
 2 -3 -4
 1 -0  1
p * l * u = 
1 2 3
1 2 4
2 1 2