Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
index 4c169aa..338e6f1 100644
--- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
+++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
@@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
+// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@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
@@ -24,7 +25,7 @@
* \param ind The array of index for the elements in @p row
* \param ncut The number of largest elements to keep
**/
-template <typename VectorV, typename VectorI, typename Index>
+template <typename VectorV, typename VectorI>
Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
{
typedef typename VectorV::RealScalar RealScalar;
@@ -66,6 +67,8 @@
* \class IncompleteLUT
* \brief Incomplete LU factorization with dual-threshold strategy
*
+ * \implsparsesolverconcept
+ *
* During the numerical factorization, two dropping rules are used :
* 1) any element whose magnitude is less than some tolerance is dropped.
* This tolerance is obtained by multiplying the input tolerance @p droptol
@@ -92,28 +95,36 @@
* alternatively, on GMANE:
* http://comments.gmane.org/gmane.comp.lib.eigen/3302
*/
-template <typename _Scalar>
-class IncompleteLUT : internal::noncopyable
+template <typename _Scalar, typename _StorageIndex = int>
+class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar, _StorageIndex> >
{
+ protected:
+ typedef SparseSolverBase<IncompleteLUT> Base;
+ using Base::m_isInitialized;
+ public:
typedef _Scalar Scalar;
+ typedef _StorageIndex StorageIndex;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar,Dynamic,1> Vector;
- typedef SparseMatrix<Scalar,RowMajor> FactorType;
- typedef SparseMatrix<Scalar,ColMajor> PermutType;
- typedef typename FactorType::Index Index;
+ typedef Matrix<StorageIndex,Dynamic,1> VectorI;
+ typedef SparseMatrix<Scalar,RowMajor,StorageIndex> FactorType;
+
+ enum {
+ ColsAtCompileTime = Dynamic,
+ MaxColsAtCompileTime = Dynamic
+ };
public:
- typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
IncompleteLUT()
: m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10),
- m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false)
+ m_analysisIsOk(false), m_factorizationIsOk(false)
{}
template<typename MatrixType>
- IncompleteLUT(const MatrixType& mat, const RealScalar& droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10)
+ explicit IncompleteLUT(const MatrixType& mat, const RealScalar& droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10)
: m_droptol(droptol),m_fillfactor(fillfactor),
- m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false)
+ m_analysisIsOk(false),m_factorizationIsOk(false)
{
eigen_assert(fillfactor != 0);
compute(mat);
@@ -146,7 +157,7 @@
*
**/
template<typename MatrixType>
- IncompleteLUT<Scalar>& compute(const MatrixType& amat)
+ IncompleteLUT& compute(const MatrixType& amat)
{
analyzePattern(amat);
factorize(amat);
@@ -157,23 +168,14 @@
void setFillfactor(int fillfactor);
template<typename Rhs, typename Dest>
- void _solve(const Rhs& b, Dest& x) const
+ void _solve_impl(const Rhs& b, Dest& x) const
{
- x = m_Pinv * b;
+ x = m_Pinv * b;
x = m_lu.template triangularView<UnitLower>().solve(x);
x = m_lu.template triangularView<Upper>().solve(x);
x = m_P * x;
}
- template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs>
- solve(const MatrixBase<Rhs>& b) const
- {
- eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
- eigen_assert(cols()==b.rows()
- && "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
- return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived());
- }
-
protected:
/** keeps off-diagonal entries; drops diagonal entries */
@@ -191,18 +193,17 @@
int m_fillfactor;
bool m_analysisIsOk;
bool m_factorizationIsOk;
- bool m_isInitialized;
ComputationInfo m_info;
- PermutationMatrix<Dynamic,Dynamic,Index> m_P; // Fill-reducing permutation
- PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // Inverse permutation
+ PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P; // Fill-reducing permutation
+ PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // Inverse permutation
};
/**
* Set control parameter droptol
* \param droptol Drop any element whose magnitude is less than this tolerance
**/
-template<typename Scalar>
-void IncompleteLUT<Scalar>::setDroptol(const RealScalar& droptol)
+template<typename Scalar, typename StorageIndex>
+void IncompleteLUT<Scalar,StorageIndex>::setDroptol(const RealScalar& droptol)
{
this->m_droptol = droptol;
}
@@ -211,52 +212,62 @@
* Set control parameter fillfactor
* \param fillfactor This is used to compute the number @p fill_in of largest elements to keep on each row.
**/
-template<typename Scalar>
-void IncompleteLUT<Scalar>::setFillfactor(int fillfactor)
+template<typename Scalar, typename StorageIndex>
+void IncompleteLUT<Scalar,StorageIndex>::setFillfactor(int fillfactor)
{
this->m_fillfactor = fillfactor;
}
-template <typename Scalar>
+template <typename Scalar, typename StorageIndex>
template<typename _MatrixType>
-void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat)
+void IncompleteLUT<Scalar,StorageIndex>::analyzePattern(const _MatrixType& amat)
{
// Compute the Fill-reducing permutation
- SparseMatrix<Scalar,ColMajor, Index> mat1 = amat;
- SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose();
- // Symmetrize the pattern
+ // Since ILUT does not perform any numerical pivoting,
+ // it is highly preferable to keep the diagonal through symmetric permutations.
+#ifndef EIGEN_MPL2_ONLY
+ // To this end, let's symmetrize the pattern and perform AMD on it.
+ SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat;
+ SparseMatrix<Scalar,ColMajor, StorageIndex> mat2 = amat.transpose();
// FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice.
// on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered...
- SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 + mat1;
- AtA.prune(keep_diag());
- internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); // Then compute the AMD ordering...
-
- m_Pinv = m_P.inverse(); // ... and the inverse permutation
+ SparseMatrix<Scalar,ColMajor, StorageIndex> AtA = mat2 + mat1;
+ AMDOrdering<StorageIndex> ordering;
+ ordering(AtA,m_P);
+ m_Pinv = m_P.inverse(); // cache the inverse permutation
+#else
+ // If AMD is not available, (MPL2-only), then let's use the slower COLAMD routine.
+ SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat;
+ COLAMDOrdering<StorageIndex> ordering;
+ ordering(mat1,m_Pinv);
+ m_P = m_Pinv.inverse();
+#endif
m_analysisIsOk = true;
m_factorizationIsOk = false;
- m_isInitialized = false;
+ m_isInitialized = true;
}
-template <typename Scalar>
+template <typename Scalar, typename StorageIndex>
template<typename _MatrixType>
-void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
+void IncompleteLUT<Scalar,StorageIndex>::factorize(const _MatrixType& amat)
{
using std::sqrt;
using std::swap;
using std::abs;
+ using internal::convert_index;
eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
Index n = amat.cols(); // Size of the matrix
m_lu.resize(n,n);
// Declare Working vectors and variables
Vector u(n) ; // real values of the row -- maximum size is n --
- VectorXi ju(n); // column position of the values in u -- maximum size is n
- VectorXi jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1
+ VectorI ju(n); // column position of the values in u -- maximum size is n
+ VectorI jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1
// Apply the fill-reducing permutation
eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
- SparseMatrix<Scalar,RowMajor, Index> mat;
+ SparseMatrix<Scalar,RowMajor, StorageIndex> mat;
mat = amat.twistedBy(m_Pinv);
// Initialization
@@ -265,7 +276,7 @@
u.fill(0);
// number of largest elements to keep in each row:
- Index fill_in = static_cast<Index> (amat.nonZeros()*m_fillfactor)/n+1;
+ Index fill_in = (amat.nonZeros()*m_fillfactor)/n + 1;
if (fill_in > n) fill_in = n;
// number of largest nonzero elements to keep in the L and the U part of the current row:
@@ -280,9 +291,9 @@
Index sizeu = 1; // number of nonzero elements in the upper part of the current row
Index sizel = 0; // number of nonzero elements in the lower part of the current row
- ju(ii) = ii;
+ ju(ii) = convert_index<StorageIndex>(ii);
u(ii) = 0;
- jr(ii) = ii;
+ jr(ii) = convert_index<StorageIndex>(ii);
RealScalar rownorm = 0;
typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii
@@ -292,9 +303,9 @@
if (k < ii)
{
// copy the lower part
- ju(sizel) = k;
+ ju(sizel) = convert_index<StorageIndex>(k);
u(sizel) = j_it.value();
- jr(k) = sizel;
+ jr(k) = convert_index<StorageIndex>(sizel);
++sizel;
}
else if (k == ii)
@@ -305,9 +316,9 @@
{
// copy the upper part
Index jpos = ii + sizeu;
- ju(jpos) = k;
+ ju(jpos) = convert_index<StorageIndex>(k);
u(jpos) = j_it.value();
- jr(k) = jpos;
+ jr(k) = convert_index<StorageIndex>(jpos);
++sizeu;
}
rownorm += numext::abs2(j_it.value());
@@ -337,7 +348,8 @@
// swap the two locations
Index j = ju(jj);
swap(ju(jj), ju(k));
- jr(minrow) = jj; jr(j) = k;
+ jr(minrow) = convert_index<StorageIndex>(jj);
+ jr(j) = convert_index<StorageIndex>(k);
swap(u(jj), u(k));
}
// Reset this location
@@ -361,8 +373,8 @@
for (; ki_it; ++ki_it)
{
Scalar prod = fact * ki_it.value();
- Index j = ki_it.index();
- Index jpos = jr(j);
+ Index j = ki_it.index();
+ Index jpos = jr(j);
if (jpos == -1) // fill-in element
{
Index newpos;
@@ -378,16 +390,16 @@
sizel++;
eigen_internal_assert(sizel<=ii);
}
- ju(newpos) = j;
+ ju(newpos) = convert_index<StorageIndex>(j);
u(newpos) = -prod;
- jr(j) = newpos;
+ jr(j) = convert_index<StorageIndex>(newpos);
}
else
u(jpos) -= prod;
}
// store the pivot element
- u(len) = fact;
- ju(len) = minrow;
+ u(len) = fact;
+ ju(len) = convert_index<StorageIndex>(minrow);
++len;
jj++;
@@ -402,7 +414,7 @@
sizel = len;
len = (std::min)(sizel, nnzL);
typename Vector::SegmentReturnType ul(u.segment(0, sizel));
- typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel));
+ typename VectorI::SegmentReturnType jul(ju.segment(0, sizel));
internal::QuickSplit(ul, jul, len);
// store the largest m_fill elements of the L part
@@ -431,39 +443,20 @@
sizeu = len + 1; // +1 to take into account the diagonal element
len = (std::min)(sizeu, nnzU);
typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
- typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
+ typename VectorI::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
internal::QuickSplit(uu, juu, len);
// store the largest elements of the U part
for(Index k = ii + 1; k < ii + len; k++)
m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
}
-
m_lu.finalize();
m_lu.makeCompressed();
m_factorizationIsOk = true;
- m_isInitialized = m_factorizationIsOk;
m_info = Success;
}
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<IncompleteLUT<_MatrixType>, Rhs>
- : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs>
-{
- typedef IncompleteLUT<_MatrixType> Dec;
- EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- dec()._solve(rhs(),dst);
- }
-};
-
-} // end namespace internal
-
} // end namespace Eigen
#endif // EIGEN_INCOMPLETE_LUT_H