Split the computation of the ILUT into two steps

Eigen/src/IterativeLinearSolvers/BiCGSTAB.h

@@ -35,9 +35,10 @@ namespace internal {
   *                approximation of Ax=b (regardless of b)
   * \param iters On input the max number of iteration, on output the number of performed iterations.
   * \param tol_error On input the tolerance error, on output an estimation of the relative error.
+  * \return false in the case of numerical issue, for example a break down of BiCGSTAB. 
   */
 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
-void bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
+bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
               const Preconditioner& precond, int& iters,
               typename Dest::RealScalar& tol_error)
 {
@@ -46,7 +47,6 @@ void bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
   typedef typename Dest::RealScalar RealScalar;
   typedef typename Dest::Scalar Scalar;
   typedef Matrix<Scalar,Dynamic,1> VectorType;
-  
   RealScalar tol = tol_error;
   int maxIters = iters;
 
@@ -70,10 +70,11 @@ void bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
 
   while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters )
   {
+//     std::cout<<i<<" : Relative residual norm " << sqrt(r.squaredNorm()/r0_sqnorm)<<"\n"; 
     Scalar rho_old = rho;
 
     rho = r0.dot(r);
-    eigen_assert((rho != Scalar(0)) && "BiCGSTAB BROKE DOWN !!!");
+    if (rho == Scalar(0)) return false; /* New search directions cannot be found */
     Scalar beta = (rho/rho_old) * (alpha / w);
     p = r + beta * (p - w * v);
     
@@ -94,6 +95,7 @@ void bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
   }
   tol_error = sqrt(r.squaredNorm()/r0_sqnorm);
   iters = i;
+  return true; 
 }
 
 }
@@ -214,17 +216,18 @@ public:
   template<typename Rhs,typename Dest>
   void _solveWithGuess(const Rhs& b, Dest& x) const
   {    
+    bool ok; 
     for(int j=0; j<b.cols(); ++j)
     {
       m_iterations = Base::maxIterations();
       m_error = Base::m_tolerance;
       
       typename Dest::ColXpr xj(x,j);
-      internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
+      ok = internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
     }
-    
+    if (ok == false) m_info = NumericalIssue; 
+    else   m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
     m_isInitialized = true;
-    m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
   }
 
   /** \internal */

Eigen/src/IterativeLinearSolvers/IncompleteLUT.h

@@ -63,12 +63,13 @@ class IncompleteLUT
   public:
     typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
     
-    IncompleteLUT() : m_droptol(NumTraits<Scalar>::dummy_precision()),m_fillfactor(50),m_isInitialized(false) {}; 
+    IncompleteLUT() : m_droptol(NumTraits<Scalar>::dummy_precision()),m_fillfactor(10),m_isInitialized(false),m_analysisIsOk(false),m_factorizationIsOk(false) {}; 
     
     template<typename MatrixType>
     IncompleteLUT(const MatrixType& mat, RealScalar droptol, int fillfactor) 
-    : m_droptol(droptol),m_fillfactor(fillfactor),m_isInitialized(false)
+    : m_droptol(droptol),m_fillfactor(fillfactor),m_isInitialized(false),m_analysisIsOk(false),m_factorizationIsOk(false)
     {
+      eigen_assert(fillfactor != 0);
       compute(mat); 
     }
     
@@ -76,15 +77,24 @@ class IncompleteLUT
     
     Index cols() const { return m_lu.cols(); }
     
+    template<typename MatrixType>
+    void analyzePattern(const MatrixType& amat)
+    {
+      /* Compute the Fill-reducing permutation */
+      SparseMatrix<Scalar,ColMajor, Index> mat1 = amat;
+      SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose();
+      SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 * mat1; /* Symmetrize the pattern */
+      AtA.prune(keep_diag());
+      internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P);  /* Then compute the AMD ordering... */
       
-/**
+      m_Pinv  = m_P.inverse(); /* ... and the inverse permutation */
- * Compute an incomplete LU factorization with dual threshold on the matrix mat
+      m_analysisIsOk = true; 
- * No pivoting is done in this version
+    }
- * 
+    
- **/
+    template<typename MatrixType>
-template<typename MatrixType>
+    void  factorize(const MatrixType& amat)
-IncompleteLUT<Scalar>& compute(const MatrixType& amat)
+    {
-{
+      eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
       int n = amat.cols();  /* Size of the matrix */
       m_lu.resize(n,n); 
       int fill_in; /* Number of largest elements to keep in each row */
@@ -99,20 +109,12 @@ IncompleteLUT<Scalar>& compute(const MatrixType& amat)
       Scalar fact, prod;
       RealScalar rownorm;
 
-  /* Compute the Fill-reducing permutation */
+      /* Apply the fill-reducing permutation */
-  SparseMatrix<Scalar,ColMajor, Index> mat1 = amat;
-  SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose();
-  SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 * mat1; /* Symmetrize the pattern */
-  AtA.prune(keep_diag());
-  internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P);  /* Then compute the AMD ordering... */
-  
-  m_Pinv  = m_P.inverse(); /* ... and the inverse permutation */
-  
-//   m_Pinv.indices().setLinSpaced(0,n);
-//   m_P.indices().setLinSpaced(0,n);
       
+      eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
       SparseMatrix<Scalar,RowMajor, Index> mat;
       mat = amat.twistedBy(m_Pinv);  
+      
       /* Initialization */
       fact = 0;
       jr.fill(-1); 
@@ -270,13 +272,28 @@ IncompleteLUT<Scalar>& compute(const MatrixType& amat)
       } /* End global for-loop */
       m_lu.finalize();
       m_lu.makeCompressed(); /* NOTE To save the extra space */
+      
+      m_factorizationIsOk = true; 
+    }
+    
+    /**
+    * Compute an incomplete LU factorization with dual threshold on the matrix mat
+    * No pivoting is done in this version
+    * 
+    **/
+    template<typename MatrixType>
+    IncompleteLUT<Scalar>& compute(const MatrixType& amat)
+    {
+      analyzePattern(amat); 
+      factorize(amat);
+      eigen_assert(m_factorizationIsOk == true); 
       m_isInitialized = true;
       return *this;
-}
+    }
 
     
     void setDroptol(RealScalar droptol); 
-    void setFill(int fillfactor); 
+    void setFillfactor(int fillfactor); 
     
     
     
@@ -303,6 +320,8 @@ protected:
     FactorType m_lu;
     RealScalar m_droptol;
     int m_fillfactor;   
+    bool m_factorizationIsOk; 
+    bool m_analysisIsOk; 
     bool m_isInitialized;
     template <typename VectorV, typename VectorI> 
     int QuickSplit(VectorV &row, VectorI &ind, int ncut); 
@@ -333,7 +352,7 @@ void IncompleteLUT<Scalar>::setDroptol(RealScalar droptol)
  * \param fillfactor  This is used to compute the  number @p fill_in of largest elements to keep on each row. 
  **/ 
 template<typename Scalar>
-void IncompleteLUT<Scalar>::setFill(int fillfactor)
+void IncompleteLUT<Scalar>::setFillfactor(int fillfactor)
 {
   this->m_fillfactor = fillfactor;   
 }