mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-04-24 02:29:33 +08:00
formatting
This commit is contained in:
Gael Guennebaud
parent
04665ef9e1
commit
cacbc5679d
@ -16,189 +16,191 @@ 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: 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.
|
||||
*
|
||||
*/
|
||||
* 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,
|
||||
Index &iters, const Index &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 > FMatrixType;
|
||||
|
||||
RealScalar tol = tol_error;
|
||||
const Index maxIters = iters;
|
||||
iters = 0;
|
||||
RealScalar tol = tol_error;
|
||||
const Index maxIters = iters;
|
||||
iters = 0;
|
||||
|
||||
const Index m = mat.rows();
|
||||
const Index m = mat.rows();
|
||||
|
||||
// residual and preconditioned residual
|
||||
VectorType p0 = rhs - mat*x;
|
||||
VectorType r0 = precond.solve(p0);
|
||||
// residual and preconditioned residual
|
||||
VectorType p0 = rhs - mat*x;
|
||||
VectorType r0 = precond.solve(p0);
|
||||
VectorType t(m), v(m), workspace(m);
|
||||
|
||||
const RealScalar r0Norm = r0.norm();
|
||||
const RealScalar r0Norm = r0.norm();
|
||||
|
||||
// is initial guess already good enough?
|
||||
if(r0Norm == 0) {
|
||||
tol_error=0;
|
||||
return true;
|
||||
}
|
||||
// is initial guess already good enough?
|
||||
if(r0Norm == 0)
|
||||
{
|
||||
tol_error = 0;
|
||||
return true;
|
||||
}
|
||||
|
||||
// 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);
|
||||
// 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);
|
||||
|
||||
// storage for Jacobi rotations
|
||||
std::vector < JacobiRotation < Scalar > > G(restart);
|
||||
// storage for Jacobi rotations
|
||||
std::vector < JacobiRotation < Scalar > > G(restart);
|
||||
|
||||
// generate first Householder vector
|
||||
// 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;
|
||||
RealScalar beta;
|
||||
r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta);
|
||||
w(0) = Scalar(beta);
|
||||
|
||||
for (Index k = 1; k <= restart; ++k) {
|
||||
for (Index k = 1; k <= restart; ++k)
|
||||
{
|
||||
++iters;
|
||||
|
||||
++iters;
|
||||
v = VectorType::Unit(m, k - 1);
|
||||
|
||||
v = VectorType::Unit(m, k - 1);
|
||||
// apply Householder reflections H_{1} ... H_{k-1} to v
|
||||
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 Householder reflections H_{1} ... H_{k-1} to v
|
||||
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;
|
||||
t=mat*v;
|
||||
v=precond.solve(t);
|
||||
|
||||
// apply matrix M to v: v = mat * v;
|
||||
t=mat*v;
|
||||
v=precond.solve(t);
|
||||
// apply Householder reflections H_{k-1} ... H_{1} to v
|
||||
for (Index 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
|
||||
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 (k <= restart) {
|
||||
|
||||
// generate new Householder vector
|
||||
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);
|
||||
v.tail(m - k).makeHouseholder(Hk_tail, tau.coeffRef(k), beta);
|
||||
|
||||
// apply Householder reflection H_{k} to v
|
||||
v.tail(m - k).applyHouseholderOnTheLeft(Hk_tail, tau.coeffRef(k), workspace.data());
|
||||
// apply Householder reflection H_{k} to v
|
||||
v.tail(m - k).applyHouseholderOnTheLeft(Hk_tail, tau.coeffRef(k), workspace.data());
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
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());
|
||||
}
|
||||
}
|
||||
|
||||
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());
|
||||
}
|
||||
}
|
||||
if (k<m && v(k) != (Scalar) 0)
|
||||
{
|
||||
// determine next Givens rotation
|
||||
G[k - 1].makeGivens(v(k - 1), v(k));
|
||||
|
||||
if (k<m && v(k) != (Scalar) 0) {
|
||||
// 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());
|
||||
}
|
||||
|
||||
// determine next Givens rotation
|
||||
G[k - 1].makeGivens(v(k - 1), v(k));
|
||||
// insert coefficients into upper matrix triangle
|
||||
H.col(k - 1).head(k) = v.head(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());
|
||||
}
|
||||
bool stop = (k==m || abs(w(k)) < tol * r0Norm || iters == maxIters);
|
||||
|
||||
// insert coefficients into upper matrix triangle
|
||||
H.col(k - 1).head(k) = v.head(k);
|
||||
if (stop || k == restart)
|
||||
{
|
||||
// solve upper triangular system
|
||||
VectorType y = w.head(k);
|
||||
H.topLeftCorner(k, k).template triangularView <Upper>().solveInPlace(y);
|
||||
|
||||
bool stop=(k==m || abs(w(k)) < tol * r0Norm || iters == maxIters);
|
||||
// use Horner-like scheme to calculate solution vector
|
||||
VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
|
||||
|
||||
if (stop || k == restart) {
|
||||
// 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());
|
||||
|
||||
// solve upper triangular system
|
||||
VectorType y = w.head(k);
|
||||
H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y);
|
||||
for (Index 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());
|
||||
}
|
||||
|
||||
// use Horner-like scheme to calculate solution vector
|
||||
VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
|
||||
x += x_new;
|
||||
|
||||
// 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());
|
||||
if (stop)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
else
|
||||
{
|
||||
k=0;
|
||||
|
||||
for (Index 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());
|
||||
}
|
||||
// reset data for restart
|
||||
p0 = rhs - mat*x;
|
||||
r0 = precond.solve(p0);
|
||||
|
||||
x += x_new;
|
||||
// clear Hessenberg matrix and Householder data
|
||||
H = FMatrixType::Zero(m, restart + 1);
|
||||
w = VectorType::Zero(restart + 1);
|
||||
tau = VectorType::Zero(restart + 1);
|
||||
|
||||
if (stop) {
|
||||
return true;
|
||||
} else {
|
||||
k=0;
|
||||
// generate first Householder vector
|
||||
r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta);
|
||||
w(0)=(Scalar) beta;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// reset data for restart
|
||||
p0 = rhs - mat*x;
|
||||
r0 = precond.solve(p0);
|
||||
|
||||
// clear Hessenberg matrix and Householder data
|
||||
H = FMatrixType::Zero(m, restart + 1);
|
||||
w = VectorType::Zero(restart + 1);
|
||||
tau = VectorType::Zero(restart + 1);
|
||||
|
||||
// generate first Householder vector
|
||||
r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta);
|
||||
w(0)=(Scalar) beta;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
return false;
|
||||
return false;
|
||||
|
||||
}
|
||||
|
||||
@ -314,8 +316,8 @@ public:
|
||||
failed = true;
|
||||
}
|
||||
m_info = failed ? NumericalIssue
|
||||
: m_error <= Base::m_tolerance ? Success
|
||||
: NoConvergence;
|
||||
: m_error <= Base::m_tolerance ? Success
|
||||
: NoConvergence;
|
||||
m_isInitialized = true;
|
||||
}
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user