gsMatrixAddons.h

Go to the documentation of this file.

inline const internal::adjugate_impl<Derived> adjugate() const;

inline void adjugateInPlace();
    
typedef BlockDiag<Derived,Dynamic> BlockDiagReturnType;
inline const BlockDiagReturnType blockDiag(Index rowFactor) const;

typedef BlockTranspose<Derived,Dynamic> BlockTransposeReturnType;
inline const BlockTransposeReturnType blockTranspose(Index rowFactor) const;

template<typename IndicesType>
const RowSelection<Derived,IndicesType> selectRows(const IndicesType & ind) const;


void gaussElim()
{
    Derived & M = derived();
    index_t piv = 0;
    const index_t nr = M.rows();
    const index_t nc = M.cols();
    for (index_t r=0; r!=nr; ++r)
    {
        piv = 0;
        while (piv!= nc && 0 == M(r, piv)) ++piv;               
        if (piv == nc ) continue;
        
        const index_t br = nr-r-1;
        const index_t bc = nc-piv-1;
        M.block(r+1, piv+1, br, bc).noalias() -=
            M.col(piv).tail(br) * M.row(r).tail(bc) / M(r, piv);
        M.col(piv).tail(br).setZero();
    }
}


inline Matrix<Scalar, Dynamic, Dynamic> cramerInverse() const
{
    const Derived & M = derived();
    eigen_assert(M.rows() == M.cols() && "Matrix is not square.");

    Matrix<Scalar, Dynamic, Dynamic> rvo1(M.rows(), M.rows());
    switch (M.rows())
    {
        case 1:
            rvo1 = M.template topLeftCorner<1, 1>().inverse();
            break;
        case 2:
            rvo1 = M.template topLeftCorner<2, 2>().inverse();
            break;
        case 3:
            rvo1 = M.template topLeftCorner<3, 3>().inverse();
            break;
        case 4:
            rvo1 = M.template topLeftCorner<4, 4>().inverse();
            break;
        default:
            gsWarn<<"Inversion by LU for matrix of size "<<M.rows()<<"\n";
            M.inverse();
            break;
    };
    return rvo1;
}

inline void cramerInverseInPlace()
{
//    derived() = cramerInverse().eval();
    derived().swap(cramerInverse().eval());
}