Eigen::SkylineInplaceLU< MatrixType > Class Template Reference

Inplace LU decomposition of a skyline matrix and associated features. More...

Public Member Functions
void	compute ()
void	computeRowMajor ()
int	flags () const
int	orderingMethod () const
RealScalar	precision () const
void	setFlags (int f)
void	setOrderingMethod (int m)
void	setPrecision (RealScalar v)
	SkylineInplaceLU (MatrixType &matrix, int flags=0)
template<typename BDerived , typename XDerived >
bool	solve (const MatrixBase< BDerived > &b, MatrixBase< XDerived > *x, const int transposed=0) const
bool	succeeded (void) const

Protected Types
typedef MatrixType::Index	Index
typedef NumTraits< typename MatrixType::Scalar >::Real	RealScalar
typedef MatrixType::Scalar	Scalar

Protected Attributes
int	m_flags
MatrixType &	m_lu
RealScalar	m_precision
int	m_status
bool	m_succeeded

Detailed Description

template<typename MatrixType>
class Eigen::SkylineInplaceLU< MatrixType >

Inplace LU decomposition of a skyline matrix and associated features.

Parameters

MatrixType

the type of the matrix of which we are computing the LU factorization

Definition at line 27 of file SkylineInplaceLU.h.

Member Typedef Documentation

Index

template<typename MatrixType >

typedef MatrixType::Index Eigen::SkylineInplaceLU< MatrixType >::Index

protected

Definition at line 30 of file SkylineInplaceLU.h.

RealScalar

template<typename MatrixType >

typedef NumTraits<typename MatrixType::Scalar>::Real Eigen::SkylineInplaceLU< MatrixType >::RealScalar

protected

Definition at line 32 of file SkylineInplaceLU.h.

Scalar

template<typename MatrixType >

typedef MatrixType::Scalar Eigen::SkylineInplaceLU< MatrixType >::Scalar

protected

Definition at line 29 of file SkylineInplaceLU.h.

Constructor & Destructor Documentation

SkylineInplaceLU()

template<typename MatrixType >

Eigen::SkylineInplaceLU< MatrixType >::SkylineInplaceLU ( MatrixType &  matrix, int  flags = 0  )

inline

Creates a LU object and compute the respective factorization of matrix using flags flags.

Definition at line 38 of file SkylineInplaceLU.h.

    : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0), m_lu(matrix) {
        m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar > ();
        m_lu.IsRowMajor ? computeRowMajor() : compute();
    }

Member Function Documentation

compute()

template<typename MatrixType >

void Eigen::SkylineInplaceLU< MatrixType >::compute

Computes/re-computes the LU factorization

Computes / recomputes the in place LU decomposition of the SkylineInplaceLU. using the default algorithm.

Definition at line 122 of file SkylineInplaceLU.h.

                                           {
    const size_t rows = m_lu.rows();
    const size_t cols = m_lu.cols();
 
    eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
    eigen_assert(!m_lu.IsRowMajor && "LU decomposition does not work with rowMajor Storage");
 
    for (Index row = 0; row < rows; row++) {
        const double pivot = m_lu.coeffDiag(row);
 
        //Lower matrix Columns update
        const Index& col = row;
        for (typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) {
            lIt.valueRef() /= pivot;
        }
 
        //Upper matrix update -> contiguous memory access
        typename MatrixType::InnerLowerIterator lIt(m_lu, col);
        for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
            typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
            typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
            const double coef = lIt.value();
 
            uItPivot += (rrow - row - 1);
 
            //update upper part  -> contiguous memory access
            for (++uItPivot; uIt && uItPivot;) {
                uIt.valueRef() -= uItPivot.value() * coef;
 
                ++uIt;
                ++uItPivot;
            }
            ++lIt;
        }
 
        //Upper matrix update -> non contiguous memory access
        typename MatrixType::InnerLowerIterator lIt3(m_lu, col);
        for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
            typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
            const double coef = lIt3.value();
 
            //update lower part ->  non contiguous memory access
            for (Index i = 0; i < rrow - row - 1; i++) {
                m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef;
                ++uItPivot;
            }
            ++lIt3;
        }
        //update diag -> contiguous
        typename MatrixType::InnerLowerIterator lIt2(m_lu, col);
        for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
 
            typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
            typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
            const double coef = lIt2.value();
 
            uItPivot += (rrow - row - 1);
            m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef;
            ++lIt2;
        }
    }
}

