1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H 11 #define EIGEN_INCOMPLETE_CHOlESKY_H 12 #include "Eigen/src/IterativeLinearSolvers/IncompleteLUT.h" 13 #include <Eigen/OrderingMethods> 14 #include <list> 15 16 namespace Eigen { 17 /** 18 * \brief Modified Incomplete Cholesky with dual threshold 19 * 20 * References : C-J. Lin and J. J. Mor, Incomplete Cholesky Factorizations with 21 * Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999 22 * 23 * \tparam _MatrixType The type of the sparse matrix. It should be a symmetric 24 * matrix. It is advised to give a row-oriented sparse matrix 25 * \tparam _UpLo The triangular part of the matrix to reference. 26 * \tparam _OrderingType 27 */ 28 29 template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> > 30 class IncompleteCholesky : internal::noncopyable 31 { 32 public: 33 typedef SparseMatrix<Scalar,ColMajor> MatrixType; 34 typedef _OrderingType OrderingType; 35 typedef typename MatrixType::RealScalar RealScalar; 36 typedef typename MatrixType::Index Index; 37 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 38 typedef Matrix<Scalar,Dynamic,1> ScalarType; 39 typedef Matrix<Index,Dynamic, 1> IndexType; 40 typedef std::vector<std::list<Index> > VectorList; 41 enum { UpLo = _UpLo }; 42 public: 43 IncompleteCholesky() : m_shift(1),m_factorizationIsOk(false) {} 44 IncompleteCholesky(const MatrixType& matrix) : m_shift(1),m_factorizationIsOk(false) 45 { 46 compute(matrix); 47 } 48 49 Index rows() const { return m_L.rows(); } 50 51 Index cols() const { return m_L.cols(); } 52 53 54 /** \brief Reports whether previous computation was successful. 55 * 56 * \returns \c Success if computation was succesful, 57 * \c NumericalIssue if the matrix appears to be negative. 58 */ 59 ComputationInfo info() const 60 { 61 eigen_assert(m_isInitialized && "IncompleteLLT is not initialized."); 62 return m_info; 63 } 64 65 /** 66 * \brief Set the initial shift parameter 67 */ 68 void setShift( Scalar shift) { m_shift = shift; } 69 70 /** 71 * \brief Computes the fill reducing permutation vector. 72 */ 73 template<typename MatrixType> 74 void analyzePattern(const MatrixType& mat) 75 { 76 OrderingType ord; 77 ord(mat.template selfadjointView<UpLo>(), m_perm); 78 m_analysisIsOk = true; 79 } 80 81 template<typename MatrixType> 82 void factorize(const MatrixType& amat); 83 84 template<typename MatrixType> 85 void compute (const MatrixType& matrix) 86 { 87 analyzePattern(matrix); 88 factorize(matrix); 89 } 90 91 template<typename Rhs, typename Dest> 92 void _solve(const Rhs& b, Dest& x) const 93 { 94 eigen_assert(m_factorizationIsOk && "factorize() should be called first"); 95 if (m_perm.rows() == b.rows()) 96 x = m_perm.inverse() * b; 97 else 98 x = b; 99 x = m_scal.asDiagonal() * x; 100 x = m_L.template triangularView<UnitLower>().solve(x); 101 x = m_L.adjoint().template triangularView<Upper>().solve(x); 102 if (m_perm.rows() == b.rows()) 103 x = m_perm * x; 104 x = m_scal.asDiagonal() * x; 105 } 106 template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs> 107 solve(const MatrixBase<Rhs>& b) const 108 { 109 eigen_assert(m_factorizationIsOk && "IncompleteLLT did not succeed"); 110 eigen_assert(m_isInitialized && "IncompleteLLT is not initialized."); 111 eigen_assert(cols()==b.rows() 112 && "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b"); 113 return internal::solve_retval<IncompleteCholesky, Rhs>(*this, b.derived()); 114 } 115 protected: 116 SparseMatrix<Scalar,ColMajor> m_L; // The lower part stored in CSC 117 ScalarType m_scal; // The vector for scaling the matrix 118 Scalar m_shift; //The initial shift parameter 119 bool m_analysisIsOk; 120 bool m_factorizationIsOk; 121 bool m_isInitialized; 122 ComputationInfo m_info; 123 PermutationType m_perm; 124 125 private: 126 template <typename IdxType, typename SclType> 127 inline void updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol); 128 }; 129 130 template<typename Scalar, int _UpLo, typename OrderingType> 131 template<typename _MatrixType> 132 void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat) 133 { 134 using std::sqrt; 135 using std::min; 136 eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 137 138 // 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 139 140 // Apply the fill-reducing permutation computed in analyzePattern() 141 if (m_perm.rows() == mat.rows() ) // To detect the null permutation 142 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>().twistedBy(m_perm); 143 else 144 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>(); 145 146 Index n = m_L.cols(); 147 Index nnz = m_L.nonZeros(); 148 Map<ScalarType> vals(m_L.valuePtr(), nnz); //values 149 Map<IndexType> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices 150 Map<IndexType> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row 151 IndexType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization 152 VectorList listCol(n); // listCol(j) is a linked list of columns to update column j 153 ScalarType curCol(n); // Store a nonzero values in each column 154 IndexType irow(n); // Row indices of nonzero elements in each column 155 156 157 // Computes the scaling factors 158 m_scal.resize(n); 159 for (int j = 0; j < n; j++) 160 { 161 m_scal(j) = m_L.col(j).norm(); 162 m_scal(j) = sqrt(m_scal(j)); 163 } 164 // Scale and compute the shift for the matrix 165 Scalar mindiag = vals[0]; 166 for (int j = 0; j < n; j++){ 167 for (int k = colPtr[j]; k < colPtr[j+1]; k++) 168 vals[k] /= (m_scal(j) * m_scal(rowIdx[k])); 169 mindiag = (min)(vals[colPtr[j]], mindiag); 170 } 171 172 if(mindiag < Scalar(0.)) m_shift = m_shift - mindiag; 173 // Apply the shift to the diagonal elements of the matrix 174 for (int j = 0; j < n; j++) 175 vals[colPtr[j]] += m_shift; 176 // jki version of the Cholesky factorization 177 for (int j=0; j < n; ++j) 178 { 179 //Left-looking factorize the column j 180 // First, load the jth column into curCol 181 Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored 182 curCol.setZero(); 183 irow.setLinSpaced(n,0,n-1); 184 for (int i = colPtr[j] + 1; i < colPtr[j+1]; i++) 185 { 186 curCol(rowIdx[i]) = vals[i]; 187 irow(rowIdx[i]) = rowIdx[i]; 188 } 189 std::list<int>::iterator k; 190 // Browse all previous columns that will update column j 191 for(k = listCol[j].begin(); k != listCol[j].end(); k++) 192 { 193 int jk = firstElt(*k); // First element to use in the column 194 jk += 1; 195 for (int i = jk; i < colPtr[*k+1]; i++) 196 { 197 curCol(rowIdx[i]) -= vals[i] * vals[jk] ; 198 } 199 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol); 200 } 201 202 // Scale the current column 203 if(RealScalar(diag) <= 0) 204 { 205 std::cerr << "\nNegative diagonal during Incomplete factorization... "<< j << "\n"; 206 m_info = NumericalIssue; 207 return; 208 } 209 RealScalar rdiag = sqrt(RealScalar(diag)); 210 vals[colPtr[j]] = rdiag; 211 for (int i = j+1; i < n; i++) 212 { 213 //Scale 214 curCol(i) /= rdiag; 215 //Update the remaining diagonals with curCol 216 vals[colPtr[i]] -= curCol(i) * curCol(i); 217 } 218 // Select the largest p elements 219 // p is the original number of elements in the column (without the diagonal) 220 int p = colPtr[j+1] - colPtr[j] - 1 ; 221 internal::QuickSplit(curCol, irow, p); 222 // Insert the largest p elements in the matrix 223 int cpt = 0; 224 for (int i = colPtr[j]+1; i < colPtr[j+1]; i++) 225 { 226 vals[i] = curCol(cpt); 227 rowIdx[i] = irow(cpt); 228 cpt ++; 229 } 230 // Get the first smallest row index and put it after the diagonal element 231 Index jk = colPtr(j)+1; 232 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); 233 } 234 m_factorizationIsOk = true; 235 m_isInitialized = true; 236 m_info = Success; 237 } 238 239 template<typename Scalar, int _UpLo, typename OrderingType> 240 template <typename IdxType, typename SclType> 241 inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol) 242 { 243 if (jk < colPtr(col+1) ) 244 { 245 Index p = colPtr(col+1) - jk; 246 Index minpos; 247 rowIdx.segment(jk,p).minCoeff(&minpos); 248 minpos += jk; 249 if (rowIdx(minpos) != rowIdx(jk)) 250 { 251 //Swap 252 std::swap(rowIdx(jk),rowIdx(minpos)); 253 std::swap(vals(jk),vals(minpos)); 254 } 255 firstElt(col) = jk; 256 listCol[rowIdx(jk)].push_back(col); 257 } 258 } 259 namespace internal { 260 261 template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs> 262 struct solve_retval<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs> 263 : solve_retval_base<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs> 264 { 265 typedef IncompleteCholesky<_Scalar, _UpLo, OrderingType> Dec; 266 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 267 268 template<typename Dest> void evalTo(Dest& dst) const 269 { 270 dec()._solve(rhs(),dst); 271 } 272 }; 273 274 } // end namespace internal 275 276 } // end namespace Eigen 277 278 #endif 279
