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

bug #1558: fix a corner case in MINRES when both v_new and w_new vanish.

Browse Source
This commit is contained in:
David Hyde 2018-07-08 22:06:38 -07:00
parent 6e654f3379
commit d908afe35f
1 changed files with 5 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -3,6 +3,7 @@
// //
// Copyright (C) 2012 Giacomo Po <gpo@ucla.edu> // Copyright (C) 2012 Giacomo Po <gpo@ucla.edu>
// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // 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 // 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 // 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"); 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;
w_new /= beta_new;
// Initialize other variables // Initialize other variables
RealScalar c(1.0); // the cosine of the Givens rotation RealScalar c(1.0); // the cosine of the Givens rotation
RealScalar c_old(1.0); RealScalar c_old(1.0);
@ -83,18 +82,18 @@ namespace Eigen {
/* 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
* the Lanczos method for the eigenproblem. IMA Journal of Applied * 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 * 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 * Systems, 2003 p.173. For the preconditioned version see
* 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 v_new /= beta_new; // overwrite v_new for next iteration
w_new /= beta_new; // overwrite w_new for next iteration
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
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
@ -102,8 +101,6 @@ namespace Eigen {
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.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
w_new /= beta_new; // overwrite w_new for next iteration
// Givens rotation // Givens rotation
const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration 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 // 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
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;
@ -286,4 +282,3 @@ namespace Eigen {
} // end namespace Eigen } // end namespace Eigen
#endif // EIGEN_MINRES_H #endif // EIGEN_MINRES_H