Improve the IncompleteLLT ... not yet robust

2 files changed, 78 insertions, 50 deletions

diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
index 746d29473..712ec3b6c 100644
--- a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
+++ b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
@@ -35,8 +35,9 @@ class IncompleteCholesky : internal::noncopyable
     typedef typename MatrixType::RealScalar RealScalar; 
     typedef typename MatrixType::Index Index; 
     typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
-    typedef Matrix<Scalar,Dynamic,1> VectorType; 
-    typedef Matrix<Index,Dynamic, 1> IndexType; 
+    typedef Matrix<Scalar,Dynamic,1> ScalarType; 
+    typedef Matrix<Index,Dynamic, 1> IndexType;
+    typedef std::vector<std::list<Index> > VectorList; 
 
   public:
     IncompleteCholesky() {}
@@ -97,6 +98,7 @@ class IncompleteCholesky : internal::noncopyable
     template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs>
     solve(const MatrixBase<Rhs>& b) const
     {
+      eigen_assert(m_factorizationIsOk && "IncompleteLLT did not succeed");
       eigen_assert(m_isInitialized && "IncompleteLLT is not initialized.");
       eigen_assert(cols()==b.rows()
                 && "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b");
@@ -110,6 +112,9 @@ class IncompleteCholesky : internal::noncopyable
     ComputationInfo m_info;
     PermutationType m_perm; 
     
+  private:
+    template <typename IdxType, typename SclType>
+    inline void updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol); 
 }; 
 
 template<typename Scalar, int _UpLo, typename OrderingType>
@@ -129,23 +134,23 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
   else
     m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
   
-  int n = mat.cols(); 
-  
-  Scalar *vals = m_L.valuePtr(); //Values 
-  Index *rowIdx = m_L.innerIndexPtr(); //Row indices 
-  Index *colPtr = m_L.outerIndexPtr(); // Pointer to the beginning of each row
-  VectorType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
-  // Initialize firstElt; 
-  for (int j = 0; j < n-1; j++) firstElt(j) = colPtr[j]+1; 
-  std::vector<std::list<Index> > listCol(n); // listCol(j) is a linked list of columns to update column j
-  VectorType curCol(n); // Store a  nonzero values in each column
-  VectorType irow(n); // Row indices of nonzero elements in each column
+  Index n = m_L.cols(); 
+  Index nnz = m_L.nonZeros();
+  Map<ScalarType> vals(m_L.valuePtr(), nnz); //values
+  Map<IndexType> rowIdx(m_L.innerIndexPtr(), nnz);  //Row indices
+  Map<IndexType> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
+  IndexType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
+  VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
+  ScalarType curCol(n); // Store a  nonzero values in each column
+  IndexType irow(n); // Row indices of nonzero elements in each column
   // jki version of the Cholesky factorization 
-  for (int j=0; j < n; j++)
+  for (int j=0; j < n; ++j)
   {
+    //Debug 
+    bool update_j = false; //This column has received updates
      //Left-looking factorize the column j 
      // First, load the jth column into curCol 
-     Scalar diag = vals[colPtr[j]];  // Lower diagonal matrix with 
+     Scalar diag = vals[colPtr[j]];  // It is assumed that only the lower part is stored
      curCol.setZero();
      irow.setLinSpaced(n,0,n-1); 
      for (int i = colPtr[j] + 1; i < colPtr[j+1]; i++)
@@ -158,48 +163,71 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
      // Browse all previous columns that will update column j
      for(k = listCol[j].begin(); k != listCol[j].end(); k++) 
      {
+       update_j = true;
        int jk = firstElt(*k); // First element to use in the column 
        Scalar a_jk = vals[jk]; 
        diag -= a_jk * a_jk; 
        jk += 1; 
-       for (int i = jk; i < colPtr[*k]; i++)
+       for (int i = jk; i < colPtr[*k+1]; i++)
        {
          curCol(rowIdx[i]) -= vals[i] * a_jk ;
        }
-       firstElt(*k) = jk; 
-       if (jk < colPtr[*k+1]) 
-       {
-         // Add this column to the updating columns list for column *k+1
-         listCol[rowIdx[jk]].push_back(*k); 
-       }
+       updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
      }
      
-     // Select the largest p elements
-     //  p is the original number of elements in the column (without the diagonal)
-     int p = colPtr[j+1] - colPtr[j] - 2 ; 
-     internal::QuickSplit(curCol, irow, p); 
-     if(RealScalar(diag) <= 0) 
-     { //FIXME We can use heuristics (Kershaw, 1978 or above reference ) to get a dynamic shift
-       m_info = NumericalIssue; 
-       return; 
-     }
-     RealScalar rdiag = sqrt(RealScalar(diag));
-     Scalar scal = Scalar(1)/rdiag; 
-     vals[colPtr[j]] = rdiag;
-     // Insert the largest p elements in the matrix and scale them meanwhile  
-     int cpt = 0; 
-     for (int i = colPtr[j]+1; i < colPtr[j+1]; i++)
-     {
-       vals[i] = curCol(cpt) * scal; 
-       rowIdx[i] = irow(cpt); 
-       cpt ++; 
+     if(update_j)
+     {  
+       // Select the largest p elements
+       //  p is the original number of elements in the column (without the diagonal)
+       int p = colPtr[j+1] - colPtr[j] - 1 ; 
+       internal::QuickSplit(curCol, irow, p); 
+       if(RealScalar(diag) <= 0) 
+       { //FIXME We can use heuristics (Kershaw, 1978 or above reference ) to get a dynamic shift
+         std::cerr << "\nNegative diagonal during Incomplete factorization...abort...\n";
+         m_info = NumericalIssue; 
+         return; 
+       }
+       RealScalar rdiag = sqrt(RealScalar(diag));
+       vals[colPtr[j]] = rdiag;
+       Scalar scal = Scalar(1)/rdiag; 
+       // Insert the largest p elements in the matrix and scale them meanwhile  
+       int cpt = 0; 
+       for (int i = colPtr[j]+1; i < colPtr[j+1]; i++)
+       {
+         vals[i] = curCol(cpt) * scal; 
+         rowIdx[i] = irow(cpt); 
+         cpt ++; 
+       }
      }
+     // Get the first smallest row index and put it after the diagonal element
+     Index jk = colPtr(j)+1;
+     updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); 
   }
   m_factorizationIsOk = true; 
   m_isInitialized = true;
   m_info = Success; 
 }
 
