diff options
Diffstat (limited to 'unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h')
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h | 278 |
1 files changed, 0 insertions, 278 deletions
diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h deleted file mode 100644 index 661c1f2e0..000000000 --- a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h +++ /dev/null @@ -1,278 +0,0 @@ -// 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> ScalarType; - typedef Matrix<Index,Dynamic, 1> IndexType; - typedef std::vector<std::list<Index> > VectorList; - enum { UpLo = _UpLo }; - public: - IncompleteCholesky() : m_shift(1),m_factorizationIsOk(false) {} - IncompleteCholesky(const MatrixType& matrix) : m_shift(1),m_factorizationIsOk(false) - { - 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 Set the initial shift parameter - */ - void setShift( Scalar shift) { m_shift = shift; } - - /** - * \brief Computes the fill reducing permutation vector. - */ - template<typename MatrixType> - void analyzePattern(const MatrixType& mat) - { - OrderingType ord; - ord(mat.template selfadjointView<UpLo>(), 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_scal.asDiagonal() * x; - 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; - x = m_scal.asDiagonal() * x; - } - template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_factorizationIsOk && "IncompleteLLT did not succeed"); - 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 - ScalarType m_scal; // The vector for scaling the matrix - Scalar m_shift; //The initial shift parameter - bool m_analysisIsOk; - bool m_factorizationIsOk; - bool m_isInitialized; - ComputationInfo m_info; - PermutationType m_perm; - - private: - template <typename IdxType, typename SclType> - inline void updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol); -}; - -template<typename Scalar, int _UpLo, typename OrderingType> -template<typename _MatrixType> -void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat) -{ - using std::sqrt; - using std::min; - eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); - - // 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() ) // To detect the null permutation - m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>().twistedBy(m_perm); - else - m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>(); - - Index n = m_L.cols(); - Index nnz = m_L.nonZeros(); - Map<ScalarType> vals(m_L.valuePtr(), nnz); //values - Map<IndexType> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices - Map<IndexType> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row - IndexType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization - VectorList listCol(n); // listCol(j) is a linked list of columns to update column j - ScalarType curCol(n); // Store a nonzero values in each column - IndexType irow(n); // Row indices of nonzero elements in each column - - - // Computes the scaling factors - m_scal.resize(n); - for (int j = 0; j < n; j++) - { - m_scal(j) = m_L.col(j).norm(); - m_scal(j) = sqrt(m_scal(j)); - } - // Scale and compute the shift for the matrix - Scalar mindiag = vals[0]; - for (int j = 0; j < n; j++){ - for (int k = colPtr[j]; k < colPtr[j+1]; k++) - vals[k] /= (m_scal(j) * m_scal(rowIdx[k])); - mindiag = (min)(vals[colPtr[j]], mindiag); - } - - if(mindiag < Scalar(0.)) m_shift = m_shift - mindiag; - // Apply the shift to the diagonal elements of the matrix - for (int j = 0; j < n; j++) - vals[colPtr[j]] += m_shift; - // 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]]; // It is assumed that only the lower part is stored - 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 - jk += 1; - for (int i = jk; i < colPtr[*k+1]; i++) - { - curCol(rowIdx[i]) -= vals[i] * vals[jk] ; - } - updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol); - } - - // Scale the current column - if(RealScalar(diag) <= 0) - { - std::cerr << "\nNegative diagonal during Incomplete factorization... "<< j << "\n"; - m_info = NumericalIssue; - return; - } - RealScalar rdiag = sqrt(RealScalar(diag)); - vals[colPtr[j]] = rdiag; - for (int i = j+1; i < n; i++) - { - //Scale - curCol(i) /= rdiag; - //Update the remaining diagonals with curCol - vals[colPtr[i]] -= curCol(i) * curCol(i); - } - // 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] - 1 ; - internal::QuickSplit(curCol, irow, p); - // Insert the largest p elements in the matrix - int cpt = 0; - for (int i = colPtr[j]+1; i < colPtr[j+1]; i++) - { - vals[i] = curCol(cpt); - rowIdx[i] = irow(cpt); - cpt ++; - } - // Get the first smallest row index and put it after the diagonal element - Index jk = colPtr(j)+1; - updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); - } - m_factorizationIsOk = true; - m_isInitialized = true; - m_info = Success; -} - -template<typename Scalar, int _UpLo, typename OrderingType> -template <typename IdxType, typename SclType> -inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol) -{ - if (jk < colPtr(col+1) ) - { - Index p = colPtr(col+1) - jk; - Index minpos; - rowIdx.segment(jk,p).minCoeff(&minpos); - minpos += jk; - if (rowIdx(minpos) != rowIdx(jk)) - { - //Swap - std::swap(rowIdx(jk),rowIdx(minpos)); - std::swap(vals(jk),vals(minpos)); - } - firstElt(col) = jk; - listCol[rowIdx(jk)].push_back(col); - } -} -namespace internal { - -template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs> -struct solve_retval<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs> - : solve_retval_base<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs> -{ - typedef IncompleteCholesky<_Scalar, _UpLo, OrderingType> 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 |