From: wolf Date: Wed, 7 Apr 1999 13:48:42 +0000 (+0000) Subject: Implement restarts for GMRES. More documentation. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f7e8acac07db046ba6f976418990e7f217c0c89b;p=dealii-svn.git Implement restarts for GMRES. More documentation. git-svn-id: https://svn.dealii.org/trunk@1073 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index 703ee111b5..85f8c3baa1 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -16,17 +16,28 @@ * Minimal Residual Method. The stopping criterion is the norm of the * residual. * - * @author Original implementation by the DEAL authors, adapted by Wolfgang Bangerth + * You have to give the number of temporary vectors to the constructor which + * are to be used to do the orthogonalization. If the number of iterations + * needed to solve the problem to the given criterion, an intermediate + * solution is computed and a restart is performed. If you don't want to + * use the restarted method, you can limit the number of iterations (stated + * in the #SolverControl# object given to the constructor) to be below + * the number of temporary vectors minus three. Note the subtraction, which + * is due to the fact that three vectors are used for other purposes, so + * the number of iterations before a restart occurs is less by three than + * the total number of temporary vectors. + * + * @author Original implementation by the DEAL authors; adapted, cleaned and documented by Wolfgang Bangerth */ template -class SolverPGMRES : public Solver { +class SolverGMRES : public Solver { public: /** * Constructor. */ - SolverPGMRES (SolverControl &cn, - VectorMemory &mem, - const unsigned int n_tmp_vectors) : + SolverGMRES (SolverControl &cn, + VectorMemory &mem, + const unsigned int n_tmp_vectors) : Solver (cn,mem), n_tmp_vectors (n_tmp_vectors) {}; @@ -38,6 +49,12 @@ class SolverPGMRES : public Solver { Vector &x, const Vector &b); + DeclException1 (ExcTooFewTmpVectors, + int, + << "The number of temporary vectors you gave (" + << arg1 << ") is too small. It should be at least 10 for " + << "any results, and much more for reasonable ones."); + protected: const unsigned int n_tmp_vectors; @@ -61,14 +78,29 @@ class SolverPGMRES : public Solver { -/* ------------------------- Inline functions ----------------------------- */ +/* --------------------- Inline and template functions ------------------- */ + + +template +SolverGMRES::SolverGMRES (SolverControl &cn, + VectorMemory &mem, + const unsigned int n_tmp_vectors) : + Solver (cn,mem), + n_tmp_vectors (n_tmp_vectors) +{ + Assert (n_tmp_vectors >= 10, ExcTooFewTmpVectors (n_tmp_vectors)); +}; + + template inline void -SolverPGMRES::givens_rotation (Vector& h, Vector& b, - Vector& ci, Vector& si, - int col) const +SolverGMRES::givens_rotation (Vector &h, + Vector &b, + Vector &ci, + Vector &si, + int col) const { for (int i=0 ; i::givens_rotation (Vector& h, Vector& b, } -// // restarted method -// template -// inline int -// SolverPGMRES::solve (Matrix& A, Vector& x, Vector& b) -// { -// int reached, kmax = mem.n()-1; -// for(int j=0;;j++) -// { -// info.usediter() = j*(kmax-1); -// reached = dgmres(A,x,b,mem,info); -// if(reached) break; -// } -// if (reached<0) return 1; -// return 0; -// } - - template -inline Solver::ReturnState -SolverPGMRES::solve (const Matrix& A, - Vector& x, - const Vector& b) +SolverGMRES::solve (const Matrix& A, + Vector & x, + const Vector& b) { // this code was written by the fathers of // DEAL. I take absolutely no guarantees @@ -122,7 +136,6 @@ SolverPGMRES::solve (const Matrix& A, // but whoever wrote this code in the first // place should get stoned, IMHO! (WB) - const unsigned int kmax = n_tmp_vectors-1; // allocate an array of n_tmp_vectors // temporary vectors from the VectorMemory // object @@ -133,22 +146,26 @@ SolverPGMRES::solve (const Matrix& A, tmp_vectors[tmp]->reinit (x.size()); }; -// WB -// int k0 = info.usediter(); - int k0 = 0; + // number of the present iteration; this + // number is not reset to zero upon a + // restart + unsigned int accumulated_iterations = 0; // matrix used for the orthogonalization // process later - FullMatrix H(kmax+1, kmax); + FullMatrix H(n_tmp_vectors, n_tmp_vectors-1); // some additional vectors, also used // in the orthogonalization - ::Vector gamma(kmax+1), ci(kmax), si(kmax), h(kmax); + ::Vector gamma(n_tmp_vectors), + ci (n_tmp_vectors-1), + si (n_tmp_vectors-1), + h (n_tmp_vectors-1); unsigned int dim; - SolverControl::State reached = SolverControl::iterate; + SolverControl::State iteration_state = SolverControl::iterate; // switch to determine whether we want a // left or a right preconditioner. at @@ -158,109 +175,132 @@ SolverPGMRES::solve (const Matrix& A, // define two aliases Vector &v = *tmp_vectors[0]; - Vector &p = *tmp_vectors[kmax]; + Vector &p = *tmp_vectors[n_tmp_vectors-1]; - if (left_precondition) - { - A.residual(p,x,b); - A.precondition(v,p); - } else { - A.residual(v,x,b); - }; - - double rho = v.l2_norm(); - gamma(0) = rho; - - v.scale (1./rho); - - // inner iteration doing at most as - // many steps as there are temporary - // vectors. the number of steps actually - // been done is propagated outside - // through the #dim# variable - for (unsigned int inner_iteration=0; - inner_iteration H1(dim+1,dim); - /* residual */ - rho = fabs(gamma(dim)); - - reached = control().check (k0+inner_iteration, rho); - }; - - /* Calculate solution */ - h.reinit(dim); - FullMatrix H1(dim+1,dim); + for (unsigned int i=0; i::solve (const Matrix& A, template double -SolverPGMRES::criterion () +SolverGMRES::criterion () { // dummy implementation. this function is // not needed for the present implementation