diff options
Diffstat (limited to 'unsupported/Eigen/src/IterativeSolvers')
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h | 10 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/DGMRES.h | 122 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/GMRES.h | 38 | ||||
| -rwxr-xr-x | unsupported/Eigen/src/IterativeSolvers/IDRS.h | 436 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/IterationController.h | 2 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/MINRES.h | 38 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/Scaling.h | 6 |
7 files changed, 530 insertions, 122 deletions
diff --git a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h index dc0093eb9..e7d70f39d 100644 --- a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h +++ b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h @@ -31,13 +31,13 @@ #ifndef EIGEN_CONSTRAINEDCG_H #define EIGEN_CONSTRAINEDCG_H -#include <Eigen/Core> +#include "../../../../Eigen/Core" namespace Eigen { namespace internal { -/** \ingroup IterativeSolvers_Module +/** \ingroup IterativeLinearSolvers_Module * Compute the pseudo inverse of the non-square matrix C such that * \f$ CINV = (C * C^T)^{-1} * C \f$ based on a conjugate gradient method. * @@ -96,10 +96,10 @@ void pseudo_inverse(const CMatrix &C, CINVMatrix &CINV) -/** \ingroup IterativeSolvers_Module +/** \ingroup IterativeLinearSolvers_Module * Constrained conjugate gradient * - * Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the contraint \f$ Cx \le f \f$ + * Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the constraint \f$ Cx \le f \f$ */ template<typename TMatrix, typename CMatrix, typename VectorX, typename VectorB, typename VectorF> @@ -158,8 +158,6 @@ void constrained_cg(const TMatrix& A, const CMatrix& C, VectorX& x, rho = r.dot(z); if (iter.finished(rho)) break; - - if (iter.noiseLevel() > 0 && transition) std::cerr << "CCG: transition\n"; if (transition || iter.first()) gamma = 0.0; else gamma = (std::max)(0.0, (rho - old_z.dot(z)) / rho_1); p = z + gamma*p; diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h index bae04fc30..5ae011b75 100644 --- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h @@ -10,7 +10,7 @@ #ifndef EIGEN_DGMRES_H #define EIGEN_DGMRES_H -#include <Eigen/Eigenvalues> +#include "../../../../Eigen/Eigenvalues" namespace Eigen { @@ -39,7 +39,6 @@ template <typename VectorType, typename IndexType> void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut) { eigen_assert(vec.size() == perm.size()); - typedef typename IndexType::Scalar Index; bool flag; for (Index k = 0; k < ncut; k++) { @@ -58,7 +57,7 @@ void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType:: } /** - * \ingroup IterativeLInearSolvers_Module + * \ingroup IterativeLinearSolvers_Module * \brief A Restarted GMRES with deflation. * This class implements a modification of the GMRES solver for * sparse linear systems. The basis is built with modified @@ -89,7 +88,7 @@ void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType:: * [1] D. NUENTSA WAKAM and F. PACULL, Memory Efficient Hybrid * Algebraic Solvers for Linear Systems Arising from Compressible * Flows, Computers and Fluids, In Press, - * http://dx.doi.org/10.1016/j.compfluid.2012.03.023 + * https://doi.org/10.1016/j.compfluid.2012.03.023 * [2] K. Burrage and J. Erhel, On the performance of various * adaptive preconditioned GMRES strategies, 5(1998), 101-121. * [3] J. Erhel, K. Burrage and B. Pohl, Restarted GMRES @@ -110,9 +109,9 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > using Base::m_tolerance; public: using Base::_solve_impl; + using Base::_solve_with_guess_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; @@ -143,44 +142,30 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > /** \internal */ template<typename Rhs,typename Dest> - void _solve_with_guess_impl(const Rhs& b, Dest& x) const - { - bool failed = false; - for(int j=0; j<b.cols(); ++j) - { - m_iterations = Base::maxIterations(); - m_error = Base::m_tolerance; - - typename Dest::ColXpr xj(x,j); - dgmres(matrix(), b.col(j), xj, Base::m_preconditioner); - } - m_info = failed ? NumericalIssue - : m_error <= Base::m_tolerance ? Success - : NoConvergence; - m_isInitialized = true; - } - - /** \internal */ - template<typename Rhs,typename Dest> - void _solve_impl(const Rhs& b, MatrixBase<Dest>& x) const + void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const { - x = b; - _solve_with_guess_impl(b,x.derived()); + EIGEN_STATIC_ASSERT(Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1, YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX); + + m_iterations = Base::maxIterations(); + m_error = Base::m_tolerance; + + dgmres(matrix(), b, x, Base::m_preconditioner); } + /** * 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 @@ -189,12 +174,12 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > /** * 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 @@ -202,27 +187,27 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > 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, StorageIndex& 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 void dgmresInitDeflation(Index& rows) const; mutable DenseMatrix m_V; // Krylov basis vectors mutable DenseMatrix m_H; // Hessenberg matrix - mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied + mutable DenseMatrix m_Hes; // Initial hessenberg matrix without Givens rotations applied mutable Index m_restart; // Maximum size of the Krylov subspace mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace 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 StorageIndex 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 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; @@ -243,18 +228,30 @@ template<typename Rhs, typename Dest> void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const { + const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); + + RealScalar normRhs = rhs.norm(); + if(normRhs <= considerAsZero) + { + x.setZero(); + m_error = 0; + return; + } + //Initialization - int n = mat.rows(); + m_isDeflInitialized = false; + 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); - //Initial residual vector and intial norm - x = precond.solve(x); + //Initial residual vector and initial norm + if(x.squaredNorm()==0) + x = precond.solve(rhs); r0 = rhs - mat * x; RealScalar beta = r0.norm(); - RealScalar normRhs = rhs.norm(); + m_error = beta/normRhs; if(m_error < m_tolerance) m_info = Success; @@ -267,8 +264,10 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts); // Compute the new residual vector for the restart - if (nbIts < m_iterations && m_info == NoConvergence) - r0 = rhs - mat * x; + if (nbIts < m_iterations && m_info == NoConvergence) { + r0 = rhs - mat * x; + beta = r0.norm(); + } } } @@ -284,7 +283,7 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh */ 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 @@ -293,8 +292,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con 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) { @@ -312,7 +311,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con // 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); @@ -328,7 +327,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con // 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()); } @@ -394,7 +393,6 @@ inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_Matr 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); @@ -418,7 +416,7 @@ inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_Matr } template< typename _MatrixType, typename _Preconditioner> -int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& 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; @@ -433,8 +431,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri // 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) @@ -442,7 +440,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri 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++; @@ -451,7 +449,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri // 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)); } @@ -462,8 +460,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri 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); } @@ -473,7 +471,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri 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); @@ -488,8 +486,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri } // 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; @@ -502,7 +500,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri } 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); diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h index 5a82b0df6..ff912094f 100644 --- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -21,7 +21,7 @@ namespace internal { * * Parameters: * \param mat matrix of linear system of equations -* \param Rhs right hand side vector 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 @@ -64,6 +64,15 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition typedef Matrix < Scalar, Dynamic, 1 > VectorType; typedef Matrix < Scalar, Dynamic, Dynamic, ColMajor> FMatrixType; + const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); + + if(rhs.norm() <= considerAsZero) + { + x.setZero(); + tol_error = 0; + return true; + } + RealScalar tol = tol_error; const Index maxIters = iters; iters = 0; @@ -307,31 +316,14 @@ public: /** \internal */ template<typename Rhs,typename Dest> - void _solve_with_guess_impl(const Rhs& b, Dest& x) const + void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const { - bool failed = false; - 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(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error)) - failed = true; - } - m_info = failed ? NumericalIssue + m_iterations = Base::maxIterations(); + m_error = Base::m_tolerance; + bool ret = internal::gmres(matrix(), b, x, Base::m_preconditioner, m_iterations, m_restart, m_error); + m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence; - m_isInitialized = true; - } - - /** \internal */ - template<typename Rhs,typename Dest> - void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const - { - x = b; - if(x.squaredNorm() == 0) return; // Check Zero right hand side - _solve_with_guess_impl(b,x.derived()); } protected: diff --git a/unsupported/Eigen/src/IterativeSolvers/IDRS.h b/unsupported/Eigen/src/IterativeSolvers/IDRS.h new file mode 100755 index 000000000..90d20fad4 --- /dev/null +++ b/unsupported/Eigen/src/IterativeSolvers/IDRS.h @@ -0,0 +1,436 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2020 Chris Schoutrop <c.e.m.schoutrop@tue.nl> +// Copyright (C) 2020 Jens Wehner <j.wehner@esciencecenter.nl> +// Copyright (C) 2020 Jan van Dijk <j.v.dijk@tue.nl> +// +// 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_IDRS_H +#define EIGEN_IDRS_H + +namespace Eigen +{ + + namespace internal + { + /** \internal Low-level Induced Dimension Reduction algoritm + \param A The matrix A + \param b The right hand side vector b + \param x On input and initial solution, on output the computed solution. + \param precond A preconditioner being able to efficiently solve for an + approximation of Ax=b (regardless of b) + \param iter On input the max number of iteration, on output the number of performed iterations. + \param relres On input the tolerance error, on output an estimation of the relative error. + \param S On input Number of the dimension of the shadow space. + \param smoothing switches residual smoothing on. + \param angle small omega lead to faster convergence at the expense of numerical stability + \param replacement switches on a residual replacement strategy to increase accuracy of residual at the expense of more Mat*vec products + \return false in the case of numerical issue, for example a break down of IDRS. + */ + template<typename Vector, typename RealScalar> + typename Vector::Scalar omega(const Vector& t, const Vector& s, RealScalar angle) + { + using numext::abs; + typedef typename Vector::Scalar Scalar; + const RealScalar ns = s.norm(); + const RealScalar nt = t.norm(); + const Scalar ts = t.dot(s); + const RealScalar rho = abs(ts / (nt * ns)); + + if (rho < angle) { + if (ts == Scalar(0)) { + return Scalar(0); + } + // Original relation for om is given by + // om = om * angle / rho; + // To alleviate potential (near) division by zero this can be rewritten as + // om = angle * (ns / nt) * (ts / abs(ts)) = angle * (ns / nt) * sgn(ts) + return angle * (ns / nt) * (ts / abs(ts)); + } + return ts / (nt * nt); + } + + template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> + bool idrs(const MatrixType& A, const Rhs& b, Dest& x, const Preconditioner& precond, + Index& iter, + typename Dest::RealScalar& relres, Index S, bool smoothing, typename Dest::RealScalar angle, bool replacement) + { + typedef typename Dest::RealScalar RealScalar; + typedef typename Dest::Scalar Scalar; + typedef Matrix<Scalar, Dynamic, 1> VectorType; + typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> DenseMatrixType; + const Index N = b.size(); + S = S < x.rows() ? S : x.rows(); + const RealScalar tol = relres; + const Index maxit = iter; + + Index replacements = 0; + bool trueres = false; + + FullPivLU<DenseMatrixType> lu_solver; + + DenseMatrixType P; + { + HouseholderQR<DenseMatrixType> qr(DenseMatrixType::Random(N, S)); + P = (qr.householderQ() * DenseMatrixType::Identity(N, S)); + } + + const RealScalar normb = b.norm(); + + if (internal::isApprox(normb, RealScalar(0))) + { + //Solution is the zero vector + x.setZero(); + iter = 0; + relres = 0; + return true; + } + // from http://homepage.tudelft.nl/1w5b5/IDRS/manual.pdf + // A peak in the residual is considered dangerously high if‖ri‖/‖b‖> C(tol/epsilon). + // With epsilon the + // relative machine precision. The factor tol/epsilon corresponds to the size of a + // finite precision number that is so large that the absolute round-off error in + // this number, when propagated through the process, makes it impossible to + // achieve the required accuracy.The factor C accounts for the accumulation of + // round-off errors. This parameter has beenset to 10−3. + // mp is epsilon/C + // 10^3 * eps is very conservative, so normally no residual replacements will take place. + // It only happens if things go very wrong. Too many restarts may ruin the convergence. + const RealScalar mp = RealScalar(1e3) * NumTraits<Scalar>::epsilon(); + + + + //Compute initial residual + const RealScalar tolb = tol * normb; //Relative tolerance + VectorType r = b - A * x; + + VectorType x_s, r_s; + + if (smoothing) + { + x_s = x; + r_s = r; + } + + RealScalar normr = r.norm(); + + if (normr <= tolb) + { + //Initial guess is a good enough solution + iter = 0; + relres = normr / normb; + return true; + } + + DenseMatrixType G = DenseMatrixType::Zero(N, S); + DenseMatrixType U = DenseMatrixType::Zero(N, S); + DenseMatrixType M = DenseMatrixType::Identity(S, S); + VectorType t(N), v(N); + Scalar om = 1.; + + //Main iteration loop, guild G-spaces: + iter = 0; + + while (normr > tolb && iter < maxit) + { + //New right hand size for small system: + VectorType f = (r.adjoint() * P).adjoint(); + + for (Index k = 0; k < S; ++k) + { + //Solve small system and make v orthogonal to P: + //c = M(k:s,k:s)\f(k:s); + lu_solver.compute(M.block(k , k , S -k, S - k )); + VectorType c = lu_solver.solve(f.segment(k , S - k )); + //v = r - G(:,k:s)*c; + v = r - G.rightCols(S - k ) * c; + //Preconditioning + v = precond.solve(v); + + //Compute new U(:,k) and G(:,k), G(:,k) is in space G_j + U.col(k) = U.rightCols(S - k ) * c + om * v; + G.col(k) = A * U.col(k ); + + //Bi-Orthogonalise the new basis vectors: + for (Index i = 0; i < k-1 ; ++i) + { + //alpha = ( P(:,i)'*G(:,k) )/M(i,i); + Scalar alpha = P.col(i ).dot(G.col(k )) / M(i, i ); + G.col(k ) = G.col(k ) - alpha * G.col(i ); + U.col(k ) = U.col(k ) - alpha * U.col(i ); + } + + //New column of M = P'*G (first k-1 entries are zero) + //M(k:s,k) = (G(:,k)'*P(:,k:s))'; + M.block(k , k , S - k , 1) = (G.col(k ).adjoint() * P.rightCols(S - k )).adjoint(); + + if (internal::isApprox(M(k,k), Scalar(0))) + { + return false; + } + + //Make r orthogonal to q_i, i = 0..k-1 + Scalar beta = f(k ) / M(k , k ); + r = r - beta * G.col(k ); + x = x + beta * U.col(k ); + normr = r.norm(); + + if (replacement && normr > tolb / mp) + { + trueres = true; + } + + //Smoothing: + if (smoothing) + { + t = r_s - r; + //gamma is a Scalar, but the conversion is not allowed + Scalar gamma = t.dot(r_s) / t.norm(); + r_s = r_s - gamma * t; + x_s = x_s - gamma * (x_s - x); + normr = r_s.norm(); + } + + if (normr < tolb || iter == maxit) + { + break; + } + + //New f = P'*r (first k components are zero) + if (k < S-1) + { + f.segment(k + 1, S - (k + 1) ) = f.segment(k + 1 , S - (k + 1)) - beta * M.block(k + 1 , k , S - (k + 1), 1); + } + }//end for + + if (normr < tolb || iter == maxit) + { + break; + } + + //Now we have sufficient vectors in G_j to compute residual in G_j+1 + //Note: r is already perpendicular to P so v = r + //Preconditioning + v = r; + v = precond.solve(v); + + //Matrix-vector multiplication: + t = A * v; + + //Computation of a new omega + om = internal::omega(t, r, angle); + + if (om == RealScalar(0.0)) + { + return false; + } + + r = r - om * t; + x = x + om * v; + normr = r.norm(); + + if (replacement && normr > tolb / mp) + { + trueres = true; + } + + //Residual replacement? + if (trueres && normr < normb) + { + r = b - A * x; + trueres = false; + replacements++; + } + + //Smoothing: + if (smoothing) + { + t = r_s - r; + Scalar gamma = t.dot(r_s) /t.norm(); + r_s = r_s - gamma * t; + x_s = x_s - gamma * (x_s - x); + normr = r_s.norm(); + } + + iter++; + + }//end while + + if (smoothing) + { + x = x_s; + } + relres=normr/normb; + return true; + } + + } // namespace internal + + template <typename _MatrixType, typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > + class IDRS; + + namespace internal + { + + template <typename _MatrixType, typename _Preconditioner> + struct traits<Eigen::IDRS<_MatrixType, _Preconditioner> > + { + typedef _MatrixType MatrixType; + typedef _Preconditioner Preconditioner; + }; + + } // namespace internal + + +/** \ingroup IterativeLinearSolvers_Module + * \brief The Induced Dimension Reduction method (IDR(s)) is a short-recurrences Krylov method for sparse square problems. + * + * This class allows to solve for A.x = b sparse linear problems. The vectors x and b can be either dense or sparse. + * he Induced Dimension Reduction method, IDR(), is a robust and efficient short-recurrence Krylov subspace method for + * solving large nonsymmetric systems of linear equations. + * + * For indefinite systems IDR(S) outperforms both BiCGStab and BiCGStab(L). Additionally, IDR(S) can handle matrices + * with complex eigenvalues more efficiently than BiCGStab. + * + * Many problems that do not converge for BiCGSTAB converge for IDR(s) (for larger values of s). And if both methods + * converge the convergence for IDR(s) is typically much faster for difficult systems (for example indefinite problems). + * + * IDR(s) is a limited memory finite termination method. In exact arithmetic it converges in at most N+N/s iterations, + * with N the system size. It uses a fixed number of 4+3s vector. In comparison, BiCGSTAB terminates in 2N iterations + * and uses 7 vectors. GMRES terminates in at most N iterations, and uses I+3 vectors, with I the number of iterations. + * Restarting GMRES limits the memory consumption, but destroys the finite termination property. + * + * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix. + * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner + * + * \implsparsesolverconcept + * + * 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. + * + * The tolerance corresponds to the relative residual error: |Ax-b|/|b| + * + * \b Performance: when using sparse matrices, best performance is achied for a row-major sparse matrix format. + * Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled. + * See \ref TopicMultiThreading for details. + * + * By default the iterations start with x=0 as an initial guess of the solution. + * One can control the start using the solveWithGuess() method. + * + * IDR(s) 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 IDRS : public IterativeSolverBase<IDRS<_MatrixType, _Preconditioner> > + { + + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef _Preconditioner Preconditioner; + + private: + typedef IterativeSolverBase<IDRS> Base; + using Base::m_error; + using Base::m_info; + using Base::m_isInitialized; + using Base::m_iterations; + using Base::matrix; + Index m_S; + bool m_smoothing; + RealScalar m_angle; + bool m_residual; + + public: + /** Default constructor. */ + IDRS(): m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {} + + /** 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. + */ + template <typename MatrixDerived> + explicit IDRS(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_S(4), m_smoothing(false), + m_angle(RealScalar(0.7)), m_residual(false) {} + + + /** \internal */ + /** Loops over the number of columns of b and does the following: + 1. sets the tolerence and maxIterations + 2. Calls the function that has the core solver routine + */ + template <typename Rhs, typename Dest> + void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const + { + m_iterations = Base::maxIterations(); + m_error = Base::m_tolerance; + + bool ret = internal::idrs(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, m_S,m_smoothing,m_angle,m_residual); + + m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence; + } + + /** Sets the parameter S, indicating the dimension of the shadow space. Default is 4*/ + void setS(Index S) + { + if (S < 1) + { + S = 4; + } + + m_S = S; + } + + /** Switches off and on smoothing. + Residual smoothing results in monotonically decreasing residual norms at + the expense of two extra vectors of storage and a few extra vector + operations. Although monotonic decrease of the residual norms is a + desirable property, the rate of convergence of the unsmoothed process and + the smoothed process is basically the same. Default is off */ + void setSmoothing(bool smoothing) + { + m_smoothing=smoothing; + } + + /** The angle must be a real scalar. In IDR(s), a value for the + iteration parameter omega must be chosen in every s+1th step. The most + natural choice is to select a value to minimize the norm of the next residual. + This corresponds to the parameter omega = 0. In practice, this may lead to + values of omega that are so small that the other iteration parameters + cannot be computed with sufficient accuracy. In such cases it is better to + increase the value of omega sufficiently such that a compromise is reached + between accurate computations and reduction of the residual norm. The + parameter angle =0.7 (”maintaining the convergence strategy”) + results in such a compromise. */ + void setAngle(RealScalar angle) + { + m_angle=angle; + } + + /** The parameter replace is a logical that determines whether a + residual replacement strategy is employed to increase the accuracy of the + solution. */ + void setResidualUpdate(bool update) + { + m_residual=update; + } + + }; + +} // namespace Eigen + +#endif /* EIGEN_IDRS_H */ diff --git a/unsupported/Eigen/src/IterativeSolvers/IterationController.h b/unsupported/Eigen/src/IterativeSolvers/IterationController.h index c9c1a4be2..a116e09e2 100644 --- a/unsupported/Eigen/src/IterativeSolvers/IterationController.h +++ b/unsupported/Eigen/src/IterativeSolvers/IterationController.h @@ -60,7 +60,7 @@ namespace Eigen { -/** \ingroup IterativeSolvers_Module +/** \ingroup IterativeLinearSolvers_Module * \class IterationController * * \brief Controls the iterations of the iterative solvers diff --git a/unsupported/Eigen/src/IterativeSolvers/MINRES.h b/unsupported/Eigen/src/IterativeSolvers/MINRES.h index 256990c1a..5db454d24 100644 --- a/unsupported/Eigen/src/IterativeSolvers/MINRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/MINRES.h @@ -3,6 +3,7 @@ // // Copyright (C) 2012 Giacomo Po <gpo@ucla.edu> // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2018 David Hyde <dabh@stanford.edu> // // 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 @@ -64,8 +65,6 @@ namespace Eigen { eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE"); RealScalar beta_new(sqrt(beta_new2)); const RealScalar beta_one(beta_new); - v_new /= beta_new; - w_new /= beta_new; // Initialize other variables RealScalar c(1.0); // the cosine of the Givens rotation RealScalar c_old(1.0); @@ -83,18 +82,18 @@ namespace Eigen { /* Note that there are 4 variants on the Lanczos algorithm. These are * described in Paige, C. C. (1972). Computational variants of * the Lanczos method for the eigenproblem. IMA Journal of Applied - * Mathematics, 10(3), 373–381. The current implementation corresponds + * Mathematics, 10(3), 373-381. The current implementation corresponds * to the case A(2,7) in the paper. It also corresponds to - * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear + * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear * Systems, 2003 p.173. For the preconditioned version see * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987). */ const RealScalar beta(beta_new); v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter -// const VectorType v_old(v); // NOT SURE IF CREATING v_old EVERY ITERATION IS EFFICIENT + v_new /= beta_new; // overwrite v_new for next iteration + w_new /= beta_new; // overwrite w_new for next iteration v = v_new; // update w = w_new; // update -// const VectorType w(w_new); // NOT SURE IF CREATING w EVERY ITERATION IS EFFICIENT v_new.noalias() = mat*w - beta*v_old; // compute v_new const RealScalar alpha = v_new.dot(w); v_new -= alpha*v; // overwrite v_new @@ -102,8 +101,6 @@ namespace Eigen { beta_new2 = v_new.dot(w_new); // compute beta_new eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE"); beta_new = sqrt(beta_new2); // compute beta_new - v_new /= beta_new; // overwrite v_new for next iteration - w_new /= beta_new; // overwrite w_new for next iteration // Givens rotation const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration @@ -117,7 +114,6 @@ namespace Eigen { // Update solution p_oold = p_old; -// const VectorType p_oold(p_old); // NOT SURE IF CREATING p_oold EVERY ITERATION IS EFFICIENT p_old = p; p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED? x += beta_one*c*eta*p; @@ -237,7 +233,7 @@ namespace Eigen { /** \internal */ template<typename Rhs,typename Dest> - void _solve_with_guess_impl(const Rhs& b, Dest& x) const + void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const { typedef typename Base::MatrixWrapper MatrixWrapper; typedef typename Base::ActualMatrixType ActualMatrixType; @@ -257,28 +253,11 @@ namespace Eigen { 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(SelfAdjointWrapper(row_mat), b.col(j), xj, - Base::m_preconditioner, m_iterations, m_error); - } - - m_isInitialized = true; + internal::minres(SelfAdjointWrapper(row_mat), b, x, + Base::m_preconditioner, m_iterations, m_error); m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; } - /** \internal */ - template<typename Rhs,typename Dest> - void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const - { - x.setZero(); - _solve_with_guess_impl(b,x.derived()); - } - protected: }; @@ -286,4 +265,3 @@ namespace Eigen { } // end namespace Eigen #endif // EIGEN_MINRES_H - diff --git a/unsupported/Eigen/src/IterativeSolvers/Scaling.h b/unsupported/Eigen/src/IterativeSolvers/Scaling.h index d113e6e90..9b3eb53e0 100644 --- a/unsupported/Eigen/src/IterativeSolvers/Scaling.h +++ b/unsupported/Eigen/src/IterativeSolvers/Scaling.h @@ -104,12 +104,18 @@ class IterScaling for (int i = 0; i < m; ++i) { Dr(i) = std::sqrt(Dr(i)); + } + for (int i = 0; i < n; ++i) + { Dc(i) = std::sqrt(Dc(i)); } // Save the scaling factors for (int i = 0; i < m; ++i) { m_left(i) /= Dr(i); + } + for (int i = 0; i < n; ++i) + { m_right(i) /= Dc(i); } // Scale the rows and the columns of the matrix |