From db1a493e01cf4dcc29264d9df7bdec72b0b9932e Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 26 Mar 1999 10:50:44 +0000 Subject: [PATCH] First cleaned-up version. git-svn-id: https://svn.dealii.org/trunk@1055 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/solver_gmres.h | 239 +++++++------------------ 1 file changed, 67 insertions(+), 172 deletions(-) diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index 183107e714..8ed068e451 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -106,115 +106,6 @@ SolverPGMRES::givens_rotation (Vector& h, Vector& b, -/* -template -inline -Solver::ReturnState -SolverPGMRES::solve (const Matrix& A, - Vector& x, - const Vector& b) -{ - // this code was written by the fathers of - // DEAL. I take absolutely no guarantees - // for any failures or airplane-explosions - // or nuclear wars or whatever resulting - // from this code. I tried to clean a bit, - // but whoever wrote this code should get - // stone, IMHO! (WB) - - int kmax = n_tmp_vectors; - FullMatrix H(kmax+1, kmax), H1(kmax+1, kmax); - - ::Vector y(kmax), b0(kmax+1); - int i,k; - - SolverControl::State conv=SolverControl::iterate; - - double rho,beta; - - // allocate an array of n_tmp_vectors - // temporary vectors from the VectorMemory - // object - vector tmp_vectors (n_tmp_vectors, 0); - for (unsigned int tmp=0; tmpreinit (x.size()); - }; - - - A.residual(*tmp_vectors[0],x,b); - - rho = tmp_vectors[0]->l2_norm(); - beta = rho; - - tmp_vectors[0]->scale (1./rho); - - for (k=0 ; k0) H.fill(H1); - - for (i=0 ; i<=k ; i++) - { - H(i,k) = *tmp_vectors[k+1] * *tmp_vectors[i]; - tmp_vectors[k+1]->add(-H(i,k),*tmp_vectors[i]); - } - - double s = tmp_vectors[k+1]->l2_norm(); - H(k+1,k) = s; - - // Re-orthogonalization - - //printf("\n"); - for (i=0 ; i<=k ; i++) - { - double htmp = *tmp_vectors[k+1] * *tmp_vectors[i]; - //printf(" %e ",htmp); - H(i,k) += htmp; - tmp_vectors[k+1]->add(-htmp,*tmp_vectors[i]); - } - //printf("\n"); - - s = tmp_vectors[k+1]->l2_norm(); - H(k+1,k) = s; - - tmp_vectors[k+1]->scale(1./s); - - // Least - Squares - - y.reinit(k+1); - b0.reinit(k+2); - b0(0) = beta; - H1 = H; - rho = H.least_squares(y,b0); - conv = control().check(k,rho); - } - - // this will miserably fail if the - // loop above was left before k=kmax-1! - for (i=0 ; i inline @@ -228,10 +119,10 @@ SolverPGMRES::solve (const Matrix& A, // for any failures or airplane-explosions // or nuclear wars or whatever resulting // from this code. I tried to clean a bit, - // but whoever wrote this code should get - // stone, IMHO! (WB) + // but whoever wrote this code in the first + // place should get stoned, IMHO! (WB) - int kmax = n_tmp_vectors-1; + const unsigned int kmax = n_tmp_vectors-1; // allocate an array of n_tmp_vectors // temporary vectors from the VectorMemory // object @@ -246,71 +137,82 @@ SolverPGMRES::solve (const Matrix& A, // int k0 = info.usediter(); int k0 = 0; - + // matrix used for the orthogonalization + // process later FullMatrix H(kmax+1, kmax); + + // some additional vectors, also used + // in the orthogonalization ::Vector gamma(kmax+1), ci(kmax), si(kmax), h(kmax); - int i,k,reached=0,dim; - int left_precondition = 1; - Vector& v = *tmp_vectors[0]; - Vector& p = *tmp_vectors[kmax]; + + unsigned int dim; + + SolverControl::State reached = SolverControl::iterate; + + // switch to determine whether we want a + // left or a right preconditioner. at + // present, left is default, but both + // ways are implemented + const bool left_precondition = true; + + // define two aliases + Vector &v = *tmp_vectors[0]; + Vector &p = *tmp_vectors[kmax]; if (left_precondition) - { - A.residual(p,x,b); - A.precondition(v,p); - } - else - { - A.residual(v,x,b); - } + { + A.residual(p,x,b); + A.precondition(v,p); + } else { + A.residual(v,x,b); + }; - double rho = sqrt(v*v); + double rho = v.l2_norm(); gamma(0) = rho; - - v.equ(1./rho,v); - for (k=0 ; k::solve (const Matrix& A, /* append vector on matrix */ - for (i=0 ; i::solve (const Matrix& A, h.reinit(dim); FullMatrix H1(dim+1,dim); - for (i=0 ; i