blob: 9ac120266a802556a47d2b59a0adcd4f803730ae [file] [log] [blame]
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_SPARSE_TRIANGULARVIEW_H
#define EIGEN_SPARSE_TRIANGULARVIEW_H
namespace Eigen {
/** \ingroup SparseCore_Module
*
* \brief Base class for a triangular part in a \b sparse matrix
*
* This class is an abstract base class of class TriangularView, and objects of type TriangularViewImpl cannot be instantiated.
* It extends class TriangularView with additional methods which are available for sparse expressions only.
*
* \sa class TriangularView, SparseMatrixBase::triangularView()
*/
template<typename MatrixType, unsigned int Mode> class TriangularViewImpl<MatrixType,Mode,Sparse>
: public SparseMatrixBase<TriangularView<MatrixType,Mode> >
{
enum { SkipFirst = ((Mode&Lower) && !(MatrixType::Flags&RowMajorBit))
|| ((Mode&Upper) && (MatrixType::Flags&RowMajorBit)),
SkipLast = !SkipFirst,
SkipDiag = (Mode&ZeroDiag) ? 1 : 0,
HasUnitDiag = (Mode&UnitDiag) ? 1 : 0
};
typedef TriangularView<MatrixType,Mode> TriangularViewType;
protected:
// dummy solve function to make TriangularView happy.
void solve() const;
typedef SparseMatrixBase<TriangularViewType> Base;
public:
EIGEN_SPARSE_PUBLIC_INTERFACE(TriangularViewType)
typedef typename MatrixType::Nested MatrixTypeNested;
typedef typename internal::remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
typedef typename internal::remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
template<typename RhsType, typename DstType>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
if(!(internal::is_same<RhsType,DstType>::value && internal::extract_data(dst) == internal::extract_data(rhs)))
dst = rhs;
this->solveInPlace(dst);
}
/** Applies the inverse of \c *this to the dense vector or matrix \a other, "in-place" */
template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const;
/** Applies the inverse of \c *this to the sparse vector or matrix \a other, "in-place" */
template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const;
};
namespace internal {
template<typename ArgType, unsigned int Mode>
struct unary_evaluator<TriangularView<ArgType,Mode>, IteratorBased>
: evaluator_base<TriangularView<ArgType,Mode> >
{
typedef TriangularView<ArgType,Mode> XprType;
protected:
typedef typename XprType::Scalar Scalar;
typedef typename XprType::StorageIndex StorageIndex;
typedef typename evaluator<ArgType>::InnerIterator EvalIterator;
enum { SkipFirst = ((Mode&Lower) && !(ArgType::Flags&RowMajorBit))
|| ((Mode&Upper) && (ArgType::Flags&RowMajorBit)),
SkipLast = !SkipFirst,
SkipDiag = (Mode&ZeroDiag) ? 1 : 0,
HasUnitDiag = (Mode&UnitDiag) ? 1 : 0
};
public:
enum {
CoeffReadCost = evaluator<ArgType>::CoeffReadCost,
Flags = XprType::Flags
};
explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()), m_arg(xpr.nestedExpression()) {}
inline Index nonZerosEstimate() const {
return m_argImpl.nonZerosEstimate();
}
class InnerIterator : public EvalIterator
{
typedef EvalIterator Base;
public:
EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& xprEval, Index outer)
: Base(xprEval.m_argImpl,outer), m_returnOne(false), m_containsDiag(Base::outer()<xprEval.m_arg.innerSize())
{
if(SkipFirst)
{
while((*this) && ((HasUnitDiag||SkipDiag) ? this->index()<=outer : this->index()<outer))
Base::operator++();
if(HasUnitDiag)
m_returnOne = m_containsDiag;
}
else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer()))
{
if((!SkipFirst) && Base::operator bool())
Base::operator++();
m_returnOne = m_containsDiag;
}
}
EIGEN_STRONG_INLINE InnerIterator& operator++()
{
if(HasUnitDiag && m_returnOne)
m_returnOne = false;
else
{
Base::operator++();
if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer()))
{
if((!SkipFirst) && Base::operator bool())
Base::operator++();
m_returnOne = m_containsDiag;
}
}
return *this;
}
EIGEN_STRONG_INLINE operator bool() const
{
if(HasUnitDiag && m_returnOne)
return true;
if(SkipFirst) return Base::operator bool();
else
{
if (SkipDiag) return (Base::operator bool() && this->index() < this->outer());
else return (Base::operator bool() && this->index() <= this->outer());
}
}
// inline Index row() const { return (ArgType::Flags&RowMajorBit ? Base::outer() : this->index()); }
// inline Index col() const { return (ArgType::Flags&RowMajorBit ? this->index() : Base::outer()); }
inline StorageIndex index() const
{
if(HasUnitDiag && m_returnOne) return internal::convert_index<StorageIndex>(Base::outer());
else return Base::index();
}
inline Scalar value() const
{
if(HasUnitDiag && m_returnOne) return Scalar(1);
else return Base::value();
}
protected:
bool m_returnOne;
bool m_containsDiag;
private:
Scalar& valueRef();
};
protected:
evaluator<ArgType> m_argImpl;
const ArgType& m_arg;
};
} // end namespace internal
template<typename Derived>
template<int Mode>
inline const TriangularView<const Derived, Mode>
SparseMatrixBase<Derived>::triangularView() const
{
return TriangularView<const Derived, Mode>(derived());
}
} // end namespace Eigen
#endif // EIGEN_SPARSE_TRIANGULARVIEW_H