Incomplete Cholesky preconditioner... not yet stable
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 224304f..5a71531 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
@@ -10,8 +10,56 @@ #ifndef EIGEN_INCOMPLETE_LUT_H #define EIGEN_INCOMPLETE_LUT_H + namespace Eigen { +namespace internal { + +/** + * Compute a quick-sort split of a vector + * On output, the vector row is permuted such that its elements satisfy + * abs(row(i)) >= abs(row(ncut)) if i<ncut + * abs(row(i)) <= abs(row(ncut)) if i>ncut + * \param row The vector of values + * \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> +int QuickSplit(VectorV &row, VectorI &ind, int ncut) +{ + typedef typename VectorV::RealScalar RealScalar; + using std::swap; + int mid; + int n = row.size(); /* length of the vector */ + int first, last ; + + ncut--; /* to fit the zero-based indices */ + first = 0; + last = n-1; + if (ncut < first || ncut > last ) return 0; + + do { + mid = first; + RealScalar abskey = std::abs(row(mid)); + for (int j = first + 1; j <= last; j++) { + if ( std::abs(row(j)) > abskey) { + ++mid; + swap(row(mid), row(j)); + swap(ind(mid), ind(j)); + } + } + /* Interchange for the pivot element */ + swap(row(mid), row(first)); + swap(ind(mid), ind(first)); + + if (mid > ncut) last = mid - 1; + else if (mid < ncut ) first = mid + 1; + } while (mid != ncut ); + + return 0; /* mid is equal to ncut */ +} + +}// end namespace internal /** * \brief Incomplete LU factorization with dual-threshold strategy * During the numerical factorization, two dropping rules are used : @@ -126,10 +174,6 @@ protected: - template <typename VectorV, typename VectorI> - int QuickSplit(VectorV &row, VectorI &ind, int ncut); - - /** keeps off-diagonal entries; drops diagonal entries */ struct keep_diag { inline bool operator() (const Index& row, const Index& col, const Scalar&) const @@ -171,51 +215,6 @@ this->m_fillfactor = fillfactor; } - -/** - * Compute a quick-sort split of a vector - * On output, the vector row is permuted such that its elements satisfy - * abs(row(i)) >= abs(row(ncut)) if i<ncut - * abs(row(i)) <= abs(row(ncut)) if i>ncut - * \param row The vector of values - * \param ind The array of index for the elements in @p row - * \param ncut The number of largest elements to keep - **/ -template <typename Scalar> -template <typename VectorV, typename VectorI> -int IncompleteLUT<Scalar>::QuickSplit(VectorV &row, VectorI &ind, int ncut) -{ - using std::swap; - int mid; - int n = row.size(); /* length of the vector */ - int first, last ; - - ncut--; /* to fit the zero-based indices */ - first = 0; - last = n-1; - if (ncut < first || ncut > last ) return 0; - - do { - mid = first; - RealScalar abskey = std::abs(row(mid)); - for (int j = first + 1; j <= last; j++) { - if ( std::abs(row(j)) > abskey) { - ++mid; - swap(row(mid), row(j)); - swap(ind(mid), ind(j)); - } - } - /* Interchange for the pivot element */ - swap(row(mid), row(first)); - swap(ind(mid), ind(first)); - - if (mid > ncut) last = mid - 1; - else if (mid < ncut ) first = mid + 1; - } while (mid != ncut ); - - return 0; /* mid is equal to ncut */ -} - template <typename Scalar> template<typename _MatrixType> void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat) @@ -400,7 +399,7 @@ len = (std::min)(sizel, nnzL); typename Vector::SegmentReturnType ul(u.segment(0, sizel)); typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel)); - QuickSplit(ul, jul, len); + internal::QuickSplit(ul, jul, len); // store the largest m_fill elements of the L part m_lu.startVec(ii); @@ -429,7 +428,7 @@ 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)); - QuickSplit(uu, juu, len); + internal::QuickSplit(uu, juu, len); // store the largest elements of the U part for(int k = ii + 1; k < ii + len; k++)
diff --git a/bench/spbench/sp_solver.cpp b/bench/spbench/sp_solver.cpp index e18f2d1..a1f4bac 100644 --- a/bench/spbench/sp_solver.cpp +++ b/bench/spbench/sp_solver.cpp
@@ -13,7 +13,7 @@ #include <Eigen/SuperLUSupport> // #include <unsupported/Eigen/src/IterativeSolvers/Scaling.h> #include <bench/BenchTimer.h> - +#include <unsupported/Eigen/IterativeSolvers> using namespace std; using namespace Eigen; @@ -26,7 +26,8 @@ VectorXd b, x, tmp; BenchTimer timer,totaltime; //SparseLU<SparseMatrix<double, ColMajor> > solver; - SuperLU<SparseMatrix<double, ColMajor> > solver; +// SuperLU<SparseMatrix<double, ColMajor> > solver; + ConjugateGradient<SparseMatrix<double, ColMajor>, Lower,IncompleteCholesky<double,Lower> > solver; ifstream matrix_file; string line; int n;
diff --git a/unsupported/Eigen/IterativeSolvers b/unsupported/Eigen/IterativeSolvers index 6c6946d..c3cc97c 100644 --- a/unsupported/Eigen/IterativeSolvers +++ b/unsupported/Eigen/IterativeSolvers
@@ -33,6 +33,7 @@ #include "../../Eigen/Jacobi" #include "../../Eigen/Householder" #include "src/IterativeSolvers/GMRES.h" +#include "src/IterativeSolvers/IncompleteCholesky.h" //#include "src/IterativeSolvers/SSORPreconditioner.h" //@}
diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h new file mode 100644 index 0000000..bdd494f --- /dev/null +++ b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
@@ -0,0 +1,221 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// 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_INCOMPLETE_CHOlESKY_H +#define EIGEN_INCOMPLETE_CHOlESKY_H +#include "Eigen/src/IterativeLinearSolvers/IncompleteLUT.h" +#include <Eigen/OrderingMethods> +#include <list> + +namespace Eigen { +/** + * \brief Modified Incomplete Cholesky with dual threshold + * + * References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with + * Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999 + * + * \tparam _MatrixType The type of the sparse matrix. It should be a symmetric + * matrix. It is advised to give a row-oriented sparse matrix + * \tparam _UpLo The triangular part of the matrix to reference. + * \tparam _OrderingType + */ + +template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> > +class IncompleteCholesky : internal::noncopyable +{ + public: + typedef SparseMatrix<Scalar,ColMajor> MatrixType; + typedef _OrderingType OrderingType; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; + typedef Matrix<Scalar,Dynamic,1> VectorType; + typedef Matrix<Index,Dynamic, 1> IndexType; + + public: + IncompleteCholesky() {} + IncompleteCholesky(const MatrixType& matrix) + { + compute(matrix); + } + + Index rows() const { return m_L.rows(); } + + Index cols() const { return m_L.cols(); } + + + /** \brief Reports whether previous computation was successful. + * + * \returns \c Success if computation was succesful, + * \c NumericalIssue if the matrix appears to be negative. + */ + ComputationInfo info() const + { + eigen_assert(m_isInitialized && "IncompleteLLT is not initialized."); + return m_info; + } + /** + * \brief Computes the fill reducing permutation vector. + */ + template<typename MatrixType> + void analyzePattern(const MatrixType& mat) + { + OrderingType ord; + ord(mat, m_perm); + m_analysisIsOk = true; + } + + template<typename MatrixType> + void factorize(const MatrixType& amat); + + template<typename MatrixType> + void compute (const MatrixType& matrix) + { + analyzePattern(matrix); + factorize(matrix); + } + + template<typename Rhs, typename Dest> + void _solve(const Rhs& b, Dest& x) const + { + eigen_assert(m_factorizationIsOk && "factorize() should be called first"); + if (m_perm.rows() == b.rows()) + x = m_perm.inverse() * b; + else + x = b; + x = m_L.template triangularView<UnitLower>().solve(x); + x = m_L.adjoint().template triangularView<Upper>().solve(x); + if (m_perm.rows() == b.rows()) + x = m_perm * x; + } + template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs> + solve(const MatrixBase<Rhs>& b) const + { + eigen_assert(m_isInitialized && "IncompleteLLT is not initialized."); + eigen_assert(cols()==b.rows() + && "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b"); + return internal::solve_retval<IncompleteCholesky, Rhs>(*this, b.derived()); + } + protected: + SparseMatrix<Scalar,ColMajor> m_L; // The lower part stored in CSC + bool m_analysisIsOk; + bool m_factorizationIsOk; + bool m_isInitialized; + ComputationInfo m_info; + PermutationType m_perm; + +}; + +template<typename Scalar, int _UpLo, typename OrderingType> +template<typename _MatrixType> +void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat) +{ + eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); + + // FIXME Stability: We should probably compute the scaling factors and the shifts that are needed to ensure an efficient LLT preconditioner. + + // Dropping strategies : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added + + // Apply the fill-reducing permutation computed in analyzePattern() + if (m_perm.rows() == mat.rows() ) + m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>().twistedBy(m_perm); + else + m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>(); + + int n = mat.cols(); + + Scalar *vals = m_L.valuePtr(); //Values + Index *rowIdx = m_L.innerIndexPtr(); //Row indices + Index *colPtr = m_L.outerIndexPtr(); // Pointer to the beginning of each row + VectorType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization + // Initialize firstElt; + for (int j = 0; j < n-1; j++) firstElt(j) = colPtr[j]+1; + std::vector<std::list<Index> > listCol(n); // listCol(j) is a linked list of columns to update column j + VectorType curCol(n); // Store a nonzero values in each column + VectorType irow(n); // Row indices of nonzero elements in each column + // jki version of the Cholesky factorization + for (int j=0; j < n; j++) + { + //Left-looking factorize the column j + // First, load the jth column into curCol + Scalar diag = vals[colPtr[j]]; // Lower diagonal matrix with + curCol.setZero(); + irow.setLinSpaced(n,0,n-1); + for (int i = colPtr[j] + 1; i < colPtr[j+1]; i++) + { + curCol(rowIdx[i]) = vals[i]; + irow(rowIdx[i]) = rowIdx[i]; + } + + std::list<int>::iterator k; + // Browse all previous columns that will update column j + for(k = listCol[j].begin(); k != listCol[j].end(); k++) + { + int jk = firstElt(*k); // First element to use in the column + Scalar a_jk = vals[jk]; + diag -= a_jk * a_jk; + jk += 1; + for (int i = jk; i < colPtr[*k]; i++) + { + curCol(rowIdx[i]) -= vals[i] * a_jk ; + } + firstElt(*k) = jk; + if (jk < colPtr[*k+1]) + { + // Add this column to the updating columns list for column *k+1 + listCol[rowIdx[jk]].push_back(*k); + } + } + + // Select the largest p elements + // p is the original number of elements in the column (without the diagonal) + int p = colPtr[j+1] - colPtr[j] - 2 ; + internal::QuickSplit(curCol, irow, p); + if(RealScalar(diag) <= 0) + { + m_info = NumericalIssue; + return; + } + RealScalar rdiag = internal::sqrt(RealScalar(diag)); + Scalar scal = Scalar(1)/rdiag; + vals[colPtr[j]] = rdiag; + // Insert the largest p elements in the matrix and scale them meanwhile + int cpt = 0; + for (int i = colPtr[j]+1; i < colPtr[j+1]; i++) + { + vals[i] = curCol(cpt) * scal; + rowIdx[i] = irow(cpt); + cpt ++; + } + } + m_factorizationIsOk = true; + m_isInitialized = true; + m_info = Success; +} + +namespace internal { + +template<typename _MatrixType, typename Rhs> +struct solve_retval<IncompleteCholesky<_MatrixType>, Rhs> + : solve_retval_base<IncompleteCholesky<_MatrixType>, Rhs> +{ + typedef IncompleteCholesky<_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 \ No newline at end of file