GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-09 22:51:51 +08:00
Code Issues Packages Projects Releases Wiki Activity

Backport MINRES fixes to 3.2

Browse Source
This commit is contained in:
Gael Guennebaud 2015-02-10 19:21:41 +01:00
parent 7b35b4cacc
commit 91953d2d37
1 changed files with 32 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
unsupported/Eigen/src/IterativeSolvers/MINRES.h
View File

@ -37,22 +37,31 @@ namespace Eigen {
typedef typename Dest::Scalar Scalar; typedef typename Dest::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,1> VectorType; typedef Matrix<Scalar,Dynamic,1> VectorType;
// Check for zero rhs
const RealScalar rhsNorm2(rhs.squaredNorm());
if(rhsNorm2 == 0)
{
x.setZero();
iters = 0;
tol_error = 0;
return;
}
// initialize // initialize
const int maxIters(iters); // initialize maxIters to iters const int maxIters(iters); // initialize maxIters to iters
const int N(mat.cols()); // the size of the matrix const int N(mat.cols()); // the size of the matrix
const RealScalar rhsNorm2(rhs.squaredNorm());
const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2) const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
// Initialize preconditioned Lanczos // Initialize preconditioned Lanczos
// VectorType v_old(N); // will be initialized inside loop VectorType v_old(N); // will be initialized inside loop
VectorType v( VectorType::Zero(N) ); //initialize v VectorType v( VectorType::Zero(N) ); //initialize v
VectorType v_new(rhs-mat*x); //initialize v_new VectorType v_new(rhs-mat*x); //initialize v_new
RealScalar residualNorm2(v_new.squaredNorm()); RealScalar residualNorm2(v_new.squaredNorm());
// VectorType w(N); // will be initialized inside loop VectorType w(N); // will be initialized inside loop
VectorType w_new(precond.solve(v_new)); // initialize w_new VectorType w_new(precond.solve(v_new)); // initialize w_new
// RealScalar beta; // will be initialized inside loop // RealScalar beta; // will be initialized inside loop
RealScalar beta_new2(v_new.dot(w_new)); RealScalar beta_new2(v_new.dot(w_new));
eigen_assert(beta_new2 >= 0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE"); eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
RealScalar beta_new(sqrt(beta_new2)); RealScalar beta_new(sqrt(beta_new2));
const RealScalar beta_one(beta_new); const RealScalar beta_one(beta_new);
v_new /= beta_new; v_new /= beta_new;
@ -62,14 +71,14 @@ namespace Eigen {
RealScalar c_old(1.0); RealScalar c_old(1.0);
RealScalar s(0.0); // the sine of the Givens rotation RealScalar s(0.0); // the sine of the Givens rotation
RealScalar s_old(0.0); // the sine of the Givens rotation RealScalar s_old(0.0); // the sine of the Givens rotation
// VectorType p_oold(N); // will be initialized in loop VectorType p_oold(N); // will be initialized in loop
VectorType p_old(VectorType::Zero(N)); // initialize p_old=0 VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
VectorType p(p_old); // initialize p=0 VectorType p(p_old); // initialize p=0
RealScalar eta(1.0); RealScalar eta(1.0);
iters = 0; // reset iters iters = 0; // reset iters
while ( iters < maxIters ){ while ( iters < maxIters )
{
// Preconditioned Lanczos // Preconditioned Lanczos
/* Note that there are 4 variants on the Lanczos algorithm. These are /* Note that there are 4 variants on the Lanczos algorithm. These are
* described in Paige, C. C. (1972). Computational variants of * described in Paige, C. C. (1972). Computational variants of
@ -81,17 +90,17 @@ namespace Eigen {
* A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987). * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
*/ */
const RealScalar beta(beta_new); 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 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 // const VectorType v_old(v); // NOT SURE IF CREATING v_old EVERY ITERATION IS EFFICIENT
v = v_new; // update v = v_new; // update
// w = w_new; // update w = w_new; // update
const VectorType w(w_new); // NOT SURE IF CREATING w EVERY ITERATION IS EFFICIENT // 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 v_new.noalias() = mat*w - beta*v_old; // compute v_new
const RealScalar alpha = v_new.dot(w); const RealScalar alpha = v_new.dot(w);
v_new -= alpha*v; // overwrite v_new v_new -= alpha*v; // overwrite v_new
w_new = precond.solve(v_new); // overwrite w_new w_new = precond.solve(v_new); // overwrite w_new
beta_new2 = v_new.dot(w_new); // compute beta_new beta_new2 = v_new.dot(w_new); // compute beta_new
eigen_assert(beta_new2 >= 0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE"); eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
beta_new = sqrt(beta_new2); // compute beta_new beta_new = sqrt(beta_new2); // compute beta_new
v_new /= beta_new; // overwrite v_new for next iteration v_new /= beta_new; // overwrite v_new for next iteration
w_new /= beta_new; // overwrite w_new for next iteration w_new /= beta_new; // overwrite w_new for next iteration
@ -107,21 +116,28 @@ namespace Eigen {
s=beta_new/r1; // new sine s=beta_new/r1; // new sine
// Update solution // Update solution
// p_oold = p_old; p_oold = p_old;
const VectorType p_oold(p_old); // NOT SURE IF CREATING p_oold EVERY ITERATION IS EFFICIENT // const VectorType p_oold(p_old); // NOT SURE IF CREATING p_oold EVERY ITERATION IS EFFICIENT
p_old = p; p_old = p;
p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED? p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED?
x += beta_one*c*eta*p; x += beta_one*c*eta*p;
/* Update the squared residual. Note that this is the estimated residual.
The real residual |Ax-b|^2 may be slightly larger */
residualNorm2 *= s*s; residualNorm2 *= s*s;
if ( residualNorm2 < threshold2){ if ( residualNorm2 < threshold2)
{
break; break;
} }
eta=-s*eta; // update eta eta=-s*eta; // update eta
iters++; // increment iteration number (for output purposes) iters++; // increment iteration number (for output purposes)
} }
tol_error = std::sqrt(residualNorm2 / rhsNorm2); // return error. Note that this is the estimated error. The real error |Ax-b|/|b| may be slightly larger
/* Compute error. Note that this is the estimated error. The real
error |Ax-b|/|b| may be slightly larger */
tol_error = std::sqrt(residualNorm2 / rhsNorm2);
} }
} }