Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/unsupported/Eigen/src/IterativeSolvers/CMakeLists.txt b/unsupported/Eigen/src/IterativeSolvers/CMakeLists.txt
deleted file mode 100644
index 7986afc..0000000
--- a/unsupported/Eigen/src/IterativeSolvers/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_IterativeSolvers_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_IterativeSolvers_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/IterativeSolvers COMPONENT Devel
- )
diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
index 9fcc8a8..4079e23 100644
--- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
@@ -39,8 +39,6 @@
void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut)
{
eigen_assert(vec.size() == perm.size());
- typedef typename IndexType::Scalar Index;
- typedef typename VectorType::Scalar Scalar;
bool flag;
for (Index k = 0; k < ncut; k++)
{
@@ -84,6 +82,8 @@
* x = solver.solve(b);
* \endcode
*
+ * DGMRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
+ *
* References :
* [1] D. NUENTSA WAKAM and F. PACULL, Memory Efficient Hybrid
* Algebraic Solvers for Linear Systems Arising from Compressible
@@ -101,16 +101,17 @@
class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
{
typedef IterativeSolverBase<DGMRES> Base;
- using Base::mp_matrix;
+ using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
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;
+ typedef typename MatrixType::StorageIndex StorageIndex;
typedef typename MatrixType::RealScalar RealScalar;
typedef _Preconditioner Preconditioner;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
@@ -133,39 +134,23 @@
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
- DGMRES(const MatrixType& A) : Base(A),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false)
- {}
+ template<typename MatrixDerived>
+ explicit DGMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
~DGMRES() {}
- /** \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<DGMRES, Rhs, Guess>
- solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
- {
- eigen_assert(m_isInitialized && "DGMRES is not initialized.");
- eigen_assert(Base::rows()==b.rows()
- && "DGMRES::solve(): invalid number of rows of the right hand side matrix b");
- return internal::solve_retval_with_guess
- <DGMRES, 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
{
bool failed = false;
- for(int j=0; j<b.cols(); ++j)
+ for(Index j=0; j<b.cols(); ++j)
{
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
typename Dest::ColXpr xj(x,j);
- dgmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner);
+ dgmres(matrix(), b.col(j), xj, Base::m_preconditioner);
}
m_info = failed ? NumericalIssue
: m_error <= Base::m_tolerance ? Success
@@ -175,25 +160,25 @@
/** \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
*/
- int restart() { return m_restart; }
+ Index restart() { return m_restart; }
/**
* Set the restart value (default is 30)
*/
- void set_restart(const int restart) { m_restart=restart; }
+ void set_restart(const Index restart) { m_restart=restart; }
/**
* Set the number of eigenvalues to deflate at each restart
*/
- void setEigenv(const int neig)
+ void setEigenv(const Index neig)
{
m_neig = neig;
if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
@@ -202,12 +187,12 @@
/**
* Get the size of the deflation subspace size
*/
- int deflSize() {return m_r; }
+ Index deflSize() {return m_r; }
/**
* Set the maximum size of the deflation subspace
*/
- void setMaxEigenv(const int maxNeig) { m_maxNeig = maxNeig; }
+ void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; }
protected:
// DGMRES algorithm
@@ -215,12 +200,12 @@
void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
// Perform one cycle of GMRES
template<typename Dest>
- int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const;
+ Index dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const;
// Compute data to use for deflation
- int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const;
+ Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const;
// Apply deflation to a vector
template<typename RhsType, typename DestType>
- int dgmresApplyDeflation(const RhsType& In, DestType& Out) const;
+ Index dgmresApplyDeflation(const RhsType& In, DestType& Out) const;
ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
// Init data for deflation
@@ -233,9 +218,9 @@
mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
- mutable int m_neig; //Number of eigenvalues to extract at each restart
- mutable int m_r; // Current number of deflated eigenvalues, size of m_U
- mutable int m_maxNeig; // Maximum number of eigenvalues to deflate
+ mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart
+ mutable Index m_r; // Current number of deflated eigenvalues, size of m_U
+ mutable Index m_maxNeig; // Maximum number of eigenvalues to deflate
mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
mutable bool m_isDeflAllocated;
mutable bool m_isDeflInitialized;
@@ -257,9 +242,9 @@
const Preconditioner& precond) const
{
//Initialization
- int n = mat.rows();
+ Index n = mat.rows();
DenseVector r0(n);
- int nbIts = 0;
+ Index nbIts = 0;
m_H.resize(m_restart+1, m_restart);
m_Hes.resize(m_restart, m_restart);
m_V.resize(n,m_restart+1);
@@ -297,7 +282,7 @@
*/
template< typename _MatrixType, typename _Preconditioner>
template<typename Dest>
-int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const
+Index DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const
{
//Initialization
DenseVector g(m_restart+1); // Right hand side of the least square problem
@@ -306,8 +291,8 @@
m_V.col(0) = r0/beta;
m_info = NoConvergence;
std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
- int it = 0; // Number of inner iterations
- int n = mat.rows();
+ Index it = 0; // Number of inner iterations
+ Index n = mat.rows();
DenseVector tv1(n), tv2(n); //Temporary vectors
while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
{
@@ -325,7 +310,7 @@
// Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
Scalar coef;
- for (int i = 0; i <= it; ++i)
+ for (Index i = 0; i <= it; ++i)
{
coef = tv1.dot(m_V.col(i));
tv1 = tv1 - coef * m_V.col(i);
@@ -341,7 +326,7 @@
// FIXME Check for happy breakdown
// Update Hessenberg matrix with Givens rotations
- for (int i = 1; i <= it; ++i)
+ for (Index i = 1; i <= it; ++i)
{
m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
}
@@ -353,7 +338,7 @@
beta = std::abs(g(it+1));
m_error = beta/normRhs;
- std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
+ // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
it++; nbIts++;
if (m_error < m_tolerance)
@@ -407,7 +392,6 @@
template< typename _MatrixType, typename _Preconditioner>
inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const
{
- typedef typename MatrixType::Index Index;
const DenseMatrix& T = schurofH.matrixT();
Index it = T.rows();
ComplexVector eig(it);
@@ -431,7 +415,7 @@
}
template< typename _MatrixType, typename _Preconditioner>
-int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const
+Index DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const
{
// First, find the Schur form of the Hessenberg matrix H
typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH;
@@ -441,13 +425,13 @@
schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU);
ComplexVector eig(it);
- Matrix<Index,Dynamic,1>perm(it);
+ Matrix<StorageIndex,Dynamic,1>perm(it);
eig = this->schurValues(schurofH);
// Reorder the absolute values of Schur values
DenseRealVector modulEig(it);
- for (int j=0; j<it; ++j) modulEig(j) = std::abs(eig(j));
- perm.setLinSpaced(it,0,it-1);
+ for (Index j=0; j<it; ++j) modulEig(j) = std::abs(eig(j));
+ perm.setLinSpaced(it,0,internal::convert_index<StorageIndex>(it-1));
internal::sortWithPermutation(modulEig, perm, neig);
if (!m_lambdaN)
@@ -455,7 +439,7 @@
m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
}
//Count the real number of extracted eigenvalues (with complex conjugates)
- int nbrEig = 0;
+ Index nbrEig = 0;
while (nbrEig < neig)
{
if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++;
@@ -464,7 +448,7 @@
// Extract the Schur vectors corresponding to the smallest Ritz values
DenseMatrix Sr(it, nbrEig);
Sr.setZero();
- for (int j = 0; j < nbrEig; j++)
+ for (Index j = 0; j < nbrEig; j++)
{
Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
}
@@ -475,8 +459,8 @@
if (m_r)
{
// Orthogonalize X against m_U using modified Gram-Schmidt
- for (int j = 0; j < nbrEig; j++)
- for (int k =0; k < m_r; k++)
+ for (Index j = 0; j < nbrEig; j++)
+ for (Index k =0; k < m_r; k++)
X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k);
}
@@ -486,7 +470,7 @@
dgmresInitDeflation(m);
DenseMatrix MX(m, nbrEig);
DenseVector tv1(m);
- for (int j = 0; j < nbrEig; j++)
+ for (Index j = 0; j < nbrEig; j++)
{
tv1 = mat * X.col(j);
MX.col(j) = precond.solve(tv1);
@@ -501,8 +485,8 @@
}
// Save X into m_U and m_MX in m_MU
- for (int j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
- for (int j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
+ for (Index j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
+ for (Index j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
// Increase the size of the invariant subspace
m_r += nbrEig;
@@ -515,28 +499,12 @@
}
template<typename _MatrixType, typename _Preconditioner>
template<typename RhsType, typename DestType>
-int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
+Index DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
{
DenseVector x1 = m_U.leftCols(m_r).transpose() * x;
y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1);
return 0;
}
-namespace internal {
-
- template<typename _MatrixType, typename _Preconditioner, typename Rhs>
-struct solve_retval<DGMRES<_MatrixType, _Preconditioner>, Rhs>
- : solve_retval_base<DGMRES<_MatrixType, _Preconditioner>, Rhs>
-{
- typedef DGMRES<_MatrixType, _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
diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
index 7ba13af..92618b1 100644
--- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
@@ -11,193 +11,197 @@
#ifndef EIGEN_GMRES_H
#define EIGEN_GMRES_H
-namespace Eigen {
+namespace Eigen {
namespace internal {
/**
- * Generalized Minimal Residual Algorithm based on the
- * Arnoldi algorithm implemented with Householder reflections.
- *
- * Parameters:
- * \param mat matrix of linear system of equations
- * \param Rhs right hand side vector of linear system of equations
- * \param x on input: initial guess, on output: solution
- * \param precond preconditioner used
- * \param iters on input: maximum number of iterations to perform
- * on output: number of iterations performed
- * \param restart number of iterations for a restart
- * \param tol_error on input: residual tolerance
- * on output: residuum achieved
- *
- * \sa IterativeMethods::bicgstab()
- *
- *
- * For references, please see:
- *
- * Saad, Y. and Schultz, M. H.
- * GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems.
- * SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869.
- *
- * Saad, Y.
- * Iterative Methods for Sparse Linear Systems.
- * Society for Industrial and Applied Mathematics, Philadelphia, 2003.
- *
- * Walker, H. F.
- * Implementations of the GMRES method.
- * Comput.Phys.Comm. 53, 1989, pp. 311 - 320.
- *
- * Walker, H. F.
- * Implementation of the GMRES Method using Householder Transformations.
- * SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163.
- *
- */
+* Generalized Minimal Residual Algorithm based on the
+* Arnoldi algorithm implemented with Householder reflections.
+*
+* Parameters:
+* \param mat matrix of linear system of equations
+* \param rhs right hand side vector of linear system of equations
+* \param x on input: initial guess, on output: solution
+* \param precond preconditioner used
+* \param iters on input: maximum number of iterations to perform
+* on output: number of iterations performed
+* \param restart number of iterations for a restart
+* \param tol_error on input: relative residual tolerance
+* on output: residuum achieved
+*
+* \sa IterativeMethods::bicgstab()
+*
+*
+* For references, please see:
+*
+* Saad, Y. and Schultz, M. H.
+* GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems.
+* SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869.
+*
+* Saad, Y.
+* Iterative Methods for Sparse Linear Systems.
+* Society for Industrial and Applied Mathematics, Philadelphia, 2003.
+*
+* Walker, H. F.
+* Implementations of the GMRES method.
+* Comput.Phys.Comm. 53, 1989, pp. 311 - 320.
+*
+* Walker, H. F.
+* Implementation of the GMRES Method using Householder Transformations.
+* SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163.
+*
+*/
template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Preconditioner & precond,
- int &iters, const int &restart, typename Dest::RealScalar & tol_error) {
+ Index &iters, const Index &restart, typename Dest::RealScalar & tol_error) {
- using std::sqrt;
- using std::abs;
+ using std::sqrt;
+ using std::abs;
- typedef typename Dest::RealScalar RealScalar;
- typedef typename Dest::Scalar Scalar;
- typedef Matrix < Scalar, Dynamic, 1 > VectorType;
- typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType;
+ typedef typename Dest::RealScalar RealScalar;
+ typedef typename Dest::Scalar Scalar;
+ typedef Matrix < Scalar, Dynamic, 1 > VectorType;
+ typedef Matrix < Scalar, Dynamic, Dynamic, ColMajor> FMatrixType;
- RealScalar tol = tol_error;
- const int maxIters = iters;
- iters = 0;
+ RealScalar tol = tol_error;
+ const Index maxIters = iters;
+ iters = 0;
- const int m = mat.rows();
+ const Index m = mat.rows();
- VectorType p0 = rhs - mat*x;
- VectorType r0 = precond.solve(p0);
-
- // is initial guess already good enough?
- if(abs(r0.norm()) < tol) {
- return true;
- }
+ // residual and preconditioned residual
+ VectorType p0 = rhs - mat*x;
+ VectorType r0 = precond.solve(p0);
- VectorType w = VectorType::Zero(restart + 1);
+ const RealScalar r0Norm = r0.norm();
- FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix
- VectorType tau = VectorType::Zero(restart + 1);
- std::vector < JacobiRotation < Scalar > > G(restart);
+ // is initial guess already good enough?
+ if(r0Norm == 0)
+ {
+ tol_error = 0;
+ return true;
+ }
- // generate first Householder vector
- VectorType e(m-1);
- RealScalar beta;
- r0.makeHouseholder(e, tau.coeffRef(0), beta);
- w(0)=(Scalar) beta;
- H.bottomLeftCorner(m - 1, 1) = e;
+ // storage for Hessenberg matrix and Householder data
+ FMatrixType H = FMatrixType::Zero(m, restart + 1);
+ VectorType w = VectorType::Zero(restart + 1);
+ VectorType tau = VectorType::Zero(restart + 1);
- for (int k = 1; k <= restart; ++k) {
+ // storage for Jacobi rotations
+ std::vector < JacobiRotation < Scalar > > G(restart);
+
+ // storage for temporaries
+ VectorType t(m), v(m), workspace(m), x_new(m);
- ++iters;
+ // generate first Householder vector
+ Ref<VectorType> H0_tail = H.col(0).tail(m - 1);
+ RealScalar beta;
+ r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta);
+ w(0) = Scalar(beta);
+
+ for (Index k = 1; k <= restart; ++k)
+ {
+ ++iters;
- VectorType v = VectorType::Unit(m, k - 1), workspace(m);
+ v = VectorType::Unit(m, k - 1);
- // apply Householder reflections H_{1} ... H_{k-1} to v
- for (int i = k - 1; i >= 0; --i) {
- v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
- }
+ // apply Householder reflections H_{1} ... H_{k-1} to v
+ // TODO: use a HouseholderSequence
+ for (Index i = k - 1; i >= 0; --i) {
+ v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
+ }
- // apply matrix M to v: v = mat * v;
- VectorType t=mat*v;
- v=precond.solve(t);
+ // apply matrix M to v: v = mat * v;
+ t.noalias() = mat * v;
+ v = precond.solve(t);
- // apply Householder reflections H_{k-1} ... H_{1} to v
- for (int i = 0; i < k; ++i) {
- v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
- }
+ // apply Householder reflections H_{k-1} ... H_{1} to v
+ // TODO: use a HouseholderSequence
+ for (Index i = 0; i < k; ++i) {
+ v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
+ }
- if (v.tail(m - k).norm() != 0.0) {
+ if (v.tail(m - k).norm() != 0.0)
+ {
+ if (k <= restart)
+ {
+ // generate new Householder vector
+ Ref<VectorType> Hk_tail = H.col(k).tail(m - k - 1);
+ v.tail(m - k).makeHouseholder(Hk_tail, tau.coeffRef(k), beta);
- if (k <= restart) {
+ // apply Householder reflection H_{k} to v
+ v.tail(m - k).applyHouseholderOnTheLeft(Hk_tail, tau.coeffRef(k), workspace.data());
+ }
+ }
- // generate new Householder vector
- VectorType e(m - k - 1);
- RealScalar beta;
- v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta);
- H.col(k).tail(m - k - 1) = e;
+ if (k > 1)
+ {
+ for (Index i = 0; i < k - 1; ++i)
+ {
+ // apply old Givens rotations to v
+ v.applyOnTheLeft(i, i + 1, G[i].adjoint());
+ }
+ }
- // apply Householder reflection H_{k} to v
- v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data());
+ if (k<m && v(k) != (Scalar) 0)
+ {
+ // determine next Givens rotation
+ G[k - 1].makeGivens(v(k - 1), v(k));
- }
- }
+ // apply Givens rotation to v and w
+ v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
+ w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
+ }
- if (k > 1) {
- for (int i = 0; i < k - 1; ++i) {
- // apply old Givens rotations to v
- v.applyOnTheLeft(i, i + 1, G[i].adjoint());
- }
- }
+ // insert coefficients into upper matrix triangle
+ H.col(k-1).head(k) = v.head(k);
- if (k<m && v(k) != (Scalar) 0) {
- // determine next Givens rotation
- G[k - 1].makeGivens(v(k - 1), v(k));
+ tol_error = abs(w(k)) / r0Norm;
+ bool stop = (k==m || tol_error < tol || iters == maxIters);
- // apply Givens rotation to v and w
- v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
- w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
+ if (stop || k == restart)
+ {
+ // solve upper triangular system
+ Ref<VectorType> y = w.head(k);
+ H.topLeftCorner(k, k).template triangularView <Upper>().solveInPlace(y);
- }
+ // use Horner-like scheme to calculate solution vector
+ x_new.setZero();
+ for (Index i = k - 1; i >= 0; --i)
+ {
+ x_new(i) += y(i);
+ // apply Householder reflection H_{i} to x_new
+ x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
+ }
- // insert coefficients into upper matrix triangle
- H.col(k - 1).head(k) = v.head(k);
+ x += x_new;
- bool stop=(k==m || abs(w(k)) < tol || iters == maxIters);
+ if(stop)
+ {
+ return true;
+ }
+ else
+ {
+ k=0;
- if (stop || k == restart) {
+ // reset data for restart
+ p0.noalias() = rhs - mat*x;
+ r0 = precond.solve(p0);
- // solve upper triangular system
- VectorType y = w.head(k);
- H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y);
+ // clear Hessenberg matrix and Householder data
+ H.setZero();
+ w.setZero();
+ tau.setZero();
- // use Horner-like scheme to calculate solution vector
- VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
+ // generate first Householder vector
+ r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta);
+ w(0) = Scalar(beta);
+ }
+ }
+ }
- // apply Householder reflection H_{k} to x_new
- x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data());
-
- for (int i = k - 2; i >= 0; --i) {
- x_new += y(i) * VectorType::Unit(m, i);
- // apply Householder reflection H_{i} to x_new
- x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
- }
-
- x += x_new;
-
- if (stop) {
- return true;
- } else {
- k=0;
-
- // reset data for a restart r0 = rhs - mat * x;
- VectorType p0=mat*x;
- VectorType p1=precond.solve(p0);
- r0 = rhs - p1;
-// r0_sqnorm = r0.squaredNorm();
- w = VectorType::Zero(restart + 1);
- H = FMatrixType::Zero(m, restart + 1);
- tau = VectorType::Zero(restart + 1);
-
- // generate first Householder vector
- RealScalar beta;
- r0.makeHouseholder(e, tau.coeffRef(0), beta);
- w(0)=(Scalar) beta;
- H.bottomLeftCorner(m - 1, 1) = e;
-
- }
-
- }
-
-
-
- }
-
- return false;
+ return false;
}
@@ -230,7 +234,7 @@
* The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
* and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
* and NumTraits<Scalar>::epsilon() for the tolerance.
- *
+ *
* This class can be used as the direct solver classes. Here is a typical usage example:
* \code
* int n = 10000;
@@ -244,29 +248,31 @@
* // update b, and solve again
* x = solver.solve(b);
* \endcode
- *
+ *
* By default the iterations start with x=0 as an initial guess of the solution.
* One can control the start using the solveWithGuess() method.
*
+ * GMRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
+ *
* \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
*/
template< typename _MatrixType, typename _Preconditioner>
class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> >
{
typedef IterativeSolverBase<GMRES> Base;
- using Base::mp_matrix;
+ using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
using Base::m_isInitialized;
-
+
private:
- int m_restart;
-
+ Index m_restart;
+
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;
@@ -276,95 +282,62 @@
GMRES() : Base(), m_restart(30) {}
/** Initialize the solver with matrix \a A for further \c Ax=b solving.
- *
+ *
* This constructor is a shortcut for the default constructor followed
* by a call to compute().
- *
+ *
* \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.
*/
- GMRES(const MatrixType& A) : Base(A), m_restart(30) {}
+ template<typename MatrixDerived>
+ explicit GMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30) {}
~GMRES() {}
-
+
/** Get the number of iterations after that a restart is performed.
*/
- int get_restart() { return m_restart; }
-
+ Index get_restart() { return m_restart; }
+
/** Set the number of iterations after that a restart is performed.
* \param restart number of iterations for a restarti, default is 30.
*/
- void set_restart(const int restart) { m_restart=restart; }
-
- /** \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<GMRES, Rhs, Guess>
- solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
- {
- eigen_assert(m_isInitialized && "GMRES is not initialized.");
- eigen_assert(Base::rows()==b.rows()
- && "GMRES::solve(): invalid number of rows of the right hand side matrix b");
- return internal::solve_retval_with_guess
- <GMRES, Rhs, Guess>(*this, b.derived(), x0);
- }
-
+ void set_restart(const Index restart) { m_restart=restart; }
+
/** \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)
+ for(Index j=0; j<b.cols(); ++j)
{
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
-
+
typename Dest::ColXpr xj(x,j);
- if(!internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
+ if(!internal::gmres(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
failed = true;
}
m_info = failed ? NumericalIssue
- : m_error <= Base::m_tolerance ? Success
- : NoConvergence;
+ : m_error <= Base::m_tolerance ? Success
+ : NoConvergence;
m_isInitialized = true;
}
/** \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:
};
-
-namespace internal {
-
- template<typename _MatrixType, typename _Preconditioner, typename Rhs>
-struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs>
- : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs>
-{
- typedef GMRES<_MatrixType, _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_GMRES_H
diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
deleted file mode 100644
index 661c1f2..0000000
--- a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
+++ /dev/null
@@ -1,278 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_INCOMPLETE_CHOlESKY_H
-#define EIGEN_INCOMPLETE_CHOlESKY_H
-#include "Eigen/src/IterativeLinearSolvers/IncompleteLUT.h"
-#include <Eigen/OrderingMethods>
-#include <list>
-
-namespace Eigen {
-/**
- * \brief Modified Incomplete Cholesky with dual threshold
- *
- * References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
- * Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999
- *
- * \tparam _MatrixType The type of the sparse matrix. It should be a symmetric
- * matrix. It is advised to give a row-oriented sparse matrix
- * \tparam _UpLo The triangular part of the matrix to reference.
- * \tparam _OrderingType
- */
-
-template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> >
-class IncompleteCholesky : internal::noncopyable
-{
- public:
- typedef SparseMatrix<Scalar,ColMajor> MatrixType;
- typedef _OrderingType OrderingType;
- typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
- typedef Matrix<Scalar,Dynamic,1> ScalarType;
- typedef Matrix<Index,Dynamic, 1> IndexType;
- typedef std::vector<std::list<Index> > VectorList;
- enum { UpLo = _UpLo };
- public:
- IncompleteCholesky() : m_shift(1),m_factorizationIsOk(false) {}
- IncompleteCholesky(const MatrixType& matrix) : m_shift(1),m_factorizationIsOk(false)
- {
- compute(matrix);
- }
-
- Index rows() const { return m_L.rows(); }
-
- Index cols() const { return m_L.cols(); }
-
-
- /** \brief Reports whether previous computation was successful.
- *
- * \returns \c Success if computation was succesful,
- * \c NumericalIssue if the matrix appears to be negative.
- */
- ComputationInfo info() const
- {
- eigen_assert(m_isInitialized && "IncompleteLLT is not initialized.");
- return m_info;
- }
-
- /**
- * \brief Set the initial shift parameter
- */
- void setShift( Scalar shift) { m_shift = shift; }
-
- /**
- * \brief Computes the fill reducing permutation vector.
- */
- template<typename MatrixType>
- void analyzePattern(const MatrixType& mat)
- {
- OrderingType ord;
- ord(mat.template selfadjointView<UpLo>(), m_perm);
- m_analysisIsOk = true;
- }
-
- template<typename MatrixType>
- void factorize(const MatrixType& amat);
-
- template<typename MatrixType>
- void compute (const MatrixType& matrix)
- {
- analyzePattern(matrix);
- factorize(matrix);
- }
-
- template<typename Rhs, typename Dest>
- void _solve(const Rhs& b, Dest& x) const
- {
- eigen_assert(m_factorizationIsOk && "factorize() should be called first");
- if (m_perm.rows() == b.rows())
- x = m_perm.inverse() * b;
- else
- x = b;
- x = m_scal.asDiagonal() * x;
- x = m_L.template triangularView<UnitLower>().solve(x);
- x = m_L.adjoint().template triangularView<Upper>().solve(x);
- if (m_perm.rows() == b.rows())
- x = m_perm * x;
- x = m_scal.asDiagonal() * x;
- }
- template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs>
- solve(const MatrixBase<Rhs>& b) const
- {
- eigen_assert(m_factorizationIsOk && "IncompleteLLT did not succeed");
- eigen_assert(m_isInitialized && "IncompleteLLT is not initialized.");
- eigen_assert(cols()==b.rows()
- && "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b");
- return internal::solve_retval<IncompleteCholesky, Rhs>(*this, b.derived());
- }
- 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;
-
- private:
- template <typename IdxType, typename SclType>
- inline void updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol);
-};
-
-template<typename Scalar, int _UpLo, typename OrderingType>
-template<typename _MatrixType>
-void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
-{
- using std::sqrt;
- using std::min;
- eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
-
- // Dropping strategies : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
-
- // Apply the fill-reducing permutation computed in analyzePattern()
- if (m_perm.rows() == mat.rows() ) // To detect the null permutation
- m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
- else
- m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
-
- Index n = m_L.cols();
- Index nnz = m_L.nonZeros();
- Map<ScalarType> vals(m_L.valuePtr(), nnz); //values
- Map<IndexType> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices
- Map<IndexType> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
- IndexType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
- VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
- ScalarType curCol(n); // Store a nonzero values in each column
- IndexType irow(n); // Row indices of nonzero elements in each column
-
-
- // Computes the scaling factors
- m_scal.resize(n);
- for (int j = 0; j < n; j++)
- {
- m_scal(j) = m_L.col(j).norm();
- m_scal(j) = sqrt(m_scal(j));
- }
- // Scale and compute the shift for the matrix
- Scalar mindiag = vals[0];
- for (int j = 0; j < n; j++){
- for (int k = colPtr[j]; k < colPtr[j+1]; k++)
- vals[k] /= (m_scal(j) * m_scal(rowIdx[k]));
- mindiag = (min)(vals[colPtr[j]], mindiag);
- }
-
- if(mindiag < Scalar(0.)) m_shift = m_shift - mindiag;
- // Apply the shift to the diagonal elements of the matrix
- for (int j = 0; j < n; j++)
- vals[colPtr[j]] += m_shift;
- // jki version of the Cholesky factorization
- for (int j=0; j < n; ++j)
- {
- //Left-looking factorize the column j
- // First, load the jth column into curCol
- Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
- curCol.setZero();
- irow.setLinSpaced(n,0,n-1);
- for (int i = colPtr[j] + 1; i < colPtr[j+1]; i++)
- {
- curCol(rowIdx[i]) = vals[i];
- irow(rowIdx[i]) = rowIdx[i];
- }
- std::list<int>::iterator k;
- // Browse all previous columns that will update column j
- for(k = listCol[j].begin(); k != listCol[j].end(); k++)
- {
- int jk = firstElt(*k); // First element to use in the column
- jk += 1;
- for (int i = jk; i < colPtr[*k+1]; i++)
- {
- curCol(rowIdx[i]) -= vals[i] * vals[jk] ;
- }
- updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
- }
-
- // Scale the current column
- if(RealScalar(diag) <= 0)
- {
- std::cerr << "\nNegative diagonal during Incomplete factorization... "<< j << "\n";
- m_info = NumericalIssue;
- return;
- }
- RealScalar rdiag = sqrt(RealScalar(diag));
- vals[colPtr[j]] = rdiag;
- for (int i = j+1; i < n; i++)
- {
- //Scale
- curCol(i) /= rdiag;
- //Update the remaining diagonals with curCol
- vals[colPtr[i]] -= curCol(i) * curCol(i);
- }
- // Select the largest p elements
- // p is the original number of elements in the column (without the diagonal)
- int p = colPtr[j+1] - colPtr[j] - 1 ;
- internal::QuickSplit(curCol, irow, p);
- // Insert the largest p elements in the matrix
- int cpt = 0;
- for (int i = colPtr[j]+1; i < colPtr[j+1]; i++)
- {
- vals[i] = curCol(cpt);
- rowIdx[i] = irow(cpt);
- cpt ++;
- }
- // Get the first smallest row index and put it after the diagonal element
- Index jk = colPtr(j)+1;
- updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
- }
- m_factorizationIsOk = true;
- m_isInitialized = true;
- m_info = Success;
-}
-
-template<typename Scalar, int _UpLo, typename OrderingType>
-template <typename IdxType, typename SclType>
-inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol)
-{
- if (jk < colPtr(col+1) )
- {
- Index p = colPtr(col+1) - jk;
- Index minpos;
- rowIdx.segment(jk,p).minCoeff(&minpos);
- minpos += jk;
- if (rowIdx(minpos) != rowIdx(jk))
- {
- //Swap
- std::swap(rowIdx(jk),rowIdx(minpos));
- std::swap(vals(jk),vals(minpos));
- }
- firstElt(col) = jk;
- listCol[rowIdx(jk)].push_back(col);
- }
-}
-namespace internal {
-
-template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs>
-struct solve_retval<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs>
- : solve_retval_base<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs>
-{
- typedef IncompleteCholesky<_Scalar, _UpLo, OrderingType> 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
diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h
index 67e7801..7d08c35 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 @@
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,43 +75,16 @@
}
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);
}
- template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs>
- solve(const MatrixBase<Rhs>& b) const
- {
- eigen_assert(m_isInitialized && "IncompleteLU is not initialized.");
- eigen_assert(cols()==b.rows()
- && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b");
- return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived());
- }
-
protected:
FactorType m_lu;
- bool m_isInitialized;
};
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<IncompleteLU<_MatrixType>, Rhs>
- : solve_retval_base<IncompleteLU<_MatrixType>, Rhs>
-{
- typedef IncompleteLU<_MatrixType> 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_INCOMPLETE_LU_H
diff --git a/unsupported/Eigen/src/IterativeSolvers/MINRES.h b/unsupported/Eigen/src/IterativeSolvers/MINRES.h
index 30f26aa..256990c 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
@@ -29,7 +29,7 @@
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 @@
}
// 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 @@
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 @@
* 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 @@
* 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 @@
{
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;
@@ -228,46 +229,41 @@
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
- MINRES(const MatrixType& A) : Base(A) {}
+ template<typename MatrixDerived>
+ explicit MINRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
/** 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);
}
@@ -277,33 +273,16 @@
/** \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
diff --git a/unsupported/Eigen/src/IterativeSolvers/Scaling.h b/unsupported/Eigen/src/IterativeSolvers/Scaling.h
index 4fd4392..d113e6e 100644
--- a/unsupported/Eigen/src/IterativeSolvers/Scaling.h
+++ b/unsupported/Eigen/src/IterativeSolvers/Scaling.h
@@ -9,6 +9,9 @@
#ifndef EIGEN_ITERSCALING_H
#define EIGEN_ITERSCALING_H
+
+namespace Eigen {
+
/**
* \ingroup IterativeSolvers_Module
* \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
@@ -41,8 +44,6 @@
*
* \sa \ref IncompleteLUT
*/
-namespace Eigen {
-using std::abs;
template<typename _MatrixType>
class IterScaling
{
@@ -71,6 +72,7 @@
*/
void compute (const MatrixType& mat)
{
+ using std::abs;
int m = mat.rows();
int n = mat.cols();
eigen_assert((m>0 && m == n) && "Please give a non - empty matrix");