![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 00005 // Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H 00012 #define EIGEN_INCOMPLETE_CHOlESKY_H 00013 00014 #include <vector> 00015 #include <list> 00016 00017 namespace Eigen { 00044 template <typename Scalar, int _UpLo = Lower, typename _OrderingType = 00045 #ifndef EIGEN_MPL2_ONLY 00046 AMDOrdering<int> 00047 #else 00048 NaturalOrdering<int> 00049 #endif 00050 > 00051 class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > 00052 { 00053 protected: 00054 typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base; 00055 using Base::m_isInitialized; 00056 public: 00057 typedef typename NumTraits<Scalar>::Real RealScalar; 00058 typedef _OrderingType OrderingType; 00059 typedef typename OrderingType::PermutationType PermutationType; 00060 typedef typename PermutationType::StorageIndex StorageIndex; 00061 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType; 00062 typedef Matrix<Scalar,Dynamic,1> VectorSx; 00063 typedef Matrix<RealScalar,Dynamic,1> VectorRx; 00064 typedef Matrix<StorageIndex,Dynamic, 1> VectorIx; 00065 typedef std::vector<std::list<StorageIndex> > VectorList; 00066 enum { UpLo = _UpLo }; 00067 enum { 00068 ColsAtCompileTime = Dynamic, 00069 MaxColsAtCompileTime = Dynamic 00070 }; 00071 public: 00072 00079 IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {} 00080 00083 template<typename MatrixType> 00084 IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false) 00085 { 00086 compute(matrix); 00087 } 00088 00090 Index rows() const { return m_L.rows(); } 00091 00093 Index cols() const { return m_L.cols(); } 00094 00095 00104 ComputationInfo info() const 00105 { 00106 eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized."); 00107 return m_info; 00108 } 00109 00112 void setInitialShift(RealScalar shift) { m_initialShift = shift; } 00113 00116 template<typename MatrixType> 00117 void analyzePattern(const MatrixType& mat) 00118 { 00119 OrderingType ord; 00120 PermutationType pinv; 00121 ord(mat.template selfadjointView<UpLo>(), pinv); 00122 if(pinv.size()>0) m_perm = pinv.inverse(); 00123 else m_perm.resize(0); 00124 m_L.resize(mat.rows(), mat.cols()); 00125 m_analysisIsOk = true; 00126 m_isInitialized = true; 00127 m_info = Success; 00128 } 00129 00137 template<typename MatrixType> 00138 void factorize(const MatrixType& mat); 00139 00146 template<typename MatrixType> 00147 void compute(const MatrixType& mat) 00148 { 00149 analyzePattern(mat); 00150 factorize(mat); 00151 } 00152 00153 // internal 00154 template<typename Rhs, typename Dest> 00155 void _solve_impl(const Rhs& b, Dest& x) const 00156 { 00157 eigen_assert(m_factorizationIsOk && "factorize() should be called first"); 00158 if (m_perm.rows() == b.rows()) x = m_perm * b; 00159 else x = b; 00160 x = m_scale.asDiagonal() * x; 00161 x = m_L.template triangularView<Lower>().solve(x); 00162 x = m_L.adjoint().template triangularView<Upper>().solve(x); 00163 x = m_scale.asDiagonal() * x; 00164 if (m_perm.rows() == b.rows()) 00165 x = m_perm.inverse() * x; 00166 } 00167 00169 const FactorType& matrixL() const { eigen_assert("m_factorizationIsOk"); return m_L; } 00170 00172 const VectorRx& scalingS() const { eigen_assert("m_factorizationIsOk"); return m_scale; } 00173 00175 const PermutationType& permutationP() const { eigen_assert("m_analysisIsOk"); return m_perm; } 00176 00177 protected: 00178 FactorType m_L; // The lower part stored in CSC 00179 VectorRx m_scale; // The vector for scaling the matrix 00180 RealScalar m_initialShift; // The initial shift parameter 00181 bool m_analysisIsOk; 00182 bool m_factorizationIsOk; 00183 ComputationInfo m_info; 00184 PermutationType m_perm; 00185 00186 private: 00187 inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol); 00188 }; 00189 00190 // Based on the following paper: 00191 // C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with 00192 // Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999 00193 // http://ftp.mcs.anl.gov/pub/tech_reports/reports/P682.pdf 00194 template<typename Scalar, int _UpLo, typename OrderingType> 00195 template<typename _MatrixType> 00196 void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat) 00197 { 00198 using std::sqrt; 00199 eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 00200 00201 // Dropping strategy : 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 00202 00203 // Apply the fill-reducing permutation computed in analyzePattern() 00204 if (m_perm.rows() == mat.rows() ) // To detect the null permutation 00205 { 00206 // The temporary is needed to make sure that the diagonal entry is properly sorted 00207 FactorType tmp(mat.rows(), mat.cols()); 00208 tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm); 00209 m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>(); 00210 } 00211 else 00212 { 00213 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>(); 00214 } 00215 00216 Index n = m_L.cols(); 00217 Index nnz = m_L.nonZeros(); 00218 Map<VectorSx> vals(m_L.valuePtr(), nnz); //values 00219 Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices 00220 Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row 00221 VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization 00222 VectorList listCol(n); // listCol(j) is a linked list of columns to update column j 00223 VectorSx col_vals(n); // Store a nonzero values in each column 00224 VectorIx col_irow(n); // Row indices of nonzero elements in each column 00225 VectorIx col_pattern(n); 00226 col_pattern.fill(-1); 00227 StorageIndex col_nnz; 00228 00229 00230 // Computes the scaling factors 00231 m_scale.resize(n); 00232 m_scale.setZero(); 00233 for (Index j = 0; j < n; j++) 00234 for (Index k = colPtr[j]; k < colPtr[j+1]; k++) 00235 { 00236 m_scale(j) += numext::abs2(vals(k)); 00237 if(rowIdx[k]!=j) 00238 m_scale(rowIdx[k]) += numext::abs2(vals(k)); 00239 } 00240 00241 m_scale = m_scale.cwiseSqrt().cwiseSqrt(); 00242 00243 for (Index j = 0; j < n; ++j) 00244 if(m_scale(j)>(std::numeric_limits<RealScalar>::min)()) 00245 m_scale(j) = RealScalar(1)/m_scale(j); 00246 else 00247 m_scale(j) = 1; 00248 00249 // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster) 00250 00251 // Scale and compute the shift for the matrix 00252 RealScalar mindiag = NumTraits<RealScalar>::highest(); 00253 for (Index j = 0; j < n; j++) 00254 { 00255 for (Index k = colPtr[j]; k < colPtr[j+1]; k++) 00256 vals[k] *= (m_scale(j)*m_scale(rowIdx[k])); 00257 eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored"); 00258 mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag); 00259 } 00260 00261 FactorType L_save = m_L; 00262 00263 RealScalar shift = 0; 00264 if(mindiag <= RealScalar(0.)) 00265 shift = m_initialShift - mindiag; 00266 00267 m_info = NumericalIssue; 00268 00269 // Try to perform the incomplete factorization using the current shift 00270 int iter = 0; 00271 do 00272 { 00273 // Apply the shift to the diagonal elements of the matrix 00274 for (Index j = 0; j < n; j++) 00275 vals[colPtr[j]] += shift; 00276 00277 // jki version of the Cholesky factorization 00278 Index j=0; 00279 for (; j < n; ++j) 00280 { 00281 // Left-looking factorization of the j-th column 00282 // First, load the j-th column into col_vals 00283 Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored 00284 col_nnz = 0; 00285 for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++) 00286 { 00287 StorageIndex l = rowIdx[i]; 00288 col_vals(col_nnz) = vals[i]; 00289 col_irow(col_nnz) = l; 00290 col_pattern(l) = col_nnz; 00291 col_nnz++; 00292 } 00293 { 00294 typename std::list<StorageIndex>::iterator k; 00295 // Browse all previous columns that will update column j 00296 for(k = listCol[j].begin(); k != listCol[j].end(); k++) 00297 { 00298 Index jk = firstElt(*k); // First element to use in the column 00299 eigen_internal_assert(rowIdx[jk]==j); 00300 Scalar v_j_jk = numext::conj(vals[jk]); 00301 00302 jk += 1; 00303 for (Index i = jk; i < colPtr[*k+1]; i++) 00304 { 00305 StorageIndex l = rowIdx[i]; 00306 if(col_pattern[l]<0) 00307 { 00308 col_vals(col_nnz) = vals[i] * v_j_jk; 00309 col_irow[col_nnz] = l; 00310 col_pattern(l) = col_nnz; 00311 col_nnz++; 00312 } 00313 else 00314 col_vals(col_pattern[l]) -= vals[i] * v_j_jk; 00315 } 00316 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol); 00317 } 00318 } 00319 00320 // Scale the current column 00321 if(numext::real(diag) <= 0) 00322 { 00323 if(++iter>=10) 00324 return; 00325 00326 // increase shift 00327 shift = numext::maxi(m_initialShift,RealScalar(2)*shift); 00328 // restore m_L, col_pattern, and listCol 00329 vals = Map<const VectorSx>(L_save.valuePtr(), nnz); 00330 rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz); 00331 colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1); 00332 col_pattern.fill(-1); 00333 for(Index i=0; i<n; ++i) 00334 listCol[i].clear(); 00335 00336 break; 00337 } 00338 00339 RealScalar rdiag = sqrt(numext::real(diag)); 00340 vals[colPtr[j]] = rdiag; 00341 for (Index k = 0; k<col_nnz; ++k) 00342 { 00343 Index i = col_irow[k]; 00344 //Scale 00345 col_vals(k) /= rdiag; 00346 //Update the remaining diagonals with col_vals 00347 vals[colPtr[i]] -= numext::abs2(col_vals(k)); 00348 } 00349 // Select the largest p elements 00350 // p is the original number of elements in the column (without the diagonal) 00351 Index p = colPtr[j+1] - colPtr[j] - 1 ; 00352 Ref<VectorSx> cvals = col_vals.head(col_nnz); 00353 Ref<VectorIx> cirow = col_irow.head(col_nnz); 00354 internal::QuickSplit(cvals,cirow, p); 00355 // Insert the largest p elements in the matrix 00356 Index cpt = 0; 00357 for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++) 00358 { 00359 vals[i] = col_vals(cpt); 00360 rowIdx[i] = col_irow(cpt); 00361 // restore col_pattern: 00362 col_pattern(col_irow(cpt)) = -1; 00363 cpt++; 00364 } 00365 // Get the first smallest row index and put it after the diagonal element 00366 Index jk = colPtr(j)+1; 00367 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); 00368 } 00369 00370 if(j==n) 00371 { 00372 m_factorizationIsOk = true; 00373 m_info = Success; 00374 } 00375 } while(m_info!=Success); 00376 } 00377 00378 template<typename Scalar, int _UpLo, typename OrderingType> 00379 inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol) 00380 { 00381 if (jk < colPtr(col+1) ) 00382 { 00383 Index p = colPtr(col+1) - jk; 00384 Index minpos; 00385 rowIdx.segment(jk,p).minCoeff(&minpos); 00386 minpos += jk; 00387 if (rowIdx(minpos) != rowIdx(jk)) 00388 { 00389 //Swap 00390 std::swap(rowIdx(jk),rowIdx(minpos)); 00391 std::swap(vals(jk),vals(minpos)); 00392 } 00393 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk); 00394 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col)); 00395 } 00396 } 00397 00398 } // end namespace Eigen 00399 00400 #endif