computeRowMajor()

template<typename MatrixType >

void Eigen::SkylineInplaceLU< MatrixType >::computeRowMajor

Definition at line 186 of file SkylineInplaceLU.h.

                                                   {
    const size_t rows = m_lu.rows();
    const size_t cols = m_lu.cols();
 
    eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
    eigen_assert(m_lu.IsRowMajor && "You're trying to apply rowMajor decomposition on a ColMajor matrix !");
 
    for (Index row = 0; row < rows; row++) {
        typename MatrixType::InnerLowerIterator llIt(m_lu, row);
 
 
        for (Index col = llIt.col(); col < row; col++) {
            if (m_lu.coeffExistLower(row, col)) {
                const double diag = m_lu.coeffDiag(col);
 
                typename MatrixType::InnerLowerIterator lIt(m_lu, row);
                typename MatrixType::InnerUpperIterator uIt(m_lu, col);
 
 
                const Index offset = lIt.col() - uIt.row();
 
 
                Index stop = offset > 0 ? col - lIt.col() : col - uIt.row();
 
                //#define VECTORIZE
#ifdef VECTORIZE
                Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
                Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
 
 
                Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal);
#else
                if (offset > 0) //Skip zero value of lIt
                    uIt += offset;
                else //Skip zero values of uIt
                    lIt += -offset;
                Scalar newCoeff = m_lu.coeffLower(row, col);
 
                for (Index k = 0; k < stop; ++k) {
                    const Scalar tmp = newCoeff;
                    newCoeff = tmp - lIt.value() * uIt.value();
                    ++lIt;
                    ++uIt;
                }
#endif
 
                m_lu.coeffRefLower(row, col) = newCoeff / diag;
            }
        }
 
        //Upper matrix update
        const Index col = row;
        typename MatrixType::InnerUpperIterator uuIt(m_lu, col);
        for (Index rrow = uuIt.row(); rrow < col; rrow++) {
 
            typename MatrixType::InnerLowerIterator lIt(m_lu, rrow);
            typename MatrixType::InnerUpperIterator uIt(m_lu, col);
            const Index offset = lIt.col() - uIt.row();
 
            Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
 
#ifdef VECTORIZE
            Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
            Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
 
            Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal);
#else
            if (offset > 0) //Skip zero value of lIt
                uIt += offset;
            else //Skip zero values of uIt
                lIt += -offset;
            Scalar newCoeff = m_lu.coeffUpper(rrow, col);
            for (Index k = 0; k < stop; ++k) {
                const Scalar tmp = newCoeff;
                newCoeff = tmp - lIt.value() * uIt.value();
 
                ++lIt;
                ++uIt;
            }
#endif
            m_lu.coeffRefUpper(rrow, col) = newCoeff;
        }
 
 
        //Diag matrix update
        typename MatrixType::InnerLowerIterator lIt(m_lu, row);
        typename MatrixType::InnerUpperIterator uIt(m_lu, row);
 
        const Index offset = lIt.col() - uIt.row();
 
 
        Index stop = offset > 0 ? lIt.size() : uIt.size();
#ifdef VECTORIZE
        Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
        Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
        Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal);
#else
        if (offset > 0) //Skip zero value of lIt
            uIt += offset;
        else //Skip zero values of uIt
            lIt += -offset;
        Scalar newCoeff = m_lu.coeffDiag(row);
        for (Index k = 0; k < stop; ++k) {
            const Scalar tmp = newCoeff;
            newCoeff = tmp - lIt.value() * uIt.value();
            ++lIt;
            ++uIt;
        }
#endif
        m_lu.coeffRefDiag(row) = newCoeff;
    }
}

flags()

template<typename MatrixType >

int Eigen::SkylineInplaceLU< MatrixType >::flags ( ) const

inline

Returns

the current flags

Definition at line 78 of file SkylineInplaceLU.h.

                      {
        return m_flags;
    }

orderingMethod()

template<typename MatrixType >

int Eigen::SkylineInplaceLU< MatrixType >::orderingMethod ( ) const

inline

Definition at line 86 of file SkylineInplaceLU.h.

                               {
        return m_flags;
    }

precision()

template<typename MatrixType >

RealScalar Eigen::SkylineInplaceLU< MatrixType >::precision ( ) const

inline

Returns

the current precision.

See also

setPrecision()

