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

added preconditioner with preconditioned-Lanczos iteration

Browse Source
This commit is contained in:
giacomo po 2012-09-01 21:59:06 +02:00
parent 5f3880c5a8
commit 751501eade
1 changed files with 97 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

134
Eigen/src/IterativeLinearSolvers/MINRES.h
View File

@ -35,73 +35,133 @@ namespace Eigen {
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;
// 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 rhsNorm2(rhs.squaredNorm());
// const RealScalar threshold(tol_error); // threshold for original convergence criterion, see below // const RealScalar threshold(tol_error); // threshold for original convergence criterion, see below
const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold
VectorType v(VectorType::Zero(N));
VectorType v_hat(rhs-mat*x);
RealScalar beta(v_hat.norm()); // VectorType v(VectorType::Zero(N));
// VectorType v_hat(rhs-mat*x);
// Compute initial residual
VectorType residual(rhs-mat*x);
// Initialize preconditioned Lanczos
VectorType v_old(N); // will be initialized inside loop
VectorType v = VectorType::Zero(N); //initialize v
VectorType v_new = residual; //initialize v_new
VectorType w(N); // will be initialized inside loop
VectorType w_new = precond.solve(v_new); // initialize w_new
RealScalar beta; // will be initialized inside loop
RealScalar beta_new = sqrt(v_new.dot(w_new));
v_new /= beta_new;
w_new /= beta_new;
// RealScalar beta(v_hat.norm());
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);
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 w(VectorType::Zero(N)); VectorType p_oold(VectorType::Zero(N)); // initialize p_oold=0
VectorType w_old(w); VectorType p_old(p_oold); // initialize p_old=0
RealScalar eta(beta); VectorType p(N); // will be initialized in loop
//RealScalar eta(beta); // CHANGE THIS
RealScalar norm_rMR=beta; RealScalar norm_rMR=beta;
const RealScalar norm_r0(beta); const RealScalar norm_r0(beta);
VectorType v_old(N), Av(N), w_oold(N); // preallocate temporaty vectors used in iteration RealScalar eta(1.0);
// VectorType v_old(N), Av(N), w_oold(N); // preallocate temporaty vectors used in iteration
RealScalar residualNorm2; // not needed for original convergnce criterion RealScalar residualNorm2; // not needed for original convergnce criterion
int n = 0; int n = 0;
while ( n < maxIters ){ while ( n < maxIters ){
// Lanczos // Preconditioned Lanczos
// VectorType v_old(v); // now pre-allocated /* Note that there are 4 variants on the Lanczos algorithm. These are
v_old = v; * described in Paige, C. C. (1972). Computational variants of
v=v_hat/beta; * the Lanczos method for the eigenproblem. IMA Journal of Applied
// VectorType Av(mat*v); // now pre-allocated * Mathematics, 10(3), 373–381. The current implementation corresonds
Av = mat*v; * to the case A(2,7) in the paper. It also corresponds to
RealScalar alpha(v.transpose()*Av); * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear
v_hat=Av-alpha*v-beta*v_old; * Systems, 2003 p.173. For the preconditioned version see
RealScalar beta_old(beta); * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
beta=v_hat.norm(); */
beta = beta_new;
v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
v = v_new; // update
w = w_new; // update
v_new.noalias() = mat*w - beta*v_old; // compute v_new
const RealScalar alpha = v_new.dot(w);
v_new -= alpha*v; // overwrite v_new
w_new = precond.solve(v_new); // overwrite w_new
beta_new = sqrt(v_new.dot(w_new)); // compute beta_new
v_new /= beta_new; // overwrite v_new
w_new /= beta_new; // overwrite w_new
// QR //
RealScalar c_oold(c_old); //
c_old=c; //
RealScalar s_oold(s_old); //
s_old=s; //
RealScalar r1_hat=c_old *alpha-c_oold*s_old *beta_old; //
RealScalar r1 =std::pow(std::pow(r1_hat,2)+std::pow(beta,2),0.5); //
RealScalar r2 =s_old *alpha+c_oold*c_old*beta_old; //
RealScalar r3 =s_oold*beta_old; // // VectorType v_old(v); // now pre-allocated
// v_old = v;
// v=v_hat/beta;
//// VectorType Av(mat*v); // now pre-allocated
// Av = mat*v;
// RealScalar alpha(v.transpose()*Av);
// v_hat=Av-alpha*v-beta*v_old;
// RealScalar beta_old(beta);
// beta=v_hat.norm();
// Givens rotation // Apply QR
c=r1_hat/r1; // RealScalar c_oold(c_old); // store old-old cosine
s=beta/r1; // c_old=c; // store old cosine
// RealScalar s_oold(s_old); // store old-old sine
// s_old=s; // store old sine
// const RealScalar r1_hat=c_old *alpha-c_oold*s_old *beta_old;
// const RealScalar r1 =std::pow(std::pow(r1_hat,2)+std::pow(beta,2),0.5);
const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration
const RealScalar r3 =s_old*beta; // s, s_old, c and c_old are still from previous iteration
// update // Compute new Givens rotation
const RealScalar r1_hat=c*alpha-c_old*s*beta;
const RealScalar r1 =std::pow(std::pow(r1_hat,2)+std::pow(beta_new,2),0.5);
c_old = c; // store for next iteration
s_old = s; // store for next iteration
c=r1_hat/r1; // new cosine
s=beta/r1; // new sine
// update w
// VectorType w_oold(w_old); // now pre-allocated // VectorType w_oold(w_old); // now pre-allocated
w_oold = w_old; p_oold = p_old;
p_old = p;
w_old=w; p=(w-r2*p_old-r3*p_oold) /r1;
w=(v-r3*w_oold-r2*w_old) /r1; // update x
x += c*eta*w; x += c*eta*p;
norm_rMR *= std::fabs(s); norm_rMR *= std::fabs(s);
eta=-s*eta;
residualNorm2 = (mat*x-rhs).squaredNorm(); // DOES mat*x NEED TO BE RECOMPUTED ???? residualNorm2 = (mat*x-rhs).squaredNorm(); // DOES mat*x NEED TO BE RECOMPUTED ????
//if(norm_rMR/norm_r0 < threshold){ // original convergence criterion, does not require "mat*x" //if(norm_rMR/norm_r0 < threshold){ // original convergence criterion, does not require "mat*x"
if ( residualNorm2 < threshold2){ if ( residualNorm2 < threshold2){
break; break;
} }
n++;
eta=-s*eta; // update eta
n++; // increment iteration
} }
tol_error = std::sqrt(residualNorm2 / rhsNorm2); // return error tol_error = std::sqrt(residualNorm2 / rhsNorm2); // return error
iters = n; // return number of iterations iters = n; // return number of iterations