Make unsupport sparse solvers use SparseSolverBase

5 files changed, 55 insertions, 26 deletions

diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
index 9fcc8a8d9..fe0bfd948 100644
--- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
@@ -108,6 +108,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
     using Base::m_isInitialized;
     using Base::m_tolerance; 
   public:
+    using Base::_solve_impl;
     typedef _MatrixType MatrixType;
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::Index Index;
@@ -138,6 +139,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
 
   ~DGMRES() {}
   
+#ifndef EIGEN_TEST_EVALUATORS
   /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
     * \a x0 as an initial solution.
     *
@@ -153,10 +155,11 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
     return internal::solve_retval_with_guess
             <DGMRES, Rhs, Guess>(*this, b.derived(), x0);
   }
+#endif
   
   /** \internal */
   template<typename Rhs,typename Dest>
-  void _solveWithGuess(const Rhs& b, Dest& x) const
+  void _solve_with_guess_impl(const Rhs& b, Dest& x) const
   {    
     bool failed = false;
     for(int j=0; j<b.cols(); ++j)
@@ -175,10 +178,10 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
 
   /** \internal */
   template<typename Rhs,typename Dest>
-  void _solve(const Rhs& b, Dest& x) const
+  void _solve_impl(const Rhs& b, MatrixBase<Dest>& x) const
   {
     x = b;
-    _solveWithGuess(b,x);
+    _solve_with_guess_impl(b,x.derived());
   }
   /** 
    * Get the restart value
@@ -522,6 +525,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x,
   return 0; 
 }
 
+#ifndef EIGEN_TEST_EVALUATORS
 namespace internal {
 
   template<typename _MatrixType, typename _Preconditioner, typename Rhs>
@@ -533,10 +537,11 @@ struct solve_retval<DGMRES<_MatrixType, _Preconditioner>, Rhs>
 
   template<typename Dest> void evalTo(Dest& dst) const
   {
-    dec()._solve(rhs(),dst);
+    dec()._solve_impl(rhs(),dst);
   }
 };
 } // end namespace internal
+#endif
 
 } // end namespace Eigen
 #endif 
diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
index 67498705b..fd76a9d2c 100644
--- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
@@ -281,6 +281,7 @@ private:
   int m_restart;
 
 public:
+  using Base::_solve_impl;
   typedef _MatrixType MatrixType;
   typedef typename MatrixType::Scalar Scalar;
   typedef typename MatrixType::Index Index;
@@ -315,6 +316,7 @@ public:
     */
   void set_restart(const int restart) { m_restart=restart; }
 
+#ifndef EIGEN_TEST_EVALUATORS
   /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
     * \a x0 as an initial solution.
     *
@@ -330,10 +332,11 @@ public:
     return internal::solve_retval_with_guess
             <GMRES, Rhs, Guess>(*this, b.derived(), x0);
   }
+#endif
 
   /** \internal */
   template<typename Rhs,typename Dest>
-  void _solveWithGuess(const Rhs& b, Dest& x) const
+  void _solve_with_guess_impl(const Rhs& b, Dest& x) const
   {
     bool failed = false;
     for(int j=0; j<b.cols(); ++j)
@@ -353,18 +356,18 @@ public:
 
   /** \internal */
   template<typename Rhs,typename Dest>
-  void _solve(const Rhs& b, Dest& x) const
+  void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const
   {
     x = b;
     if(x.squaredNorm() == 0) return; // Check Zero right hand side
-    _solveWithGuess(b,x);
+    _solve_with_guess_impl(b,x.derived());
   }
 
 protected:
 
 };
 
-
+#ifndef EIGEN_TEST_EVALUATORS
 namespace internal {
 
   template<typename _MatrixType, typename _Preconditioner, typename Rhs>
@@ -376,12 +379,12 @@ struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs>
 
   template<typename Dest> void evalTo(Dest& dst) const
   {
-    dec()._solve(rhs(),dst);
+    dec()._solve_impl(rhs(),dst);
   }
 };
 
 } // end namespace internal
-
+#endif
 } // end namespace Eigen
 
 #endif // EIGEN_GMRES_H
diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
index 661c1f2e0..1ee1c89b2 100644
--- a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
+++ b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
@@ -27,8 +27,11 @@ namespace Eigen {
  */
 
 template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> >
-class IncompleteCholesky : internal::noncopyable
+class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
 {
+  protected:
+    typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
+    using Base::m_isInitialized;
   public:
     typedef SparseMatrix<Scalar,ColMajor> MatrixType;
     typedef _OrderingType OrderingType;
@@ -89,7 +92,7 @@ class IncompleteCholesky : internal::noncopyable
     }
     
     template<typename Rhs, typename Dest>
-    void _solve(const Rhs& b, Dest& x) const
+    void _solve_impl(const Rhs& b, Dest& x) const
     {
       eigen_assert(m_factorizationIsOk && "factorize() should be called first");
       if (m_perm.rows() == b.rows())
@@ -103,6 +106,8 @@ class IncompleteCholesky : internal::noncopyable
         x = m_perm * x;
       x = m_scal.asDiagonal() * x;
     }
+    
+#ifndef EIGEN_TEST_EVALUATORS
     template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs>
     solve(const MatrixBase<Rhs>& b) const
     {
@@ -112,13 +117,14 @@ class IncompleteCholesky : internal::noncopyable
                 && "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b");
       return internal::solve_retval<IncompleteCholesky, Rhs>(*this, b.derived());
     }
