From: kanschat Date: Wed, 3 Nov 1999 15:09:41 +0000 (+0000) Subject: It works! X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0c8fe1165f6f0fd0c1d7f7256786382c4e096d96;p=dealii-svn.git It works! git-svn-id: https://svn.dealii.org/trunk@1820 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 index 401675c176..79cbaf0afd 100644 --- a/deal.II/lac/include/lac/solver_qmrs.h +++ b/deal.II/lac/include/lac/solver_qmrs.h @@ -4,7 +4,7 @@ #define __lac__solver_qmrs_H /*---------------------------- solver_qmrs.h ---------------------------*/ - +//TODO: Check for programming errors!!! #include #include @@ -21,6 +21,10 @@ * Freund/Nachtigal: Software for simplified Lanczos and QMR * algorithms, Appl. Num. Math. 19 (1995), pp. 319-341 * + * This version is for right preconditioning only, since then only the + * preconditioner is used: left preconditioning seems to require the + * inverse. + * * 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 @@ -116,7 +120,19 @@ class SolverQMRS : public Solver Vector *Vq; Vector *Vt; Vector *Vd; + /** + * Iteration vector. + */ + Vector *Vx; + /** + * RHS vector. + */ + const Vector *Vb; + /** + * Pointer to the matrix to be inverted. + */ + const Matrix* MA; /** * Within the iteration loop, the * square of the residual vector is @@ -132,6 +148,17 @@ class SolverQMRS : public Solver * Breakdown threshold. */ AdditionalData additional_data; + private: + /** + * The iteration loop itself. + */ + template + typename Solver::ReturnState + iterate(const Preconditioner& precondition); + /** + * The current iteration step. + */ + unsigned int step; }; @@ -162,12 +189,10 @@ template template typename Solver::ReturnState SolverQMRS::solve (const Matrix &A, - Vector &x, - const Vector &b, - const Preconditioner& precondition) + Vector &x, + const Vector &b, + const Preconditioner& precondition) { - SolverControl::State conv=SolverControl::iterate; - deallog.push("QMRS"); // Memory allocation @@ -176,64 +201,121 @@ SolverQMRS::solve (const Matrix &A, Vq = memory.alloc(); Vt = memory.alloc(); Vd = memory.alloc(); + + MA = &A; + Vx = &x; + Vb = &b; + // resize the vectors, but do not set + // the values since they'd be overwritten + // soon anyway. + Vv->reinit(x.size(), true); + Vp->reinit(x.size(), true); + Vq->reinit(x.size(), true); + Vt->reinit(x.size(), true); + + step = 0; + + ReturnState state = breakdown; + + do + { + if (step) + deallog << "Restart step " << step << endl; + state = iterate(precondition); + } + while (state == breakdown); + + // Deallocate Memory + + memory.free(Vv); + memory.free(Vp); + memory.free(Vq); + memory.free(Vt); + memory.free(Vd); + + // Output + + deallog.pop(); + + return state; +}; + +template +template +typename Solver::ReturnState +SolverQMRS::iterate(const Preconditioner& precondition) +{ +/* Remark: the matrix A in the article is the preconditioned matrix. + * Therefore, we have to precondition x before we compute the first residual. + * In step 1 we replace p by q to avoid one preconditioning step. + * There are still two steps left, making this algorithm expensive. + */ + + SolverControl::State state = SolverControl::iterate; + // 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()); - + Vector& x = *Vx; + const Vector& b = *Vb; + + const Matrix& A = *MA; + int it=0; double tau, rho, theta=0, sigma, alpha, psi, theta_old, rho_old, beta; + double res; + + d.reinit(x); - 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; - }; + // Apply right preconditioning to x + precondition(q,x); + // Preconditioned residual + res = A.residual(v, q, b); + + if (control().check(step, res) == SolverControl::success) + return ReturnState(success); - // Step 0 - p.equ(1.,v); + p = v; + precondition(q,p); tau = v.norm_sqr(); + //deallog << "tau:" << tau << endl; rho = q*v; + //deallog << "rho:" << rho << endl; - while (conv == SolverControl::iterate) + + while (state == SolverControl::iterate) { - it++; + step++; it++; // Step 1 - A.vmult(t,p); + A.vmult(t,q); // Step 2 sigma = q*t; + +//TODO: Find a really good breakdown criterion +// The absolute one detects breakdown instead of convergence if (fabs(sigma) < additional_data.breakdown) return ReturnState(breakdown); // Step 3 alpha = rho/sigma; + //deallog << "alpha:" << alpha << endl; + v.add(-alpha,t); // Step 4 theta_old = theta; theta = v*v/tau; - psi = 1.*(1.+theta); + psi = 1./(1.+theta); tau *= theta*psi; + //deallog << "psi:" << psi << endl; + //deallog << "theta:" << theta << endl; + //deallog << "tau:" << tau << endl; + d.sadd(psi*theta_old, psi*alpha, p); x.add(d); // Step 5 @@ -241,9 +323,11 @@ SolverQMRS::solve (const Matrix &A, res = A.residual(q,x,b); else res = sqrt((it+1)*tau); - conv = control().check(it,res); - if (conv) break; - + state = control().check(step,res); + if (state == SolverControl::success) + return ReturnState(success); + else if (state == SolverControl::failure) + return ReturnState(exceeded); // Step 6 if (fabs(rho) < additional_data.breakdown) return ReturnState(breakdown); @@ -253,29 +337,11 @@ SolverQMRS::solve (const Matrix &A, rho = q*v; beta = rho/rho_old; - p.sadd(beta,1.,v); + p.sadd(beta,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; -}; - + } + return ReturnState(exceeded); +}