diff options
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h | 11 | ||||
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/IncompleteLUT.h | 21 | ||||
| -rw-r--r-- | test/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | test/mpl2only.cpp | 20 |
4 files changed, 45 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; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9684c90e8..c8a8ba6f4 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -255,6 +255,7 @@ ei_add_test(special_numbers) ei_add_test(rvalue_types) ei_add_test(dense_storage) ei_add_test(ctorleak) +ei_add_test(mpl2only) # # ei_add_test(denseLM) diff --git a/test/mpl2only.cpp b/test/mpl2only.cpp new file mode 100644 index 000000000..5ef0d2b2e --- /dev/null +++ b/test/mpl2only.cpp @@ -0,0 +1,20 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2015 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 +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#define EIGEN_MPL2_ONLY +#include <Eigen/Dense> +#include <Eigen/SparseCore> +#include <Eigen/SparseLU> +#include <Eigen/SparseQR> +#include <Eigen/IterativeLinearSolvers> + +int main() +{ + return 0; +} |