+#endif
+
   protected:
     SparseMatrix<Scalar,ColMajor> m_L;  // The lower part stored in CSC
     ScalarType m_scal; // The vector for scaling the matrix 
     Scalar m_shift; //The initial shift parameter
     bool m_analysisIsOk; 
     bool m_factorizationIsOk; 
-    bool m_isInitialized;
     ComputationInfo m_info;
     PermutationType m_perm; 
     
@@ -256,6 +262,8 @@ inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const Idx
     listCol[rowIdx(jk)].push_back(col);
   }
 }
+
+#ifndef EIGEN_TEST_EVALUATORS
 namespace internal {
 
 template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs>
@@ -272,6 +280,7 @@ struct solve_retval<IncompleteCholesky<_Scalar,  _UpLo, OrderingType>, Rhs>
 };
 
 } // end namespace internal
+#endif
 
 } // end namespace Eigen 
 
diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h
index 67e780181..e86f65644 100644
--- a/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h
+++ b/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h
@@ -13,8 +13,12 @@
 namespace Eigen { 
 
 template <typename _Scalar>
-class IncompleteLU
+class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> >
 {
+  protected:
+    typedef SparseSolverBase<IncompleteLU<_Scalar> > Base;
+    using Base::m_isInitialized;
+    
     typedef _Scalar Scalar;
     typedef Matrix<Scalar,Dynamic,1> Vector;
     typedef typename Vector::Index Index;
@@ -23,10 +27,10 @@ class IncompleteLU
   public:
     typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
 
-    IncompleteLU() : m_isInitialized(false) {}
+    IncompleteLU() {}
 
     template<typename MatrixType>
-    IncompleteLU(const MatrixType& mat) : m_isInitialized(false)
+    IncompleteLU(const MatrixType& mat)
     {
       compute(mat);
     }
@@ -71,12 +75,13 @@ class IncompleteLU
     }
 
     template<typename Rhs, typename Dest>
-    void _solve(const Rhs& b, Dest& x) const
+    void _solve_impl(const Rhs& b, Dest& x) const
     {
       x = m_lu.template triangularView<UnitLower>().solve(b);
       x = m_lu.template triangularView<Upper>().solve(x);
     }
 
+#ifndef EIGEN_TEST_EVALUATORS
     template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs>
     solve(const MatrixBase<Rhs>& b) const
     {
@@ -85,12 +90,13 @@ class IncompleteLU
                 && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b");
       return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived());
     }
+#endif
 
   protected:
     FactorType m_lu;
-    bool m_isInitialized;
 };
 
+#ifndef EIGEN_TEST_EVALUATORS
 namespace internal {
 
 template<typename _MatrixType, typename Rhs>
@@ -102,11 +108,12 @@ struct solve_retval<IncompleteLU<_MatrixType>, Rhs>
 
   template<typename Dest> void evalTo(Dest& dst) const
   {
-    dec()._solve(rhs(),dst);
+    dec()._solve_impl(rhs(),dst);
   }
 };
 
 } // end namespace internal
+#endif
 
 } // end namespace Eigen
 
diff --git a/unsupported/Eigen/src/IterativeSolvers/MINRES.h b/unsupported/Eigen/src/IterativeSolvers/MINRES.h
index 98f9ecc17..28d5c692d 100644
--- a/unsupported/Eigen/src/IterativeSolvers/MINRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/MINRES.h
@@ -2,7 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2012 Giacomo Po <gpo@ucla.edu>
-// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -217,6 +217,7 @@ namespace Eigen {
         using Base::m_info;
         using Base::m_isInitialized;
     public:
+        using Base::_solve_impl;
         typedef _MatrixType MatrixType;
         typedef typename MatrixType::Scalar Scalar;
         typedef typename MatrixType::Index Index;
@@ -244,7 +245,8 @@ namespace Eigen {
         
         /** Destructor. */
         ~MINRES(){}
-		
+
+#ifndef EIGEN_TEST_EVALUATORS
         /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
          * \a x0 as an initial solution.
          *
@@ -260,10 +262,11 @@ namespace Eigen {
             return internal::solve_retval_with_guess
             <MINRES, Rhs, Guess>(*this, b.derived(), x0);
         }
+#endif
         
         /** \internal */
         template<typename Rhs,typename Dest>
-        void _solveWithGuess(const Rhs& b, Dest& x) const
+        void _solve_with_guess_impl(const Rhs& b, Dest& x) const
         {
             m_iterations = Base::maxIterations();
             m_error = Base::m_tolerance;
@@ -284,16 +287,17 @@ namespace Eigen {
         
         /** \internal */
         template<typename Rhs,typename Dest>
-        void _solve(const Rhs& b, Dest& x) const
+        void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const
         {
             x.setZero();
-            _solveWithGuess(b,x);
+            _solve_with_guess_impl(b,x.derived());
         }
         
     protected:
         
     };
     
+#ifndef EIGEN_TEST_EVALUATORS
     namespace internal {
         
         template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
@@ -305,12 +309,13 @@ namespace Eigen {
             
             template<typename Dest> void evalTo(Dest& dst) const
             {
-                dec()._solve(rhs(),dst);
+                dec()._solve_impl(rhs(),dst);
             }
         };
         
     } // end namespace internal
-    
+#endif
+
 } // end namespace Eigen
 
 #endif // EIGEN_MINRES_H