update

1 files changed, 31 insertions, 53 deletions

diff --git a/eigen/unsupported/Eigen/src/IterativeSolvers/MINRES.h b/eigen/unsupported/Eigen/src/IterativeSolvers/MINRES.h
index 670f274..256990c 100644
--- a/eigen/unsupported/Eigen/src/IterativeSolvers/MINRES.h
+++ b/eigen/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
@@ -29,7 +29,7 @@ namespace Eigen {
         template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
         EIGEN_DONT_INLINE
         void minres(const MatrixType& mat, const Rhs& rhs, Dest& x,
-                    const Preconditioner& precond, int& iters,
+                    const Preconditioner& precond, Index& iters,
                     typename Dest::RealScalar& tol_error)
         {
             using std::sqrt;
@@ -48,8 +48,8 @@ namespace Eigen {
             }
             
             // initialize
-            const int maxIters(iters);  // initialize maxIters to iters
-            const int N(mat.cols());    // the size of the matrix
+            const Index maxIters(iters);  // initialize maxIters to iters
+            const Index N(mat.cols());    // the size of the matrix
             const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
             
             // Initialize preconditioned Lanczos
@@ -144,7 +144,6 @@ namespace Eigen {
     
     template< typename _MatrixType, int _UpLo=Lower,
     typename _Preconditioner = IdentityPreconditioner>
-//    typename _Preconditioner = IdentityPreconditioner<typename _MatrixType::Scalar> > // preconditioner must be positive definite
     class MINRES;
     
     namespace internal {
@@ -166,8 +165,8 @@ namespace Eigen {
      * The vectors x and b can be either dense or sparse.
      *
      * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
-     * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
-     *               or Upper. Default is Lower.
+     * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower,
+     *               Upper, or Lower|Upper in which the full matrix entries will be considered. Default is Lower.
      * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
      *
      * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
@@ -192,6 +191,8 @@ namespace Eigen {
      * By default the iterations start with x=0 as an initial guess of the solution.
      * One can control the start using the solveWithGuess() method.
      *
+     * MINRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
+     *
      * \sa class ConjugateGradient, BiCGSTAB, SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
      */
     template< typename _MatrixType, int _UpLo, typename _Preconditioner>
@@ -199,15 +200,15 @@ namespace Eigen {
     {
         
         typedef IterativeSolverBase<MINRES> Base;
-        using Base::mp_matrix;
+        using Base::matrix;
         using Base::m_error;
         using Base::m_iterations;
         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;
         typedef typename MatrixType::RealScalar RealScalar;
         typedef _Preconditioner Preconditioner;
         
@@ -233,42 +234,36 @@ namespace Eigen {
         
         /** Destructor. */
         ~MINRES(){}
-		
-        /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
-         * \a x0 as an initial solution.
-         *
-         * \sa compute()
-         */
-        template<typename Rhs,typename Guess>
-        inline const internal::solve_retval_with_guess<MINRES, Rhs, Guess>
-        solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
-        {
-            eigen_assert(m_isInitialized && "MINRES is not initialized.");
-            eigen_assert(Base::rows()==b.rows()
-                         && "MINRES::solve(): invalid number of rows of the right hand side matrix b");
-            return internal::solve_retval_with_guess
-            <MINRES, Rhs, Guess>(*this, b.derived(), x0);
-        }
-        
+
         /** \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
         {
+            typedef typename Base::MatrixWrapper MatrixWrapper;
+            typedef typename Base::ActualMatrixType ActualMatrixType;
+            enum {
+              TransposeInput  =   (!MatrixWrapper::MatrixFree)
+                              &&  (UpLo==(Lower|Upper))
+                              &&  (!MatrixType::IsRowMajor)
+                              &&  (!NumTraits<Scalar>::IsComplex)
+            };
+            typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type RowMajorWrapper;
+            EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
             typedef typename internal::conditional<UpLo==(Lower|Upper),
-                                                   const MatrixType&,
-                                                   SparseSelfAdjointView<const MatrixType, UpLo>
-                                                  >::type MatrixWrapperType;
-                                          
+                                                  RowMajorWrapper,
+                                                  typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
+                                            >::type SelfAdjointWrapper;
+
             m_iterations = Base::maxIterations();
             m_error = Base::m_tolerance;
-            
+            RowMajorWrapper row_mat(matrix());
             for(int j=0; j<b.cols(); ++j)
             {
                 m_iterations = Base::maxIterations();
                 m_error = Base::m_tolerance;
                 
                 typename Dest::ColXpr xj(x,j);
-                internal::minres(MatrixWrapperType(*mp_matrix), b.col(j), xj,
+                internal::minres(SelfAdjointWrapper(row_mat), b.col(j), xj,
                                  Base::m_preconditioner, m_iterations, m_error);
             }
             
@@ -278,33 +273,16 @@ 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:
         
     };
-    
-    namespace internal {
-        
-        template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
-        struct solve_retval<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs>
-        : solve_retval_base<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs>
-        {
-            typedef MINRES<_MatrixType,_UpLo,_Preconditioner> Dec;
-            EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-            
-            template<typename Dest> void evalTo(Dest& dst) const
-            {
-                dec()._solve(rhs(),dst);
-            }
-        };
-        
-    } // end namespace internal
-    
+
 } // end namespace Eigen
 
 #endif // EIGEN_MINRES_H