Definition at line 61 of file SkylineInplaceLU.h.

                                 {
        return m_precision;
    }

setFlags()

template<typename MatrixType >

void Eigen::SkylineInplaceLU< MatrixType >::setFlags ( int  f )

inline

Sets the flags. Possible values are:

• CompleteFactorization
• IncompleteFactorization
• MemoryEfficient
• one of the ordering methods
• etc...

See also

flags()

Definition at line 73 of file SkylineInplaceLU.h.

                         {
        m_flags = f;
    }

setOrderingMethod()

template<typename MatrixType >

void Eigen::SkylineInplaceLU< MatrixType >::setOrderingMethod ( int  m )

inline

Definition at line 82 of file SkylineInplaceLU.h.

                                  {
        m_flags = m;
    }

setPrecision()

template<typename MatrixType >

void Eigen::SkylineInplaceLU< MatrixType >::setPrecision ( RealScalar  v )

inline

Sets the relative threshold value used to prune zero coefficients during the decomposition.

Setting a value greater than zero speeds up computation, and yields to an incomplete factorization with fewer non zero coefficients. Such approximate factors are especially useful to initialize an iterative solver.

Note that the exact meaning of this parameter might depends on the actual backend. Moreover, not all backends support this feature.

See also

precision()

Definition at line 54 of file SkylineInplaceLU.h.

                                    {
        m_precision = v;
    }

solve()

template<typename MatrixType >

template<typename BDerived , typename XDerived >

bool Eigen::SkylineInplaceLU< MatrixType >::solve	(	const MatrixBase< BDerived > &	b,
		MatrixBase< XDerived > *	x,
		const int	transposed = 0
	)		const

Returns

the lower triangular matrix L

the upper triangular matrix U

Computes *x = U^-1 L^-1 b

If transpose is set to SvTranspose or SvAdjoint, the solution of the transposed/adjoint system is computed instead.

Not all backends implement the solution of the transposed or adjoint system.

Definition at line 309 of file SkylineInplaceLU.h.

                                                                                                                           {
    const size_t rows = m_lu.rows();
    const size_t cols = m_lu.cols();
 
 
    for (Index row = 0; row < rows; row++) {
        x->coeffRef(row) = b.coeff(row);
        Scalar newVal = x->coeff(row);
        typename MatrixType::InnerLowerIterator lIt(m_lu, row);
 
        Index col = lIt.col();
        while (lIt.col() < row) {
 
            newVal -= x->coeff(col++) * lIt.value();
            ++lIt;
        }
 
        x->coeffRef(row) = newVal;
    }
 
 
    for (Index col = rows - 1; col > 0; col--) {
        x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col);
 
        const Scalar x_col = x->coeff(col);
 
        typename MatrixType::InnerUpperIterator uIt(m_lu, col);
        uIt += uIt.size()-1;
 
 
        while (uIt) {
            x->coeffRef(uIt.row()) -= x_col * uIt.value();
            //TODO : introduce --operator
            uIt += -1;
        }
 
 
    }
    x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0);
 
    return true;
}

succeeded()

template<typename MatrixType >

bool Eigen::SkylineInplaceLU< MatrixType >::succeeded ( void  ) const

inline

Returns

true if the factorization succeeded

Definition at line 105 of file SkylineInplaceLU.h.

                                      {
        return m_succeeded;
    }

Member Data Documentation

m_flags

template<typename MatrixType >

int Eigen::SkylineInplaceLU< MatrixType >::m_flags

protected

Definition at line 111 of file SkylineInplaceLU.h.

m_lu

template<typename MatrixType >

MatrixType& Eigen::SkylineInplaceLU< MatrixType >::m_lu

protected

Definition at line 114 of file SkylineInplaceLU.h.

m_precision

template<typename MatrixType >

RealScalar Eigen::SkylineInplaceLU< MatrixType >::m_precision

protected

Definition at line 110 of file SkylineInplaceLU.h.

m_status

template<typename MatrixType >

int Eigen::SkylineInplaceLU< MatrixType >::m_status

mutableprotected

Definition at line 112 of file SkylineInplaceLU.h.

m_succeeded

template<typename MatrixType >

bool Eigen::SkylineInplaceLU< MatrixType >::m_succeeded

protected

Definition at line 113 of file SkylineInplaceLU.h.

---

The documentation for this class was generated from the following file:

• SkylineInplaceLU.h