include/deal.II/lac/solver_qmrs.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: solver_qmrs.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__solver_qmrs_h
00013 #define __deal2__solver_qmrs_h
00014 
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/lac/solver.h>
00017 #include <deal.II/lac/solver_control.h>
00018 #include <deal.II/base/logstream.h>
00019 #include <cmath>
00020 #include <deal.II/base/subscriptor.h>
00021 
00022 #include <cmath>
00023 
00024 DEAL_II_NAMESPACE_OPEN
00025 
00028 
00061 template <class VECTOR = Vector<double> >
00062 class SolverQMRS : public Solver<VECTOR>
00063 {
00064   public:
00085     struct AdditionalData
00086     {
00094         AdditionalData(bool exact_residual = false,
00095                        double breakdown=1.e-16) :
00096                         exact_residual(exact_residual),
00097                         breakdown(breakdown)
00098           {}
00099 
00103         bool exact_residual;
00104 
00108         double breakdown;
00109     };
00110 
00114     SolverQMRS (SolverControl &cn,
00115                 VectorMemory<VECTOR> &mem,
00116                 const AdditionalData &data=AdditionalData());
00117 
00123     SolverQMRS (SolverControl        &cn,
00124                 const AdditionalData &data=AdditionalData());
00125 
00130     template<class MATRIX, class PRECONDITIONER>
00131     void
00132     solve (const MATRIX         &A,
00133            VECTOR               &x,
00134            const VECTOR         &b,
00135            const PRECONDITIONER &precondition);
00136 
00146     virtual void print_vectors(const unsigned int step,
00147                                const VECTOR& x,
00148                                const VECTOR& r,
00149                                const VECTOR& d) const;
00150    protected:
00155     virtual double criterion();
00156 
00163     VECTOR *Vv;
00164     VECTOR *Vp;
00165     VECTOR *Vq;
00166     VECTOR *Vt;
00167     VECTOR *Vd;
00171     VECTOR *Vx;
00175     const VECTOR *Vb;
00176 
00187     double res2;
00191     AdditionalData additional_data;
00192   private:
00196     template<class MATRIX, class PRECONDITIONER>
00197     bool
00198     iterate(const MATRIX& A, const PRECONDITIONER& precondition);
00202     unsigned int step;
00203 };
00204 
00206 /*------------------------- Implementation ----------------------------*/
00207 
00208 #ifndef DOXYGEN
00209 
00210 template<class VECTOR>
00211 SolverQMRS<VECTOR>::SolverQMRS(SolverControl &cn,
00212                                VectorMemory<VECTOR> &mem,
00213                                const AdditionalData &data)
00214                 :
00215                 Solver<VECTOR>(cn,mem),
00216                 additional_data(data)
00217 {}
00218 
00219 
00220 
00221 template<class VECTOR>
00222 SolverQMRS<VECTOR>::SolverQMRS(SolverControl &cn,
00223                                const AdditionalData &data)
00224                 :
00225                 Solver<VECTOR>(cn),
00226                 additional_data(data)
00227 {}
00228 
00229 
00230 
00231 template<class VECTOR>
00232 double
00233 SolverQMRS<VECTOR>::criterion()
00234 {
00235   return std::sqrt(res2);
00236 }
00237 
00238 
00239 
00240 template<class VECTOR>
00241 void
00242 SolverQMRS<VECTOR>::print_vectors(const unsigned int,
00243                                   const VECTOR&,
00244                                   const VECTOR&,
00245                                   const VECTOR&) const
00246 {}
00247 
00248 
00249 
00250 template<class VECTOR>
00251 template<class MATRIX, class PRECONDITIONER>
00252 void
00253 SolverQMRS<VECTOR>::solve (const MATRIX         &A,
00254                            VECTOR               &x,
00255                            const VECTOR         &b,
00256                            const PRECONDITIONER &precondition)
00257 {
00258   deallog.push("QMRS");
00259 
00260                                    // Memory allocation
00261   Vv  = this->memory.alloc();
00262   Vp  = this->memory.alloc();
00263   Vq  = this->memory.alloc();
00264   Vt  = this->memory.alloc();
00265   Vd  = this->memory.alloc();
00266 
00267   Vx = &x;
00268   Vb = &b;
00269                                    // resize the vectors, but do not set
00270                                    // the values since they'd be overwritten
00271                                    // soon anyway.
00272   Vv->reinit(x, true);
00273   Vp->reinit(x, true);
00274   Vq->reinit(x, true);
00275   Vt->reinit(x, true);
00276 
00277   step = 0;
00278 
00279   bool state;
00280 
00281   do
00282     {
00283       if (step)
00284         deallog << "Restart step " << step << std::endl;
00285       state = iterate(A, precondition);
00286     }
00287   while (state);
00288 
00289                                    // Deallocate Memory
00290   this->memory.free(Vv);
00291   this->memory.free(Vp);
00292   this->memory.free(Vq);
00293   this->memory.free(Vt);
00294   this->memory.free(Vd);
00295 
00296                                    // Output
00297   deallog.pop();
00298 
00299                                    // in case of failure: throw
00300                                    // exception
00301   if (this->control().last_check() != SolverControl::success)
00302     throw SolverControl::NoConvergence (this->control().last_step(),
00303                                         this->control().last_value());
00304                                    // otherwise exit as normal
00305 }
00306 
00307 
00308 
00309 template<class VECTOR>
00310 template<class MATRIX, class PRECONDITIONER>
00311 bool
00312 SolverQMRS<VECTOR>::iterate(const MATRIX         &A,
00313                             const PRECONDITIONER &precondition)
00314 {
00315 /* Remark: the matrix A in the article is the preconditioned matrix.
00316  * Therefore, we have to precondition x before we compute the first residual.
00317  * In step 1 we replace p by q to avoid one preconditioning step.
00318  * There are still two steps left, making this algorithm expensive.
00319  */
00320 
00321   SolverControl::State state = SolverControl::iterate;
00322 
00323                                    // define some aliases for simpler access
00324   VECTOR& v  = *Vv;
00325   VECTOR& p  = *Vp;
00326   VECTOR& q  = *Vq;
00327   VECTOR& t  = *Vt;
00328   VECTOR& d  = *Vd;
00329   VECTOR& x  = *Vx;
00330   const VECTOR& b = *Vb;
00331 
00332   int  it=0;
00333 
00334   double tau, rho, theta=0, sigma, alpha, psi, theta_old, rho_old, beta;
00335   double res;
00336 
00337   d.reinit(x);
00338 
00339                                    // Apply right preconditioning to x
00340   precondition.vmult(q,x);
00341                                    // Preconditioned residual
00342   A.vmult(v,q);
00343   v.sadd(-1.,1.,b);
00344   res = v.l2_norm();
00345 
00346   if (this->control().check(step, res) == SolverControl::success)
00347     return false;
00348 
00349   p = v;
00350 
00351   precondition.vmult(q,p);
00352 
00353   tau = v.norm_sqr();
00354   rho = q*v;
00355 
00356   while (state == SolverControl::iterate)
00357     {
00358       step++; it++;
00359                                        // Step 1
00360       A.vmult(t,q);
00361                                        // Step 2
00362       sigma = q*t;
00363 
00364 //TODO:[?] Find a really good breakdown criterion. The absolute one detects breakdown instead of convergence
00365       if (std::fabs(sigma/rho) < additional_data.breakdown)
00366         return true;
00367                                        // Step 3
00368       alpha = rho/sigma;
00369 
00370       v.add(-alpha,t);
00371                                        // Step 4
00372       theta_old = theta;
00373       theta = v*v/tau;
00374       psi = 1./(1.+theta);
00375       tau *= theta*psi;
00376 
00377       d.sadd(psi*theta_old, psi*alpha, p);
00378       x.add(d);
00379 
00380       print_vectors(step,x,v,d);
00381                                        // Step 5
00382       if (additional_data.exact_residual)
00383         {
00384           A.vmult(q,x);
00385           q.sadd(-1.,1.,b);
00386           res = q.l2_norm();
00387         }
00388       else
00389         res = std::sqrt((it+1)*tau);
00390       state = this->control().check(step,res);
00391       if ((state == SolverControl::success)
00392       || (state == SolverControl::failure))
00393         return false;
00394                                        // Step 6
00395       if (std::fabs(rho) < additional_data.breakdown)
00396         return true;
00397                                        // Step 7
00398       rho_old = rho;
00399       precondition.vmult(q,v);
00400       rho = q*v;
00401 
00402       beta = rho/rho_old;
00403       p.sadd(beta,v);
00404       precondition.vmult(q,p);
00405     }
00406   return false;
00407 }
00408 
00409 #endif // DOXYGEN
00410 
00411 DEAL_II_NAMESPACE_CLOSE
00412 
00413 #endif
00414 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3