From: wolf Date: Fri, 26 Mar 1999 10:14:50 +0000 (+0000) Subject: Check in an initial working version of the preconditioned GMRES solver. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fb5af3cd839d9f4ed8d8c8b46bd0ee7b657c3fb1;p=dealii-svn.git Check in an initial working version of the preconditioned GMRES solver. git-svn-id: https://svn.dealii.org/trunk@1051 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 new file mode 100644 index 0000000000..183107e714 --- /dev/null +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -0,0 +1,394 @@ +/*---------------------------- solver_pgmres.h ---------------------------*/ +/* $Id$ */ +#ifndef __solver_pgmres_H +#define __solver_pgmres_H +/*---------------------------- solver_pgmres.h ---------------------------*/ + +#include +#include +#include +#include + + + +/** + * Implementation of the Restarted Preconditioned Direct Generalized + * Minimal Residual Method. The stopping criterion is the norm of the + * residual. + * + * @author Original implementation by the DEAL authors, adapted by Wolfgang Bangerth + */ +template +class SolverPGMRES : public Solver { + public: + /** + * Constructor. + */ + SolverPGMRES (SolverControl &cn, + VectorMemory &mem, + const unsigned int n_tmp_vectors) : + Solver (cn,mem), + n_tmp_vectors (n_tmp_vectors) + {}; + + /** + * Solver method. + */ + virtual ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b); + + protected: + const unsigned int n_tmp_vectors; + + /** + * Implementation of the computation of + * the norm of the residual. + */ + virtual double criterion(); + + /** + * Transformation of an upper + * Hessenberg matrix into + * tridiagonal structure by givens + * rotation of the last column + */ + void givens_rotation (Vector &h, Vector &b, + Vector &ci, Vector &si, + int col) const; +}; + + + + +/* ------------------------- Inline functions ----------------------------- */ + +template +inline +void +SolverPGMRES::givens_rotation (Vector& h, Vector& b, + Vector& ci, Vector& si, + int col) const +{ + for (int i=0 ; i +// 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) +{ + // 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 +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-1; + // 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()); + }; + +// WB +// int k0 = info.usediter(); + int k0 = 0; + + + FullMatrix H(kmax+1, kmax); + ::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]; + + if (left_precondition) + { + A.residual(p,x,b); + A.precondition(v,p); + } + else + { + A.residual(v,x,b); + } + + double rho = sqrt(v*v); + gamma(0) = rho; + + v.equ(1./rho,v); + + for (k=0 ; k H1(dim+1,dim); + + for (i=0 ; i +double +SolverPGMRES::criterion () +{ + // dummy implementation. this function is + // not needed for the present implementation + // of gmres + Assert (false, ExcInternalError()); + return 0; +}; + + + +/*---------------------------- solver_pgmres.h ---------------------------*/ +/* end of #ifndef __solver_pgmres_H */ +#endif +/*---------------------------- solver_pgmres.h ---------------------------*/