Allows Lower|Upper as a template argument of CG and MINRES: in this case the full matrix will be considered.

5 files changed, 23 insertions, 8 deletions

diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
index 3e024bda1..4857dd9e9 100644
--- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
+++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
@@ -113,8 +113,8 @@ struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
   * The matrix A must be selfadjoint. The matrix A and the vectors x and b can be either dense or sparse.
   *
   * \tparam _MatrixType the type of the 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()
@@ -197,6 +197,10 @@ public:
   template<typename Rhs,typename Dest>
   void _solve_with_guess_impl(const Rhs& b, Dest& x) const
   {
+    typedef typename internal::conditional<UpLo==(Lower|Upper),
+                                           Ref<const MatrixType>&,
+                                           SparseSelfAdjointView<const Ref<const MatrixType>, UpLo>
+                                          >::type MatrixWrapperType;
     m_iterations = Base::maxIterations();
     m_error = Base::m_tolerance;
 
@@ -206,8 +210,7 @@ public:
       m_error = Base::m_tolerance;
 
       typename Dest::ColXpr xj(x,j);
-      internal::conjugate_gradient(mp_matrix.template selfadjointView<UpLo>(), b.col(j), xj,
-                                   Base::m_preconditioner, m_iterations, m_error);
+      internal::conjugate_gradient(MatrixWrapperType(mp_matrix), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
     }
 
     m_isInitialized = true;
diff --git a/test/conjugate_gradient.cpp b/test/conjugate_gradient.cpp
index 869051b31..019cc4d64 100644
--- a/test/conjugate_gradient.cpp
+++ b/test/conjugate_gradient.cpp
@@ -12,13 +12,15 @@
 
 template<typename T> void test_conjugate_gradient_T()
 {
-  ConjugateGradient<SparseMatrix<T>, Lower> cg_colmajor_lower_diag;
-  ConjugateGradient<SparseMatrix<T>, Upper> cg_colmajor_upper_diag;
+  ConjugateGradient<SparseMatrix<T>, Lower      > cg_colmajor_lower_diag;
+  ConjugateGradient<SparseMatrix<T>, Upper      > cg_colmajor_upper_diag;
+  ConjugateGradient<SparseMatrix<T>, Lower|Upper> cg_colmajor_loup_diag;
   ConjugateGradient<SparseMatrix<T>, Lower, IdentityPreconditioner> cg_colmajor_lower_I;
   ConjugateGradient<SparseMatrix<T>, Upper, IdentityPreconditioner> cg_colmajor_upper_I;
 
   CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_lower_diag)  );
   CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_upper_diag)  );
+  CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_loup_diag)   );
   CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_lower_I)     );
   CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_upper_I)     );
 }
diff --git a/test/sparse_solver.h b/test/sparse_solver.h
index 887ff02a9..6cf99e0b0 100644
--- a/test/sparse_solver.h
+++ b/test/sparse_solver.h
@@ -195,7 +195,10 @@ int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typena
   dA = dM * dM.adjoint();
   
   halfA.resize(size,size);
-  halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M);
+  if(Solver::UpLo==(Lower|Upper))
+    halfA = A;
+  else
+    halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M);
   
   return size;
 }
diff --git a/unsupported/Eigen/src/IterativeSolvers/MINRES.h b/unsupported/Eigen/src/IterativeSolvers/MINRES.h
index a34902001..65cffc255 100644
--- a/unsupported/Eigen/src/IterativeSolvers/MINRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/MINRES.h
@@ -250,6 +250,11 @@ namespace Eigen {
         template<typename Rhs,typename Dest>
         void _solve_with_guess_impl(const Rhs& b, Dest& x) const
         {
+            typedef typename internal::conditional<UpLo==(Lower|Upper),
+                                                   Ref<const MatrixType>&,
+                                                   SparseSelfAdjointView<const Ref<const MatrixType>, UpLo>
+                                                  >::type MatrixWrapperType;
+                                          
             m_iterations = Base::maxIterations();
             m_error = Base::m_tolerance;
             
@@ -259,7 +264,7 @@ namespace Eigen {
                 m_error = Base::m_tolerance;
                 
                 typename Dest::ColXpr xj(x,j);
-                internal::minres(mp_matrix.template selfadjointView<UpLo>(), b.col(j), xj,
+                internal::minres(MatrixWrapperType(mp_matrix), b.col(j), xj,
                                  Base::m_preconditioner, m_iterations, m_error);
             }
             
diff --git a/unsupported/test/minres.cpp b/unsupported/test/minres.cpp
index 81b762c37..f6e526bbd 100644
--- a/unsupported/test/minres.cpp
+++ b/unsupported/test/minres.cpp
@@ -21,6 +21,7 @@ template<typename T> void test_minres_T()
   // Diagonal preconditioner
   MINRES<SparseMatrix<T>, Lower, DiagonalPreconditioner<T> > minres_colmajor_lower_diag;
   MINRES<SparseMatrix<T>, Upper, DiagonalPreconditioner<T> > minres_colmajor_upper_diag;
+  MINRES<SparseMatrix<T>, Upper, DiagonalPreconditioner<T> > minres_colmajor_uplo_diag;
   
   // call tests for SPD matrix
   CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_lower_I) );
@@ -28,6 +29,7 @@ template<typename T> void test_minres_T()
     
   CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_lower_diag)  );
   CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_upper_diag)  );
+  CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_uplo_diag)  );
     
   // TO DO: symmetric semi-definite matrix
   // TO DO: symmetric indefinite matrix