Incomplete Cholesky preconditioner... not yet stable

1 files changed, 50 insertions, 51 deletions

diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
index 224304f0e..5a71531cd 100644
--- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
+++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
@@ -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++)