00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_INCOMPLETE_LUT_H
00026 #define EIGEN_INCOMPLETE_LUT_H
00027
00028 namespace Eigen {
00029
00053 template <typename _Scalar>
00054 class IncompleteLUT : internal::noncopyable
00055 {
00056 typedef _Scalar Scalar;
00057 typedef typename NumTraits<Scalar>::Real RealScalar;
00058 typedef Matrix<Scalar,Dynamic,1> Vector;
00059 typedef SparseMatrix<Scalar,RowMajor> FactorType;
00060 typedef SparseMatrix<Scalar,ColMajor> PermutType;
00061 typedef typename FactorType::Index Index;
00062
00063 public:
00064 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
00065
00066 IncompleteLUT()
00067 : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10),
00068 m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false)
00069 {}
00070
00071 template<typename MatrixType>
00072 IncompleteLUT(const MatrixType& mat, RealScalar droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10)
00073 : m_droptol(droptol),m_fillfactor(fillfactor),
00074 m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false)
00075 {
00076 eigen_assert(fillfactor != 0);
00077 compute(mat);
00078 }
00079
00080 Index rows() const { return m_lu.rows(); }
00081
00082 Index cols() const { return m_lu.cols(); }
00083
00089 ComputationInfo info() const
00090 {
00091 eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
00092 return m_info;
00093 }
00094
00095 template<typename MatrixType>
00096 void analyzePattern(const MatrixType& amat);
00097
00098 template<typename MatrixType>
00099 void factorize(const MatrixType& amat);
00100
00106 template<typename MatrixType>
00107 IncompleteLUT<Scalar>& compute(const MatrixType& amat)
00108 {
00109 analyzePattern(amat);
00110 factorize(amat);
00111 eigen_assert(m_factorizationIsOk == true);
00112 m_isInitialized = true;
00113 return *this;
00114 }
00115
00116 void setDroptol(RealScalar droptol);
00117 void setFillfactor(int fillfactor);
00118
00119 template<typename Rhs, typename Dest>
00120 void _solve(const Rhs& b, Dest& x) const
00121 {
00122 x = m_Pinv * b;
00123 x = m_lu.template triangularView<UnitLower>().solve(x);
00124 x = m_lu.template triangularView<Upper>().solve(x);
00125 x = m_P * x;
00126 }
00127
00128 template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs>
00129 solve(const MatrixBase<Rhs>& b) const
00130 {
00131 eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
00132 eigen_assert(cols()==b.rows()
00133 && "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
00134 return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived());
00135 }
00136
00137 protected:
00138
00139 template <typename VectorV, typename VectorI>
00140 int QuickSplit(VectorV &row, VectorI &ind, int ncut);
00141
00142
00144 struct keep_diag {
00145 inline bool operator() (const Index& row, const Index& col, const Scalar&) const
00146 {
00147 return row!=col;
00148 }
00149 };
00150
00151 protected:
00152
00153 FactorType m_lu;
00154 RealScalar m_droptol;
00155 int m_fillfactor;
00156 bool m_analysisIsOk;
00157 bool m_factorizationIsOk;
00158 bool m_isInitialized;
00159 ComputationInfo m_info;
00160 PermutationMatrix<Dynamic,Dynamic,Index> m_P;
00161 PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;
00162 };
00163
00168 template<typename Scalar>
00169 void IncompleteLUT<Scalar>::setDroptol(RealScalar droptol)
00170 {
00171 this->m_droptol = droptol;
00172 }
00173
00178 template<typename Scalar>
00179 void IncompleteLUT<Scalar>::setFillfactor(int fillfactor)
00180 {
00181 this->m_fillfactor = fillfactor;
00182 }
00183
00184
00194 template <typename Scalar>
00195 template <typename VectorV, typename VectorI>
00196 int IncompleteLUT<Scalar>::QuickSplit(VectorV &row, VectorI &ind, int ncut)
00197 {
00198 using std::swap;
00199 int mid;
00200 int n = row.size();
00201 int first, last ;
00202
00203 ncut--;
00204 first = 0;
00205 last = n-1;
00206 if (ncut < first || ncut > last ) return 0;
00207
00208 do {
00209 mid = first;
00210 RealScalar abskey = std::abs(row(mid));
00211 for (int j = first + 1; j <= last; j++) {
00212 if ( std::abs(row(j)) > abskey) {
00213 ++mid;
00214 swap(row(mid), row(j));
00215 swap(ind(mid), ind(j));
00216 }
00217 }
00218
00219 swap(row(mid), row(first));
00220 swap(ind(mid), ind(first));
00221
00222 if (mid > ncut) last = mid - 1;
00223 else if (mid < ncut ) first = mid + 1;
00224 } while (mid != ncut );
00225
00226 return 0;
00227 }
00228
00229 template <typename Scalar>
00230 template<typename _MatrixType>
00231 void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat)
00232 {
00233
00234 SparseMatrix<Scalar,ColMajor, Index> mat1 = amat;
00235 SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose();
00236
00237
00238
00239 SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 + mat1;
00240 AtA.prune(keep_diag());
00241 internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P);
00242
00243 m_Pinv = m_P.inverse();
00244
00245 m_analysisIsOk = true;
00246 }
00247
00248 template <typename Scalar>
00249 template<typename _MatrixType>
00250 void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
00251 {
00252 using std::sqrt;
00253 using std::swap;
00254 using std::abs;
00255
00256 eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
00257 int n = amat.cols();
00258 m_lu.resize(n,n);
00259
00260 Vector u(n) ;
00261 VectorXi ju(n);
00262 VectorXi jr(n);
00263
00264
00265 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00266 SparseMatrix<Scalar,RowMajor, Index> mat;
00267 mat = amat.twistedBy(m_Pinv);
00268
00269
00270 jr.fill(-1);
00271 ju.fill(0);
00272 u.fill(0);
00273
00274
00275 int fill_in = static_cast<int> (amat.nonZeros()*m_fillfactor)/n+1;
00276 if (fill_in > n) fill_in = n;
00277
00278
00279 int nnzL = fill_in/2;
00280 int nnzU = nnzL;
00281 m_lu.reserve(n * (nnzL + nnzU + 1));
00282
00283
00284 for (int ii = 0; ii < n; ii++)
00285 {
00286
00287
00288 int sizeu = 1;
00289 int sizel = 0;
00290 ju(ii) = ii;
00291 u(ii) = 0;
00292 jr(ii) = ii;
00293 RealScalar rownorm = 0;
00294
00295 typename FactorType::InnerIterator j_it(mat, ii);
00296 for (; j_it; ++j_it)
00297 {
00298 int k = j_it.index();
00299 if (k < ii)
00300 {
00301
00302 ju(sizel) = k;
00303 u(sizel) = j_it.value();
00304 jr(k) = sizel;
00305 ++sizel;
00306 }
00307 else if (k == ii)
00308 {
00309 u(ii) = j_it.value();
00310 }
00311 else
00312 {
00313
00314 int jpos = ii + sizeu;
00315 ju(jpos) = k;
00316 u(jpos) = j_it.value();
00317 jr(k) = jpos;
00318 ++sizeu;
00319 }
00320 rownorm += internal::abs2(j_it.value());
00321 }
00322
00323
00324 if(rownorm==0)
00325 {
00326 m_info = NumericalIssue;
00327 return;
00328 }
00329
00330 rownorm = sqrt(rownorm);
00331
00332
00333 int jj = 0;
00334 int len = 0;
00335 while (jj < sizel)
00336 {
00337
00338
00339 int k;
00340 int minrow = ju.segment(jj,sizel-jj).minCoeff(&k);
00341 k += jj;
00342 if (minrow != ju(jj))
00343 {
00344
00345 int j = ju(jj);
00346 swap(ju(jj), ju(k));
00347 jr(minrow) = jj; jr(j) = k;
00348 swap(u(jj), u(k));
00349 }
00350
00351 jr(minrow) = -1;
00352
00353
00354 typename FactorType::InnerIterator ki_it(m_lu, minrow);
00355 while (ki_it && ki_it.index() < minrow) ++ki_it;
00356 eigen_internal_assert(ki_it && ki_it.col()==minrow);
00357 Scalar fact = u(jj) / ki_it.value();
00358
00359
00360 if(abs(fact) <= m_droptol)
00361 {
00362 jj++;
00363 continue;
00364 }
00365
00366
00367 ++ki_it;
00368 for (; ki_it; ++ki_it)
00369 {
00370 Scalar prod = fact * ki_it.value();
00371 int j = ki_it.index();
00372 int jpos = jr(j);
00373 if (jpos == -1)
00374 {
00375 int newpos;
00376 if (j >= ii)
00377 {
00378 newpos = ii + sizeu;
00379 sizeu++;
00380 eigen_internal_assert(sizeu<=n);
00381 }
00382 else
00383 {
00384 newpos = sizel;
00385 sizel++;
00386 eigen_internal_assert(sizel<=ii);
00387 }
00388 ju(newpos) = j;
00389 u(newpos) = -prod;
00390 jr(j) = newpos;
00391 }
00392 else
00393 u(jpos) -= prod;
00394 }
00395
00396 u(len) = fact;
00397 ju(len) = minrow;
00398 ++len;
00399
00400 jj++;
00401 }
00402
00403
00404 for(int k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
00405
00406
00407
00408
00409 sizel = len;
00410 len = (std::min)(sizel, nnzL);
00411 typename Vector::SegmentReturnType ul(u.segment(0, sizel));
00412 typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel));
00413 QuickSplit(ul, jul, len);
00414
00415
00416 m_lu.startVec(ii);
00417 for(int k = 0; k < len; k++)
00418 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
00419
00420
00421
00422 if (u(ii) == Scalar(0))
00423 u(ii) = sqrt(m_droptol) * rownorm;
00424 m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
00425
00426
00427
00428 len = 0;
00429 for(int k = 1; k < sizeu; k++)
00430 {
00431 if(abs(u(ii+k)) > m_droptol * rownorm )
00432 {
00433 ++len;
00434 u(ii + len) = u(ii + k);
00435 ju(ii + len) = ju(ii + k);
00436 }
00437 }
00438 sizeu = len + 1;
00439 len = (std::min)(sizeu, nnzU);
00440 typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
00441 typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
00442 QuickSplit(uu, juu, len);
00443
00444
00445 for(int k = ii + 1; k < ii + len; k++)
00446 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
00447 }
00448
00449 m_lu.finalize();
00450 m_lu.makeCompressed();
00451
00452 m_factorizationIsOk = true;
00453 m_info = Success;
00454 }
00455
00456 namespace internal {
00457
00458 template<typename _MatrixType, typename Rhs>
00459 struct solve_retval<IncompleteLUT<_MatrixType>, Rhs>
00460 : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs>
00461 {
00462 typedef IncompleteLUT<_MatrixType> Dec;
00463 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00464
00465 template<typename Dest> void evalTo(Dest& dst) const
00466 {
00467 dec()._solve(rhs(),dst);
00468 }
00469 };
00470
00471 }
00472
00473 }
00474
00475 #endif // EIGEN_INCOMPLETE_LUT_H
00476