// 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/.
/** \brief Computes a permutation vector to have a sorted sequence * \param vec The vector to reorder. * \param perm gives the sorted sequence on output. Must be initialized with 0..n-1 * \param ncut Put the ncut smallest elements at the end of the vector * WARNING This is an expensive sort, so should be used only * for small size vectors * TODO Use modified QuickSplit or std::nth_element to get the smallest values
*/ template <typename VectorType, typename IndexType> void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut)
{
eigen_assert(vec.size() == perm.size()); bool flag; for (Index k = 0; k < ncut; k++)
{
flag = false; for (Index j = 0; j < vec.size()-1; j++)
{ if ( vec(perm(j)) < vec(perm(j+1)) )
{
std::swap(perm(j),perm(j+1));
flag = true;
} if (!flag) break; // The vector is in sorted order
}
}
}
} /** * \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 * Gram-Schmidt. At each restart, a few approximated eigenvectors * corresponding to the smallest eigenvalues are used to build a * preconditioner for the next cycle. This preconditioner * for deflation can be combined with any other preconditioner, * the IncompleteLUT for instance. The preconditioner is applied * at right of the matrix and the combination is multiplicative. * * \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 * Typical usage : * \code * SparseMatrix<double> A; * VectorXd x, b; * //Fill A and b ... * DGMRES<SparseMatrix<double> > solver; * solver.set_restart(30); // Set restarting value * solver.setEigenv(1); // Set the number of eigenvalues to deflate * solver.compute(A); * 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 * Flows, Computers and Fluids, In Press, * 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 * preconditioned by deflation,J. Computational and Applied * Mathematics, 69(1996), 303-318.
*
*/ template< typename _MatrixType, typename _Preconditioner> class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
{ typedef IterativeSolverBase<DGMRES> Base; 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; using Base::_solve_with_guess_impl; typedef _MatrixType MatrixType; typedeftypename MatrixType::Scalar Scalar; typedeftypename MatrixType::StorageIndex StorageIndex; typedeftypename MatrixType::RealScalar RealScalar; typedef _Preconditioner Preconditioner; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; typedef Matrix<RealScalar,Dynamic,1> DenseRealVector; typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector;
/** 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 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) {}
/** * Get the restart value
*/
Index restart() { return m_restart; }
/** * Set the restart value (default is 30)
*/ void set_restart(const Index restart) { m_restart=restart; }
/** * Set the number of eigenvalues to deflate at each restart
*/ void setEigenv(const Index neig)
{
m_neig = neig; if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
}
/** * Get the size of the deflation subspace size
*/
Index deflSize() {return m_r; }
/** * Set the maximum size of the deflation subspace
*/ void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; }
protected: // DGMRES algorithm template<typename Rhs, typename Dest> void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const; // Perform one cycle of GMRES template<typename Dest>
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
Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const; // Apply deflation to a vector template<typename RhsType, typename DestType>
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 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 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 mutablebool m_isDeflAllocated; mutablebool m_isDeflInitialized;
//Adaptive strategy mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed mutablebool m_force; // Force the use of deflation at each restart
}; /** * \brief Perform several cycles of restarted GMRES with modified Gram Schmidt, * * A right preconditioner is used combined with deflation. *
*/ template< typename _MatrixType, typename _Preconditioner> 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)();
// Iterative process while (nbIts < m_iterations && m_info == NoConvergence)
{
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;
beta = r0.norm();
}
}
}
/** * \brief Perform one restart cycle of DGMRES * \param mat The coefficient matrix * \param precond The preconditioner * \param x the new approximated solution * \param r0 The initial residual vector * \param beta The norm of the residual computed so far * \param normRhs The norm of the right hand side vector * \param nbIts The number of iterations
*/ template< typename _MatrixType, typename _Preconditioner> template<typename Dest>
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
g.setZero();
g(0) = Scalar(beta);
m_V.col(0) = r0/beta;
m_info = NoConvergence;
std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
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)
{ // Apply preconditioner(s) at right if (m_isDeflInitialized )
{
dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
tv2 = precond.solve(tv1);
} else
{
tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
}
tv1 = mat * tv2;
// Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
Scalar coef; for (Index i = 0; i <= it; ++i)
{
coef = tv1.dot(m_V.col(i));
tv1 = tv1 - coef * m_V.col(i);
m_H(i,it) = coef;
m_Hes(i,it) = coef;
} // Normalize the vector
coef = tv1.norm();
m_V.col(it+1) = tv1/coef;
m_H(it+1, it) = coef; // m_Hes(it+1,it) = coef;
// FIXME Check for happy breakdown
// Update Hessenberg matrix with Givens rotations for (Index i = 1; i <= it; ++i)
{
m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
} // Compute the new plane rotation
gr[it].makeGivens(m_H(it, it), m_H(it+1,it)); // Apply the new rotation
m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
g.applyOnTheLeft(it,it+1, gr[it].adjoint());
if (m_error < m_tolerance)
{ // The method has converged
m_info = Success; break;
}
}
// Compute the new coefficients by solving the least square problem // it++; //FIXME Check first if the matrix is singular ... zero diagonal
DenseVector nrs(m_restart);
nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it));
// Form the new solution if (m_isDeflInitialized)
{
tv1 = m_V.leftCols(it) * nrs;
dgmresApplyDeflation(tv1, tv2);
x = x + precond.solve(tv2);
} else
x = x + precond.solve(m_V.leftCols(it) * nrs);
// Go for a new cycle and compute data for deflation if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
dgmresComputeDeflationData(mat, precond, it, m_neig); return 0;
// Reorder the absolute values of Schur values
DenseRealVector modulEig(it); 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)
{
m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
} //Count the real number of extracted eigenvalues (with complex conjugates)
Index nbrEig = 0; while (nbrEig < neig)
{ if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; else nbrEig += 2;
} // Extract the Schur vectors corresponding to the smallest Ritz values
DenseMatrix Sr(it, nbrEig);
Sr.setZero(); for (Index j = 0; j < nbrEig; j++)
{
Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
}
// Form the Schur vectors of the initial matrix using the Krylov basis
DenseMatrix X;
X = m_V.leftCols(it) * Sr; if (m_r)
{ // Orthogonalize X against m_U using modified Gram-Schmidt 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);
}
// Compute m_MX = A * M^-1 * X
Index m = m_V.rows(); if (!m_isDeflAllocated)
dgmresInitDeflation(m);
DenseMatrix MX(m, nbrEig);
DenseVector tv1(m); for (Index j = 0; j < nbrEig; j++)
{
tv1 = mat * X.col(j);
MX.col(j) = precond.solve(tv1);
}
// Save X into m_U and m_MX in m_MU 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;
// Factorize m_T into m_luT
m_luT.compute(m_T.topLeftCorner(m_r, m_r));
//FIXME CHeck if the factorization was correctly done (nonsingular matrix)
m_isDeflInitialized = true; return 0;
} template<typename _MatrixType, typename _Preconditioner> template<typename RhsType, typename DestType>
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;
}
} // end namespace Eigen #endif
¤ Dauer der Verarbeitung: 0.16 Sekunden
(vorverarbeitet)
¤
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
