Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 1 | // This file is part of Eigen, a lightweight C++ template library |
| 2 | // for linear algebra. |
| 3 | // |
| 4 | // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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_LUT_H |
| 11 | #define EIGEN_INCOMPLETE_LUT_H |
| 12 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 13 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 14 | namespace Eigen { |
| 15 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 16 | namespace internal { |
| 17 | |
| 18 | /** \internal |
| 19 | * Compute a quick-sort split of a vector |
| 20 | * On output, the vector row is permuted such that its elements satisfy |
| 21 | * abs(row(i)) >= abs(row(ncut)) if i<ncut |
| 22 | * abs(row(i)) <= abs(row(ncut)) if i>ncut |
| 23 | * \param row The vector of values |
| 24 | * \param ind The array of index for the elements in @p row |
| 25 | * \param ncut The number of largest elements to keep |
| 26 | **/ |
| 27 | template <typename VectorV, typename VectorI, typename Index> |
| 28 | Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) |
| 29 | { |
| 30 | typedef typename VectorV::RealScalar RealScalar; |
| 31 | using std::swap; |
| 32 | using std::abs; |
| 33 | Index mid; |
| 34 | Index n = row.size(); /* length of the vector */ |
| 35 | Index first, last ; |
| 36 | |
| 37 | ncut--; /* to fit the zero-based indices */ |
| 38 | first = 0; |
| 39 | last = n-1; |
| 40 | if (ncut < first || ncut > last ) return 0; |
| 41 | |
| 42 | do { |
| 43 | mid = first; |
| 44 | RealScalar abskey = abs(row(mid)); |
| 45 | for (Index j = first + 1; j <= last; j++) { |
| 46 | if ( abs(row(j)) > abskey) { |
| 47 | ++mid; |
| 48 | swap(row(mid), row(j)); |
| 49 | swap(ind(mid), ind(j)); |
| 50 | } |
| 51 | } |
| 52 | /* Interchange for the pivot element */ |
| 53 | swap(row(mid), row(first)); |
| 54 | swap(ind(mid), ind(first)); |
| 55 | |
| 56 | if (mid > ncut) last = mid - 1; |
| 57 | else if (mid < ncut ) first = mid + 1; |
| 58 | } while (mid != ncut ); |
| 59 | |
| 60 | return 0; /* mid is equal to ncut */ |
| 61 | } |
| 62 | |
| 63 | }// end namespace internal |
| 64 | |
| 65 | /** \ingroup IterativeLinearSolvers_Module |
| 66 | * \class IncompleteLUT |
| 67 | * \brief Incomplete LU factorization with dual-threshold strategy |
| 68 | * |
| 69 | * During the numerical factorization, two dropping rules are used : |
| 70 | * 1) any element whose magnitude is less than some tolerance is dropped. |
| 71 | * This tolerance is obtained by multiplying the input tolerance @p droptol |
| 72 | * by the average magnitude of all the original elements in the current row. |
| 73 | * 2) After the elimination of the row, only the @p fill largest elements in |
| 74 | * the L part and the @p fill largest elements in the U part are kept |
| 75 | * (in addition to the diagonal element ). Note that @p fill is computed from |
| 76 | * the input parameter @p fillfactor which is used the ratio to control the fill_in |
| 77 | * relatively to the initial number of nonzero elements. |
| 78 | * |
| 79 | * The two extreme cases are when @p droptol=0 (to keep all the @p fill*2 largest elements) |
| 80 | * and when @p fill=n/2 with @p droptol being different to zero. |
| 81 | * |
| 82 | * References : Yousef Saad, ILUT: A dual threshold incomplete LU factorization, |
| 83 | * Numerical Linear Algebra with Applications, 1(4), pp 387-402, 1994. |
| 84 | * |
| 85 | * NOTE : The following implementation is derived from the ILUT implementation |
| 86 | * in the SPARSKIT package, Copyright (C) 2005, the Regents of the University of Minnesota |
| 87 | * released under the terms of the GNU LGPL: |
| 88 | * http://www-users.cs.umn.edu/~saad/software/SPARSKIT/README |
| 89 | * However, Yousef Saad gave us permission to relicense his ILUT code to MPL2. |
| 90 | * See the Eigen mailing list archive, thread: ILUT, date: July 8, 2012: |
| 91 | * http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/07/msg00064.html |
| 92 | * alternatively, on GMANE: |
| 93 | * http://comments.gmane.org/gmane.comp.lib.eigen/3302 |
| 94 | */ |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 95 | template <typename _Scalar> |
| 96 | class IncompleteLUT : internal::noncopyable |
| 97 | { |
| 98 | typedef _Scalar Scalar; |
| 99 | typedef typename NumTraits<Scalar>::Real RealScalar; |
| 100 | typedef Matrix<Scalar,Dynamic,1> Vector; |
| 101 | typedef SparseMatrix<Scalar,RowMajor> FactorType; |
| 102 | typedef SparseMatrix<Scalar,ColMajor> PermutType; |
| 103 | typedef typename FactorType::Index Index; |
| 104 | |
| 105 | public: |
| 106 | typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; |
| 107 | |
| 108 | IncompleteLUT() |
| 109 | : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10), |
| 110 | m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false) |
| 111 | {} |
| 112 | |
| 113 | template<typename MatrixType> |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 114 | IncompleteLUT(const MatrixType& mat, const RealScalar& droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 115 | : m_droptol(droptol),m_fillfactor(fillfactor), |
| 116 | m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false) |
| 117 | { |
| 118 | eigen_assert(fillfactor != 0); |
| 119 | compute(mat); |
| 120 | } |
| 121 | |
| 122 | Index rows() const { return m_lu.rows(); } |
| 123 | |
| 124 | Index cols() const { return m_lu.cols(); } |
| 125 | |
| 126 | /** \brief Reports whether previous computation was successful. |
| 127 | * |
| 128 | * \returns \c Success if computation was succesful, |
| 129 | * \c NumericalIssue if the matrix.appears to be negative. |
| 130 | */ |
| 131 | ComputationInfo info() const |
| 132 | { |
| 133 | eigen_assert(m_isInitialized && "IncompleteLUT is not initialized."); |
| 134 | return m_info; |
| 135 | } |
| 136 | |
| 137 | template<typename MatrixType> |
| 138 | void analyzePattern(const MatrixType& amat); |
| 139 | |
| 140 | template<typename MatrixType> |
| 141 | void factorize(const MatrixType& amat); |
| 142 | |
| 143 | /** |
| 144 | * Compute an incomplete LU factorization with dual threshold on the matrix mat |
| 145 | * No pivoting is done in this version |
| 146 | * |
| 147 | **/ |
| 148 | template<typename MatrixType> |
| 149 | IncompleteLUT<Scalar>& compute(const MatrixType& amat) |
| 150 | { |
| 151 | analyzePattern(amat); |
| 152 | factorize(amat); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 153 | m_isInitialized = m_factorizationIsOk; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 154 | return *this; |
| 155 | } |
| 156 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 157 | void setDroptol(const RealScalar& droptol); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 158 | void setFillfactor(int fillfactor); |
| 159 | |
| 160 | template<typename Rhs, typename Dest> |
| 161 | void _solve(const Rhs& b, Dest& x) const |
| 162 | { |
| 163 | x = m_Pinv * b; |
| 164 | x = m_lu.template triangularView<UnitLower>().solve(x); |
| 165 | x = m_lu.template triangularView<Upper>().solve(x); |
| 166 | x = m_P * x; |
| 167 | } |
| 168 | |
| 169 | template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs> |
| 170 | solve(const MatrixBase<Rhs>& b) const |
| 171 | { |
| 172 | eigen_assert(m_isInitialized && "IncompleteLUT is not initialized."); |
| 173 | eigen_assert(cols()==b.rows() |
| 174 | && "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b"); |
| 175 | return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived()); |
| 176 | } |
| 177 | |
| 178 | protected: |
| 179 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 180 | /** keeps off-diagonal entries; drops diagonal entries */ |
| 181 | struct keep_diag { |
| 182 | inline bool operator() (const Index& row, const Index& col, const Scalar&) const |
| 183 | { |
| 184 | return row!=col; |
| 185 | } |
| 186 | }; |
| 187 | |
| 188 | protected: |
| 189 | |
| 190 | FactorType m_lu; |
| 191 | RealScalar m_droptol; |
| 192 | int m_fillfactor; |
| 193 | bool m_analysisIsOk; |
| 194 | bool m_factorizationIsOk; |
| 195 | bool m_isInitialized; |
| 196 | ComputationInfo m_info; |
| 197 | PermutationMatrix<Dynamic,Dynamic,Index> m_P; // Fill-reducing permutation |
| 198 | PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // Inverse permutation |
| 199 | }; |
| 200 | |
| 201 | /** |
| 202 | * Set control parameter droptol |
| 203 | * \param droptol Drop any element whose magnitude is less than this tolerance |
| 204 | **/ |
| 205 | template<typename Scalar> |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 206 | void IncompleteLUT<Scalar>::setDroptol(const RealScalar& droptol) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 207 | { |
| 208 | this->m_droptol = droptol; |
| 209 | } |
| 210 | |
| 211 | /** |
| 212 | * Set control parameter fillfactor |
| 213 | * \param fillfactor This is used to compute the number @p fill_in of largest elements to keep on each row. |
| 214 | **/ |
| 215 | template<typename Scalar> |
| 216 | void IncompleteLUT<Scalar>::setFillfactor(int fillfactor) |
| 217 | { |
| 218 | this->m_fillfactor = fillfactor; |
| 219 | } |
| 220 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 221 | template <typename Scalar> |
| 222 | template<typename _MatrixType> |
| 223 | void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat) |
| 224 | { |
| 225 | // Compute the Fill-reducing permutation |
| 226 | SparseMatrix<Scalar,ColMajor, Index> mat1 = amat; |
| 227 | SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose(); |
| 228 | // Symmetrize the pattern |
| 229 | // FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice. |
| 230 | // on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered... |
| 231 | SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 + mat1; |
| 232 | AtA.prune(keep_diag()); |
| 233 | internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); // Then compute the AMD ordering... |
| 234 | |
| 235 | m_Pinv = m_P.inverse(); // ... and the inverse permutation |
| 236 | |
| 237 | m_analysisIsOk = true; |
| 238 | } |
| 239 | |
| 240 | template <typename Scalar> |
| 241 | template<typename _MatrixType> |
| 242 | void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) |
| 243 | { |
| 244 | using std::sqrt; |
| 245 | using std::swap; |
| 246 | using std::abs; |
| 247 | |
| 248 | eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix"); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 249 | Index n = amat.cols(); // Size of the matrix |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 250 | m_lu.resize(n,n); |
| 251 | // Declare Working vectors and variables |
| 252 | Vector u(n) ; // real values of the row -- maximum size is n -- |
| 253 | VectorXi ju(n); // column position of the values in u -- maximum size is n |
| 254 | VectorXi jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1 |
| 255 | |
| 256 | // Apply the fill-reducing permutation |
| 257 | eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); |
| 258 | SparseMatrix<Scalar,RowMajor, Index> mat; |
| 259 | mat = amat.twistedBy(m_Pinv); |
| 260 | |
| 261 | // Initialization |
| 262 | jr.fill(-1); |
| 263 | ju.fill(0); |
| 264 | u.fill(0); |
| 265 | |
| 266 | // number of largest elements to keep in each row: |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 267 | Index fill_in = static_cast<Index> (amat.nonZeros()*m_fillfactor)/n+1; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 268 | if (fill_in > n) fill_in = n; |
| 269 | |
| 270 | // number of largest nonzero elements to keep in the L and the U part of the current row: |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 271 | Index nnzL = fill_in/2; |
| 272 | Index nnzU = nnzL; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 273 | m_lu.reserve(n * (nnzL + nnzU + 1)); |
| 274 | |
| 275 | // global loop over the rows of the sparse matrix |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 276 | for (Index ii = 0; ii < n; ii++) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 277 | { |
| 278 | // 1 - copy the lower and the upper part of the row i of mat in the working vector u |
| 279 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 280 | Index sizeu = 1; // number of nonzero elements in the upper part of the current row |
| 281 | Index sizel = 0; // number of nonzero elements in the lower part of the current row |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 282 | ju(ii) = ii; |
| 283 | u(ii) = 0; |
| 284 | jr(ii) = ii; |
| 285 | RealScalar rownorm = 0; |
| 286 | |
| 287 | typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii |
| 288 | for (; j_it; ++j_it) |
| 289 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 290 | Index k = j_it.index(); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 291 | if (k < ii) |
| 292 | { |
| 293 | // copy the lower part |
| 294 | ju(sizel) = k; |
| 295 | u(sizel) = j_it.value(); |
| 296 | jr(k) = sizel; |
| 297 | ++sizel; |
| 298 | } |
| 299 | else if (k == ii) |
| 300 | { |
| 301 | u(ii) = j_it.value(); |
| 302 | } |
| 303 | else |
| 304 | { |
| 305 | // copy the upper part |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 306 | Index jpos = ii + sizeu; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 307 | ju(jpos) = k; |
| 308 | u(jpos) = j_it.value(); |
| 309 | jr(k) = jpos; |
| 310 | ++sizeu; |
| 311 | } |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 312 | rownorm += numext::abs2(j_it.value()); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 313 | } |
| 314 | |
| 315 | // 2 - detect possible zero row |
| 316 | if(rownorm==0) |
| 317 | { |
| 318 | m_info = NumericalIssue; |
| 319 | return; |
| 320 | } |
| 321 | // Take the 2-norm of the current row as a relative tolerance |
| 322 | rownorm = sqrt(rownorm); |
| 323 | |
| 324 | // 3 - eliminate the previous nonzero rows |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 325 | Index jj = 0; |
| 326 | Index len = 0; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 327 | while (jj < sizel) |
| 328 | { |
| 329 | // In order to eliminate in the correct order, |
| 330 | // we must select first the smallest column index among ju(jj:sizel) |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 331 | Index k; |
| 332 | Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 333 | k += jj; |
| 334 | if (minrow != ju(jj)) |
| 335 | { |
| 336 | // swap the two locations |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 337 | Index j = ju(jj); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 338 | swap(ju(jj), ju(k)); |
| 339 | jr(minrow) = jj; jr(j) = k; |
| 340 | swap(u(jj), u(k)); |
| 341 | } |
| 342 | // Reset this location |
| 343 | jr(minrow) = -1; |
| 344 | |
| 345 | // Start elimination |
| 346 | typename FactorType::InnerIterator ki_it(m_lu, minrow); |
| 347 | while (ki_it && ki_it.index() < minrow) ++ki_it; |
| 348 | eigen_internal_assert(ki_it && ki_it.col()==minrow); |
| 349 | Scalar fact = u(jj) / ki_it.value(); |
| 350 | |
| 351 | // drop too small elements |
| 352 | if(abs(fact) <= m_droptol) |
| 353 | { |
| 354 | jj++; |
| 355 | continue; |
| 356 | } |
| 357 | |
| 358 | // linear combination of the current row ii and the row minrow |
| 359 | ++ki_it; |
| 360 | for (; ki_it; ++ki_it) |
| 361 | { |
| 362 | Scalar prod = fact * ki_it.value(); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 363 | Index j = ki_it.index(); |
| 364 | Index jpos = jr(j); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 365 | if (jpos == -1) // fill-in element |
| 366 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 367 | Index newpos; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 368 | if (j >= ii) // dealing with the upper part |
| 369 | { |
| 370 | newpos = ii + sizeu; |
| 371 | sizeu++; |
| 372 | eigen_internal_assert(sizeu<=n); |
| 373 | } |
| 374 | else // dealing with the lower part |
| 375 | { |
| 376 | newpos = sizel; |
| 377 | sizel++; |
| 378 | eigen_internal_assert(sizel<=ii); |
| 379 | } |
| 380 | ju(newpos) = j; |
| 381 | u(newpos) = -prod; |
| 382 | jr(j) = newpos; |
| 383 | } |
| 384 | else |
| 385 | u(jpos) -= prod; |
| 386 | } |
| 387 | // store the pivot element |
| 388 | u(len) = fact; |
| 389 | ju(len) = minrow; |
| 390 | ++len; |
| 391 | |
| 392 | jj++; |
| 393 | } // end of the elimination on the row ii |
| 394 | |
| 395 | // reset the upper part of the pointer jr to zero |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 396 | for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 397 | |
| 398 | // 4 - partially sort and insert the elements in the m_lu matrix |
| 399 | |
| 400 | // sort the L-part of the row |
| 401 | sizel = len; |
| 402 | len = (std::min)(sizel, nnzL); |
| 403 | typename Vector::SegmentReturnType ul(u.segment(0, sizel)); |
| 404 | typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel)); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 405 | internal::QuickSplit(ul, jul, len); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 406 | |
| 407 | // store the largest m_fill elements of the L part |
| 408 | m_lu.startVec(ii); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 409 | for(Index k = 0; k < len; k++) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 410 | m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); |
| 411 | |
| 412 | // store the diagonal element |
| 413 | // apply a shifting rule to avoid zero pivots (we are doing an incomplete factorization) |
| 414 | if (u(ii) == Scalar(0)) |
| 415 | u(ii) = sqrt(m_droptol) * rownorm; |
| 416 | m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii); |
| 417 | |
| 418 | // sort the U-part of the row |
| 419 | // apply the dropping rule first |
| 420 | len = 0; |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 421 | for(Index k = 1; k < sizeu; k++) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 422 | { |
| 423 | if(abs(u(ii+k)) > m_droptol * rownorm ) |
| 424 | { |
| 425 | ++len; |
| 426 | u(ii + len) = u(ii + k); |
| 427 | ju(ii + len) = ju(ii + k); |
| 428 | } |
| 429 | } |
| 430 | sizeu = len + 1; // +1 to take into account the diagonal element |
| 431 | len = (std::min)(sizeu, nnzU); |
| 432 | typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1)); |
| 433 | typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1)); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 434 | internal::QuickSplit(uu, juu, len); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 435 | |
| 436 | // store the largest elements of the U part |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 437 | for(Index k = ii + 1; k < ii + len; k++) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 438 | m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); |
| 439 | } |
| 440 | |
| 441 | m_lu.finalize(); |
| 442 | m_lu.makeCompressed(); |
| 443 | |
| 444 | m_factorizationIsOk = true; |
| 445 | m_info = Success; |
| 446 | } |
| 447 | |
| 448 | namespace internal { |
| 449 | |
| 450 | template<typename _MatrixType, typename Rhs> |
| 451 | struct solve_retval<IncompleteLUT<_MatrixType>, Rhs> |
| 452 | : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs> |
| 453 | { |
| 454 | typedef IncompleteLUT<_MatrixType> Dec; |
| 455 | EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) |
| 456 | |
| 457 | template<typename Dest> void evalTo(Dest& dst) const |
| 458 | { |
| 459 | dec()._solve(rhs(),dst); |
| 460 | } |
| 461 | }; |
| 462 | |
| 463 | } // end namespace internal |
| 464 | |
| 465 | } // end namespace Eigen |
| 466 | |
| 467 | #endif // EIGEN_INCOMPLETE_LUT_H |