+template<typename Scalar, int _UpLo, typename OrderingType>
+template <typename IdxType, typename SclType>
+inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol)
+{
+  if (jk < colPtr(col+1) )
+  {
+    Index p = colPtr(col+1) - jk;
+    Index minpos; 
+    rowIdx.segment(jk,p).minCoeff(&minpos);
+    minpos += jk;
+    if (minpos != rowIdx(jk))
+    {
+      //Swap
+      std::swap(rowIdx(jk),rowIdx(minpos));
+      std::swap(vals(jk),vals(minpos));
+    }
+    firstElt(col) = jk;
+    listCol[rowIdx(jk)].push_back(col);
+  }
+}
 namespace internal {
 
 template<typename _MatrixType, typename Rhs>
diff --git a/unsupported/Eigen/src/IterativeSolvers/Scaling.h b/unsupported/Eigen/src/IterativeSolvers/Scaling.h
index fdef0aca3..ed4ee48c8 100644
--- a/unsupported/Eigen/src/IterativeSolvers/Scaling.h
+++ b/unsupported/Eigen/src/IterativeSolvers/Scaling.h
@@ -7,8 +7,8 @@
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-#ifndef EIGEN_SCALING_H
-#define EIGEN_SCALING_H
+#ifndef EIGEN_ITERSCALING_H
+#define EIGEN_ITERSCALING_H
 /**
   * \ingroup IterativeSolvers_Module
   * \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
@@ -24,7 +24,7 @@
   * VectorXd x(n), b(n);
   * SparseMatrix<double> A;
   * // fill A and b;
-  * Scaling<SparseMatrix<double> > scal; 
+  * IterScaling<SparseMatrix<double> > scal; 
   * // Compute the left and right scaling vectors. The matrix is equilibrated at output
   * scal.computeRef(A); 
   * // Scale the right hand side
@@ -41,10 +41,10 @@
   * 
   * \sa \ref IncompleteLUT 
   */
+namespace Eigen {
 using std::abs; 
-using namespace Eigen;
 template<typename _MatrixType>
-class Scaling
+class IterScaling
 {
   public:
     typedef _MatrixType MatrixType; 
@@ -52,15 +52,15 @@ class Scaling
     typedef typename MatrixType::Index Index;
     
   public:
-    Scaling() { init(); }
+    IterScaling() { init(); }
     
-    Scaling(const MatrixType& matrix)
+    IterScaling(const MatrixType& matrix)
     {
       init();
       compute(matrix);
     }
     
-    ~Scaling() { }
+    ~IterScaling() { }
     
     /** 
      * Compute the left and right diagonal matrices to scale the input matrix @p mat
@@ -181,5 +181,5 @@ class Scaling
     double m_tol; 
     int m_maxits; // Maximum number of iterations allowed
 };
-
+}
 #endif
\ No newline at end of file