From: cvs Date: Wed, 13 Oct 1999 20:14:49 +0000 (+0000) Subject: First version of QMRS, do not rely on it! X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c1d1ef55baeda15b8e536d44995aedad5cad8721;p=dealii-svn.git First version of QMRS, do not rely on it! git-svn-id: https://svn.dealii.org/trunk@1768 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/solver_qmrs.h b/deal.II/lac/include/lac/solver_qmrs.h new file mode 100644 index 0000000000..382682343e --- /dev/null +++ b/deal.II/lac/include/lac/solver_qmrs.h @@ -0,0 +1,243 @@ +/*---------------------------- solver_qmrs.h ---------------------------*/ +/* $Id$ */ +#ifndef __lac__solver_qmrs_H +#define __lac__solver_qmrs_H +/*---------------------------- solver_qmrs.h ---------------------------*/ + + + +#include +#include +#include +#include + + +/** + * QMRS method. + * + * The QMRS method is supposed to solve symmetric indefinite linear + * systems with symmetric, not necessarily definite preconditioners. + * This version of QMRS is adapted from + * Freund/Nachtigal: Software for simplified Lanczos and QMR + * algorithms, Appl. Num. Math. 19 (1995), pp. 319-341 + * + * Like all other solver classes, this class has a local structure called + * #AdditionalData# which is used to pass additional parameters to the + * solver, like damping parameters or the number of temporary vectors. We + * use this additional structure instead of passing these values directly + * to the constructor because this makes the use of the #SolverSelector# and + * other classes much easier and guarantees that these will continue to + * work even if number or type of the additional parameters for a certain + * solver changes. + * + * However, since the QMRS method does not need additional data, the respective + * structure is empty and does not offer any functionality. The constructor + * has a default argument, so you may call it without the additional + * parameter. + * + * @author Guido Kanschat, 1999 + */ +template +class SolverQMRS : public Solver +{ + public: + /** + * Standardized data struct to + * pipe additional data to the + * solver. This solver does not + * need additional data. + */ + struct AdditionalData {}; + + /** + * Constructor. + */ + SolverQMRS (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data=AdditionalData()); + + /** + * Solver method. + */ + template + typename Solver::ReturnState + solve (const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner& precondition); + + protected: + /** + * Implementation of the computation of + * the norm of the residual. + */ + virtual long double criterion(); + + /** + * Temporary vectors, allocated through + * the #VectorMemory# object at the start + * of the actual solution process and + * deallocated at the end. + */ + Vector *Vv; + Vector *Vp; + Vector *Vq; + Vector *Vt; + Vector *Vd; + + /** + * Within the iteration loop, the + * square of the residual vector is + * stored in this variable. The + * function #criterion# uses this + * variable to compute the convergence + * value, which in this class is the + * norm of the residual vector and thus + * the square root of the #res2# value. + */ + long double res2; +}; + + + + +/*------------------------- Implementation ----------------------------*/ + + +template +SolverQMRS::SolverQMRS(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &) : + Solver(cn,mem) {}; + + +template +long double +SolverQMRS::criterion() +{ + return sqrt(res2); +}; + + + +template +template +typename Solver::ReturnState +SolverQMRS::solve (const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner& precondition) +{ + SolverControl::State conv=SolverControl::iterate; + + deallog.push("QMRS"); + + // Memory allocation + Vv = memory.alloc(); + Vp = memory.alloc(); + Vq = memory.alloc(); + Vt = memory.alloc(); + Vd = memory.alloc(); + // define some aliases for simpler access + Vector& v = *Vv; + Vector& p = *Vp; + Vector& q = *Vq; + Vector& t = *Vt; + Vector& d = *Vd; + // resize the vectors, but do not set + // the values since they'd be overwritten + // soon anyway. + v.reinit(x.size(), true); + p.reinit(x.size(), true); + q.reinit(x.size(), true); + t.reinit(x.size(), true); + // This vector wants to be zero. + d.reinit(x.size()); + + int it=0; + + double tau, rho, theta=0, sigma, alpha, psi, theta_old, rho_old, beta; + + double res = A.residual(v,x,b); + conv = control().check(0,res); + if (conv) + { + memory.free(Vv); + memory.free(Vp); + memory.free(Vq); + memory.free(Vt); + memory.free(Vd); + deallog.pop(); + return success; + }; + + // Step 0 + p.equ(1.,v); + precondition(q,p); + + tau = v.norm_sqr(); + rho = q*v; + + while (conv == SolverControl::iterate) + { + it++; + // Step 1 + A.vmult(t,p); + // Step 2 + sigma = q*t; +// if (fabs(sigma) < ??) +//TODO: Breakdown criteria here and below + + // Step 3 + alpha = rho/sigma; + v.add(-alpha,t); + // Step 4 + theta_old = theta; + theta = v*v/tau; + psi = 1.*(1.+theta); + tau *= theta*psi; + + d.sadd(psi*theta_old, psi*alpha, p); + x.add(d); + // Step 5 + res = sqrt((it+1)*tau); + conv = control().check(it,res); + if (conv) break; + // Step 6 +// if (fabs(rho) < ??) + // Step 7 + rho_old = rho; + precondition(q,v); + rho = q*v; + + beta = rho/rho_old; + p.sadd(beta,1.,v); + precondition(q,p); + }; + + + // Deallocate Memory + + memory.free(Vv); + memory.free(Vp); + memory.free(Vq); + memory.free(Vt); + memory.free(Vd); + + // Output + + deallog.pop(); + + if (conv == SolverControl::failure) + return exceeded; + else + return success; +}; + + + + +/*---------------------------- solver_qmrs.h ---------------------------*/ +/* end of #ifndef __solver_qmrs_H */ +#endif +/*---------------------------- solver_qmrs.h ---------------------------*/