GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

Idrs refactoring

Browse Source
This commit is contained in:
Jens Wehner 2021-12-02 23:32:07 +00:00 committed by Rasmus Munk Larsen
parent f63c6dd1f9
commit 4ee2e9b340
1 changed files with 323 additions and 364 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

141
unsupported/Eigen/src/IterativeSolvers/IDRS.h
View File

@ -9,17 +9,14 @@
// Public License v. 2.0. If a copy of the MPL was not distributed // 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/. // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_IDRS_H #ifndef EIGEN_IDRS_H
#define EIGEN_IDRS_H #define EIGEN_IDRS_H
#include "./InternalHeaderCheck.h" #include "./InternalHeaderCheck.h"
namespace Eigen namespace Eigen {
{
namespace internal namespace internal {
{
/** \internal Low-level Induced Dimension Reduction algorithm /** \internal Low-level Induced Dimension Reduction algorithm
\param A The matrix A \param A The matrix A
\param b The right hand side vector b \param b The right hand side vector b
@ -31,16 +28,15 @@ namespace Eigen
\param S On input Number of the dimension of the shadow space. \param S On input Number of the dimension of the shadow space.
\param smoothing switches residual smoothing on. \param smoothing switches residual smoothing on.
\param angle small omega lead to faster convergence at the expense of numerical stability \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 \param replacement switches on a residual replacement strategy to increase accuracy of residual at the
\return false in the case of numerical issue, for example a break down of IDRS. 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> template <typename Vector, typename RealScalar>
typename Vector::Scalar omega(const Vector& t, const Vector& s, RealScalar angle) typename Vector::Scalar omega(const Vector& t, const Vector& s, RealScalar angle) {
{
using numext::abs; using numext::abs;
typedef typename Vector::Scalar Scalar; typedef typename Vector::Scalar Scalar;
const RealScalar ns = s.norm(); const RealScalar ns = s.StableNorm();
const RealScalar nt = t.norm(); const RealScalar nt = t.StableNorm();
const Scalar ts = t.dot(s); const Scalar ts = t.dot(s);
const RealScalar rho = abs(ts / (nt * ns)); const RealScalar rho = abs(ts / (nt * ns));
@ -58,10 +54,9 @@ namespace Eigen
} }
template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
bool idrs(const MatrixType& A, const Rhs& b, Dest& x, const Preconditioner& precond, bool idrs(const MatrixType& A, const Rhs& b, Dest& x, const Preconditioner& precond, Index& iter,
Index& iter, typename Dest::RealScalar& relres, Index S, bool smoothing, typename Dest::RealScalar angle,
typename Dest::RealScalar& relres, Index S, bool smoothing, typename Dest::RealScalar angle, bool replacement) bool replacement) {
{
typedef typename Dest::RealScalar RealScalar; typedef typename Dest::RealScalar RealScalar;
typedef typename Dest::Scalar Scalar; typedef typename Dest::Scalar Scalar;
typedef Matrix<Scalar, Dynamic, 1> VectorType; typedef Matrix<Scalar, Dynamic, 1> VectorType;
@ -82,10 +77,9 @@ namespace Eigen
P = (qr.householderQ() * DenseMatrixType::Identity(N, S)); P = (qr.householderQ() * DenseMatrixType::Identity(N, S));
} }
const RealScalar normb = b.norm(); const RealScalar normb = b.StableNorm();
if (internal::isApprox(normb, RealScalar(0))) if (internal::isApprox(normb, RealScalar(0))) {
{
// Solution is the zero vector // Solution is the zero vector
x.setZero(); x.setZero();
iter = 0; iter = 0;
@ -105,24 +99,20 @@ namespace Eigen
// It only happens if things go very wrong. Too many restarts may ruin the convergence. // It only happens if things go very wrong. Too many restarts may ruin the convergence.
const RealScalar mp = RealScalar(1e3) * NumTraits<Scalar>::epsilon(); const RealScalar mp = RealScalar(1e3) * NumTraits<Scalar>::epsilon();
// Compute initial residual // Compute initial residual
const RealScalar tolb = tol * normb; // Relative tolerance const RealScalar tolb = tol * normb; // Relative tolerance
VectorType r = b - A * x; VectorType r = b - A * x;
VectorType x_s, r_s; VectorType x_s, r_s;
if (smoothing) if (smoothing) {
{
x_s = x; x_s = x;
r_s = r; r_s = r;
} }
RealScalar normr = r.norm(); RealScalar normr = r.StableNorm();
if (normr <= tolb) if (normr <= tolb) {
{
// Initial guess is a good enough solution // Initial guess is a good enough solution
iter = 0; iter = 0;
relres = normr / normb; relres = normr / normb;
@ -138,13 +128,11 @@ namespace Eigen
// Main iteration loop, guild G-spaces: // Main iteration loop, guild G-spaces:
iter = 0; iter = 0;
while (normr > tolb && iter < maxit) while (normr > tolb && iter < maxit) {
{
// New right hand size for small system: // New right hand size for small system:
VectorType f = (r.adjoint() * P).adjoint(); VectorType f = (r.adjoint() * P).adjoint();
for (Index k = 0; k < S; ++k) for (Index k = 0; k < S; ++k) {
{
// Solve small system and make v orthogonal to P: // Solve small system and make v orthogonal to P:
// c = M(k:s,k:s)\f(k:s); // c = M(k:s,k:s)\f(k:s);
lu_solver.compute(M.block(k, k, S - k, S - k)); lu_solver.compute(M.block(k, k, S - k, S - k));
@ -159,8 +147,7 @@ namespace Eigen
G.col(k) = A * U.col(k); G.col(k) = A * U.col(k);
// Bi-Orthogonalise the new basis vectors: // Bi-Orthogonalise the new basis vectors:
for (Index i = 0; i < k-1 ; ++i) for (Index i = 0; i < k - 1; ++i) {
{
// alpha = ( P(:,i)'*G(:,k) )/M(i,i); // alpha = ( P(:,i)'*G(:,k) )/M(i,i);
Scalar alpha = P.col(i).dot(G.col(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); G.col(k) = G.col(k) - alpha * G.col(i);
@ -171,8 +158,7 @@ namespace Eigen
// M(k:s,k) = (G(:,k)'*P(:,k:s))'; // 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(); M.block(k, k, S - k, 1) = (G.col(k).adjoint() * P.rightCols(S - k)).adjoint();
if (internal::isApprox(M(k,k), Scalar(0))) if (internal::isApprox(M(k, k), Scalar(0))) {
{
return false; return false;
} }
@ -180,38 +166,33 @@ namespace Eigen
Scalar beta = f(k) / M(k, k); Scalar beta = f(k) / M(k, k);
r = r - beta * G.col(k); r = r - beta * G.col(k);
x = x + beta * U.col(k); x = x + beta * U.col(k);
normr = r.norm(); normr = r.StableNorm();
if (replacement && normr > tolb / mp) if (replacement && normr > tolb / mp) {
{
trueres = true; trueres = true;
} }
// Smoothing: // Smoothing:
if (smoothing) if (smoothing) {
{
t = r_s - r; t = r_s - r;
// gamma is a Scalar, but the conversion is not allowed // gamma is a Scalar, but the conversion is not allowed
Scalar gamma = t.dot(r_s) / t.norm(); Scalar gamma = t.dot(r_s) / t.StableNorm();
r_s = r_s - gamma * t; r_s = r_s - gamma * t;
x_s = x_s - gamma * (x_s - x); x_s = x_s - gamma * (x_s - x);
normr = r_s.norm(); normr = r_s.StableNorm();
} }
if (normr < tolb || iter == maxit) if (normr < tolb || iter == maxit) {
{
break; break;
} }
// New f = P'*r (first k components are zero) // New f = P'*r (first k components are zero)
if (k < S-1) 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); 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 } // end for
if (normr < tolb || iter == maxit) if (normr < tolb || iter == maxit) {
{
break; break;
} }
@ -227,44 +208,39 @@ namespace Eigen
// Computation of a new omega // Computation of a new omega
om = internal::omega(t, r, angle); om = internal::omega(t, r, angle);
if (om == RealScalar(0.0)) if (om == RealScalar(0.0)) {
{
return false; return false;
} }
r = r - om * t; r = r - om * t;
x = x + om * v; x = x + om * v;
normr = r.norm(); normr = r.StableNorm();
if (replacement && normr > tolb / mp) if (replacement && normr > tolb / mp) {
{
trueres = true; trueres = true;
} }
// Residual replacement? // Residual replacement?
if (trueres && normr < normb) if (trueres && normr < normb) {
{
r = b - A * x; r = b - A * x;
trueres = false; trueres = false;
replacements++; replacements++;
} }
// Smoothing: // Smoothing:
if (smoothing) if (smoothing) {
{
t = r_s - r; t = r_s - r;
Scalar gamma = t.dot(r_s) /t.norm(); Scalar gamma = t.dot(r_s) / t.StableNorm();
r_s = r_s - gamma * t; r_s = r_s - gamma * t;
x_s = x_s - gamma * (x_s - x); x_s = x_s - gamma * (x_s - x);
normr = r_s.norm(); normr = r_s.StableNorm();
} }
iter++; iter++;
} // end while } // end while
if (smoothing) if (smoothing) {
{
x = x_s; x = x_s;
} }
relres = normr / normb; relres = normr / normb;
@ -276,21 +252,19 @@ namespace Eigen
template <typename MatrixType_, typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> > template <typename MatrixType_, typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
class IDRS; class IDRS;
namespace internal namespace internal {
{
template <typename MatrixType_, typename Preconditioner_> template <typename MatrixType_, typename Preconditioner_>
struct traits<Eigen::IDRS<MatrixType_, Preconditioner_> > struct traits<Eigen::IDRS<MatrixType_, Preconditioner_> > {
{
typedef MatrixType_ MatrixType; typedef MatrixType_ MatrixType;
typedef Preconditioner_ Preconditioner; typedef Preconditioner_ Preconditioner;
}; };
} // namespace internal } // namespace internal
/** \ingroup IterativeLinearSolvers_Module /** \ingroup IterativeLinearSolvers_Module
* \brief The Induced Dimension Reduction method (IDR(s)) is a short-recurrences Krylov method for sparse square problems. * \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. * 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 * he Induced Dimension Reduction method, IDR(), is a robust and efficient short-recurrence Krylov subspace method for
@ -330,9 +304,7 @@ namespace Eigen
* \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
*/ */
template <typename MatrixType_, typename Preconditioner_> template <typename MatrixType_, typename Preconditioner_>
class IDRS : public IterativeSolverBase<IDRS<MatrixType_, Preconditioner_> > class IDRS : public IterativeSolverBase<IDRS<MatrixType_, Preconditioner_> > {
{
public: public:
typedef MatrixType_ MatrixType; typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
@ -366,9 +338,8 @@ namespace Eigen
matrix A, or modify a copy of A. matrix A, or modify a copy of A.
*/ */
template <typename MatrixDerived> template <typename MatrixDerived>
explicit IDRS(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_S(4), m_smoothing(false), explicit IDRS(const EigenBase<MatrixDerived>& A)
m_angle(RealScalar(0.7)), m_residual(false) {} : Base(A.derived()), m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {}
/** \internal */ /** \internal */
/** Loops over the number of columns of b and does the following: /** Loops over the number of columns of b and does the following:
@ -376,21 +347,19 @@ namespace Eigen
2. Calls the function that has the core solver routine 2. Calls the function that has the core solver routine
*/ */
template <typename Rhs, typename Dest> template <typename Rhs, typename Dest>
void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const {
{
m_iterations = Base::maxIterations(); m_iterations = Base::maxIterations();
m_error = Base::m_tolerance; 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); 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; 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*/ /** Sets the parameter S, indicating the dimension of the shadow space. Default is 4*/
void setS(Index S) void setS(Index S) {
{ if (S < 1) {
if (S < 1)
{
S = 4; S = 4;
} }
@ -403,10 +372,7 @@ namespace Eigen
operations. Although monotonic decrease of the residual norms is a operations. Although monotonic decrease of the residual norms is a
desirable property, the rate of convergence of the unsmoothed process and desirable property, the rate of convergence of the unsmoothed process and
the smoothed process is basically the same. Default is off */ the smoothed process is basically the same. Default is off */
void setSmoothing(bool smoothing) void setSmoothing(bool smoothing) { m_smoothing = smoothing; }
{
m_smoothing=smoothing;
}
/** The angle must be a real scalar. In IDR(s), a value for the /** 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 iteration parameter omega must be chosen in every s+1th step. The most
@ -418,19 +384,12 @@ namespace Eigen
between accurate computations and reduction of the residual norm. The between accurate computations and reduction of the residual norm. The
parameter angle =0.7 (”maintaining the convergence strategy”) parameter angle =0.7 (”maintaining the convergence strategy”)
results in such a compromise. */ results in such a compromise. */
void setAngle(RealScalar angle) void setAngle(RealScalar angle) { m_angle = angle; }
{
m_angle=angle;
}
/** The parameter replace is a logical that determines whether a /** The parameter replace is a logical that determines whether a
residual replacement strategy is employed to increase the accuracy of the residual replacement strategy is employed to increase the accuracy of the
solution. */ solution. */
void setResidualUpdate(bool update) void setResidualUpdate(bool update) { m_residual = update; }
{
m_residual=update;
}
}; };
} // namespace Eigen } // namespace Eigen