diff options
Diffstat (limited to 'Eigen/src/IterativeLinearSolvers')
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h | 11 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/IncompleteLUT.h | 21 |
2 files changed, 24 insertions, 8 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h index 1e2e9f9b9..8f549af82 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h @@ -24,7 +24,8 @@ namespace Eigen { * \tparam _MatrixType The type of the sparse matrix. It is advised to give a row-oriented sparse matrix * \tparam _UpLo The triangular part that will be used for the computations. It can be Lower * or Upper. Default is Lower. - * \tparam _OrderingType The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> + * \tparam _OrderingType The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<int>, + * unless EIGEN_MPL2_ONLY is defined, in which case the default is NaturalOrdering<int>. * * \implsparsesolverconcept * @@ -38,7 +39,13 @@ namespace Eigen { * \f$ \sigma \f$ is the initial shift value as returned and set by setInitialShift() method. The default value is \f$ \sigma = 10^{-3} \f$. * */ -template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> > +template <typename Scalar, int _UpLo = Lower, typename _OrderingType = +#ifndef EIGEN_MPL2_ONLY +AMDOrdering<int> +#else +NaturalOrdering<int> +#endif +> class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > { protected: diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 10b9fcc18..519472377 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -168,7 +168,7 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar, _StorageInd template<typename Rhs, typename Dest> 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; @@ -221,16 +221,25 @@ template<typename _MatrixType> void IncompleteLUT<Scalar,StorageIndex>::analyzePattern(const _MatrixType& amat) { // Compute the Fill-reducing permutation + // 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(); - // Symmetrize the pattern // 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, StorageIndex> AtA = mat2 + mat1; - AtA.prune(keep_diag()); - internal::minimum_degree_ordering<Scalar, StorageIndex>(AtA, m_P); // Then compute the AMD ordering... - - m_Pinv = m_P.inverse(); // ... and the inverse permutation + 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; |