[Eigen] Day 11. Gauss-jordan elimination with pivoting

//gauss-jordan elimination with pivoting.
inline void gauJoElim(mxd::r _m, mxd::r _b){
 
    i::t  idx = 0;
    vxd::t pivot_( _m.rows() );
 
    pivot_.setZero();
    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 b0      = _b.block(  0        ,  0      ,  idx                    ,  _b.cols()   );
        mxd::b b1      = _b.block(  idx      ,  0      ,  1                      ,  _b.cols()   );
        mxd::b b2      = _b.block(  idx + 1  ,  0      ,  _b.rows() - (idx + 1)  ,  _b.cols()   );
 
        mxd::b m_piv   = _m.block(  idx      ,  0      , _m .rows() - (idx)      ,  _m.cols()    );    //matrix to be pivoted
        mxd::b b_piv   = _b.block(  idx      ,  0      , _b .rows() - (idx)      ,  _b.cols()    );
 
        pivot_[idx]    = pivot_m(alpha11.value() , a21);
        m_piv          = pivotM_ip(pivot_.row(idx).value(), pivot_.rows() - idx) * m_piv;
        b_piv          = pivotM_ip(pivot_.row(idx).value(), pivot_.rows() - idx) * b_piv;
 
        a01 = a01/alpha11.value();
        a21 = a21/alpha11.value();
 
        A02 = A02 - a01 * a12;
        A22 = A22 - a21 * a12;
 
        b0  = b0 - a01 * b1;
        b2  = b2 - a21 * b1;
 
        a01 = a01.Zero(  a01.rows(), a01.cols()  );
        a21 = a21.Zero(  a21.rows(), a21.cols()  );
 
 
        idx++;
    }
 
    divGauJoElim(_m,_b);
 
}
 
 

main :

    mxd::t m(3,3);
 
    m <<    1, 2, 3,
            1, 2, 4,
            2, 1, 2;
 
 
 
    mxd::t m2(3,1);
    m2 << 6, 7, 5;
 
    gauJoElim(m, m2);
    co(m);
    co(m2);
 

result :

output = 
1 0 0
0 1 0
0 0 1
output = 
1
2
2