From eb168ef8ed29ac71eae6c10480a0239b5783730c Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Mon, 27 Feb 2012 14:10:26 +0100
Subject: [PATCH] add analyzePattern/factorize API to iterative solvers and
 basic preconditioners

---
 .../BasicPreconditioners.h                    | 22 ++++++++++-
 .../IterativeLinearSolvers/IncompleteLUT.h    |  8 ++--
 .../IterativeSolverBase.h                     | 39 ++++++++++++++++++-
 3 files changed, 62 insertions(+), 7 deletions(-)

diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
index b1da2f8ce..b8b7641a5 100644
--- a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
+++ b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
@@ -62,9 +62,15 @@ class DiagonalPreconditioner
 
     Index rows() const { return m_invdiag.size(); }
     Index cols() const { return m_invdiag.size(); }
-
+    
     template<typename MatrixType>
-    DiagonalPreconditioner& compute(const MatrixType& mat)
+    DiagonalPreconditioner& analyzePattern(const MatrixType& )
+    {
+      return *this;
+    }
+    
+    template<typename MatrixType>
+    DiagonalPreconditioner& factorize(const MatrixType& mat)
     {
       m_invdiag.resize(mat.cols());
       for(int j=0; j<mat.outerSize(); ++j)
@@ -79,6 +85,12 @@ class DiagonalPreconditioner
       m_isInitialized = true;
       return *this;
     }
+    
+    template<typename MatrixType>
+    DiagonalPreconditioner& compute(const MatrixType& mat)
+    {
+      return factorize(mat);
+    }
 
     template<typename Rhs, typename Dest>
     void _solve(const Rhs& b, Dest& x) const
@@ -130,6 +142,12 @@ class IdentityPreconditioner
 
     template<typename MatrixType>
     IdentityPreconditioner(const MatrixType& ) {}
+    
+    template<typename MatrixType>
+    IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
+    
+    template<typename MatrixType>
+    IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
 
     template<typename MatrixType>
     IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
index 3f136cc4a..fd2a188e5 100644
--- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
+++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
@@ -97,10 +97,10 @@ class IncompleteLUT
     void factorize(const MatrixType& amat);
     
     /**
-    * Compute an incomplete LU factorization with dual threshold on the matrix mat
-    * No pivoting is done in this version
-    * 
-    **/
+      * 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)
     {
diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h
index bcf1669e2..70d86f60b 100644
--- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h
+++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h
@@ -70,6 +70,39 @@ public:
   }
 
   ~IterativeSolverBase() {}
+  
+  /** Initializes the iterative solver for the sparcity pattern of the matrix \a A for further solving \c Ax=b problems.
+    *
+    * Currently, this function mostly call analyzePattern on the preconditioner. In the future
+    * we might, for instance, implement column reodering for faster matrix vector products.
+    */
+  Derived& analyzePattern(const MatrixType& A)
+  {
+    m_preconditioner.analyzePattern(A);
+    m_isInitialized = true;
+    m_analysisIsOk = true;
+    m_info = Success;
+    return derived();
+  }
+  
+  /** Initializes the iterative solver with the numerical values of the matrix \a A for further solving \c Ax=b problems.
+    *
+    * Currently, this function mostly call factorize on the preconditioner.
+    *
+    * \warning this class stores a reference to the matrix A as well as some
+    * precomputed values that depend on it. Therefore, if \a A is changed
+    * this class becomes invalid. Call compute() to update it with the new
+    * matrix A, or modify a copy of A.
+    */
+  Derived& factorize(const MatrixType& A)
+  {
+    eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 
+    mp_matrix = &A;
+    m_preconditioner.factorize(A);
+    m_factorizationIsOk = true;
+    m_info = Success;
+    return derived();
+  }
 
   /** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems.
     *
@@ -86,6 +119,8 @@ public:
     mp_matrix = &A;
     m_preconditioner.compute(A);
     m_isInitialized = true;
+    m_analysisIsOk = true;
+    m_factorizationIsOk = true;
     m_info = Success;
     return derived();
   }
@@ -194,6 +229,8 @@ protected:
   void init()
   {
     m_isInitialized = false;
+    m_analysisIsOk = false;
+    m_factorizationIsOk = false;
     m_maxIterations = -1;
     m_tolerance = NumTraits<Scalar>::epsilon();
   }
@@ -206,7 +243,7 @@ protected:
   mutable RealScalar m_error;
   mutable int m_iterations;
   mutable ComputationInfo m_info;
-  mutable bool m_isInitialized;
+  mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk;
 };
 
 namespace internal {