devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Incomplete Cholesky preconditioner... not yet stable

Browse Source
This commit is contained in:
Desire NUENTSA
2012-09-11 12:12:19 +02:00
parent 504edbddb1
commit 45672e724e
4 changed files with 275 additions and 53 deletions
Download Patch File Download Diff File Expand all files Collapse all files

101
Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
View File

@@ -10,8 +10,56 @@
#ifndef EIGEN_INCOMPLETE_LUT_H
#define EIGEN_INCOMPLETE_LUT_H
namespace Eigen {
namespace internal {
/**
* Compute a quick-sort split of a vector
* On output, the vector row is permuted such that its elements satisfy
* abs(row(i)) >= abs(row(ncut)) if i<ncut
* abs(row(i)) <= abs(row(ncut)) if i>ncut
* \param row The vector of values
* \param ind The array of index for the elements in @p row
* \param ncut The number of largest elements to keep
**/
template <typename VectorV, typename VectorI>
int QuickSplit(VectorV &row, VectorI &ind, int ncut)
{
typedef typename VectorV::RealScalar RealScalar;
using std::swap;
int mid;
int n = row.size(); /* length of the vector */
int first, last ;
ncut--; /* to fit the zero-based indices */
first = 0;
last = n-1;
if (ncut < first || ncut > last ) return 0;
do {
mid = first;
RealScalar abskey = std::abs(row(mid));
for (int j = first + 1; j <= last; j++) {
if ( std::abs(row(j)) > abskey) {
++mid;
swap(row(mid), row(j));
swap(ind(mid), ind(j));
}
}
/* Interchange for the pivot element */
swap(row(mid), row(first));
swap(ind(mid), ind(first));
if (mid > ncut) last = mid - 1;
else if (mid < ncut ) first = mid + 1;
} while (mid != ncut );
return 0; /* mid is equal to ncut */
}
}// end namespace internal
/**
* \brief Incomplete LU factorization with dual-threshold strategy
* During the numerical factorization, two dropping rules are used :
@@ -126,10 +174,6 @@ class IncompleteLUT : internal::noncopyable
protected:
template <typename VectorV, typename VectorI>
int QuickSplit(VectorV &row, VectorI &ind, int ncut);
/** keeps off-diagonal entries; drops diagonal entries */
struct keep_diag {
inline bool operator() (const Index& row, const Index& col, const Scalar&) const
@@ -171,51 +215,6 @@ void IncompleteLUT<Scalar>::setFillfactor(int fillfactor)
this->m_fillfactor = fillfactor;
}
/**
* Compute a quick-sort split of a vector
* On output, the vector row is permuted such that its elements satisfy
* abs(row(i)) >= abs(row(ncut)) if i<ncut
* abs(row(i)) <= abs(row(ncut)) if i>ncut
* \param row The vector of values
* \param ind The array of index for the elements in @p row
* \param ncut The number of largest elements to keep
**/
template <typename Scalar>
template <typename VectorV, typename VectorI>
int IncompleteLUT<Scalar>::QuickSplit(VectorV &row, VectorI &ind, int ncut)
{
using std::swap;
int mid;
int n = row.size(); /* length of the vector */
int first, last ;
ncut--; /* to fit the zero-based indices */
first = 0;
last = n-1;
if (ncut < first || ncut > last ) return 0;
do {
mid = first;
RealScalar abskey = std::abs(row(mid));
for (int j = first + 1; j <= last; j++) {
if ( std::abs(row(j)) > abskey) {
++mid;
swap(row(mid), row(j));
swap(ind(mid), ind(j));
}
}
/* Interchange for the pivot element */
swap(row(mid), row(first));
swap(ind(mid), ind(first));
if (mid > ncut) last = mid - 1;
else if (mid < ncut ) first = mid + 1;
} while (mid != ncut );
return 0; /* mid is equal to ncut */
}
template <typename Scalar>
template<typename _MatrixType>
void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat)
@@ -400,7 +399,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
len = (std::min)(sizel, nnzL);
typename Vector::SegmentReturnType ul(u.segment(0, sizel));
typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel));
QuickSplit(ul, jul, len);
internal::QuickSplit(ul, jul, len);
// store the largest m_fill elements of the L part
m_lu.startVec(ii);
@@ -429,7 +428,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
len = (std::min)(sizeu, nnzU);
typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
QuickSplit(uu, juu, len);
internal::QuickSplit(uu, juu, len);
// store the largest elements of the U part
for(int k = ii + 1; k < ii + len; k++)