include/deal.II/lac/solver_richardson.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: solver_richardson.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2010, 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_richardson_h
00013 #define __deal2__solver_richardson_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/logstream.h>
00018 #include <deal.II/lac/solver.h>
00019 #include <deal.II/lac/solver_control.h>
00020 #include <deal.II/base/subscriptor.h>
00021 
00022 DEAL_II_NAMESPACE_OPEN
00023 
00026 
00049 template <class VECTOR = Vector<double> >
00050 class SolverRichardson : public Solver<VECTOR>
00051 {
00052   public:
00058     struct AdditionalData
00059     {
00065         AdditionalData (const double omega                       = 1,
00066                         const bool   use_preconditioned_residual = false);
00067 
00071         double omega;
00072 
00076         bool use_preconditioned_residual;
00077 
00078     };
00079 
00083     SolverRichardson (SolverControl        &cn,
00084                       VectorMemory<VECTOR> &mem,
00085                       const AdditionalData &data=AdditionalData());
00086 
00092     SolverRichardson (SolverControl        &cn,
00093                       const AdditionalData &data=AdditionalData());
00094 
00098     virtual ~SolverRichardson ();
00099 
00104     template<class MATRIX, class PRECONDITIONER>
00105     void
00106     solve (const MATRIX         &A,
00107            VECTOR               &x,
00108            const VECTOR         &b,
00109            const PRECONDITIONER &precondition);
00110 
00114     template<class MATRIX, class PRECONDITIONER>
00115     void
00116     Tsolve (const MATRIX         &A,
00117             VECTOR               &x,
00118             const VECTOR         &b,
00119             const PRECONDITIONER &precondition);
00120 
00125     void set_omega (const double om=1.);
00126 
00136     virtual void print_vectors(const unsigned int step,
00137                                const VECTOR& x,
00138                                const VECTOR& r,
00139                                const VECTOR& d) const;
00140 
00141   protected:
00146     virtual typename VECTOR::value_type criterion();
00147 
00154     VECTOR* Vr;
00164     VECTOR* Vd;
00165 
00169     AdditionalData additional_data;
00170 
00181     typename VECTOR::value_type res;
00182 };
00183 
00185 /*----------------- Implementation of the Richardson Method ------------------*/
00186 
00187 #ifndef DOXYGEN
00188 
00189 template <class VECTOR>
00190 inline
00191 SolverRichardson<VECTOR>::AdditionalData::
00192 AdditionalData (const double omega,
00193                 const bool   use_preconditioned_residual)
00194                 :
00195                 omega(omega),
00196                 use_preconditioned_residual(use_preconditioned_residual)
00197 {}
00198 
00199 
00200 template <class VECTOR>
00201 SolverRichardson<VECTOR>::SolverRichardson(SolverControl &cn,
00202                                            VectorMemory<VECTOR> &mem,
00203                                            const AdditionalData &data)
00204                 :
00205                 Solver<VECTOR> (cn,mem),
00206                 additional_data(data)
00207 {}
00208 
00209 
00210 
00211 template <class VECTOR>
00212 SolverRichardson<VECTOR>::SolverRichardson(SolverControl &cn,
00213                                            const AdditionalData &data)
00214                 :
00215                 Solver<VECTOR> (cn),
00216                 additional_data(data)
00217 {}
00218 
00219 
00220 
00221 template <class VECTOR>
00222 SolverRichardson<VECTOR>::~SolverRichardson()
00223 {}
00224 
00225 
00226 template <class VECTOR>
00227 template <class MATRIX, class PRECONDITIONER>
00228 void
00229 SolverRichardson<VECTOR>::solve (const MATRIX         &A,
00230                                  VECTOR               &x,
00231                                  const VECTOR         &b,
00232                                  const PRECONDITIONER &precondition)
00233 {
00234   SolverControl::State conv=SolverControl::iterate;
00235 
00236                                    // Memory allocation
00237   Vr  = this->memory.alloc(); VECTOR& r  = *Vr; r.reinit(x);
00238   Vd  = this->memory.alloc(); VECTOR& d  = *Vd; d.reinit(x);
00239 
00240   deallog.push("Richardson");
00241 
00242   try
00243     {
00244                                        // Main loop
00245       for(int iter=0; conv==SolverControl::iterate; iter++)
00246         {
00247                                            // Do not use residual,
00248                                            // but do it in 2 steps
00249           A.vmult(r,x);
00250           r.sadd(-1.,1.,b);
00251           precondition.vmult(d,r);
00252 
00253                                            // The required norm of the
00254                                            // (preconditioned)
00255                                            // residual is computed in
00256                                            // criterion() and stored
00257                                            // in res.
00258           conv = this->control().check (iter, criterion());
00259 //        conv = this->control().check (iter, std::sqrt(A.matrix_norm_square(r)));
00260           if (conv != SolverControl::iterate)
00261             break;
00262 
00263           x.add(additional_data.omega,d);
00264           print_vectors(iter,x,r,d);
00265         }
00266     }
00267   catch (...)
00268     {
00269       this->memory.free(Vr);
00270       this->memory.free(Vd);
00271       deallog.pop();
00272       throw;
00273     }
00274                                    // Deallocate Memory
00275   this->memory.free(Vr);
00276   this->memory.free(Vd);
00277   deallog.pop();
00278 
00279                                    // in case of failure: throw
00280                                    // exception
00281   if (this->control().last_check() != SolverControl::success)
00282     throw SolverControl::NoConvergence (this->control().last_step(),
00283                                         this->control().last_value());
00284                                    // otherwise exit as normal
00285 }
00286 
00287 
00288 template <class VECTOR>
00289 template <class MATRIX, class PRECONDITIONER>
00290 void
00291 SolverRichardson<VECTOR>::Tsolve (const MATRIX         &A,
00292                                   VECTOR               &x,
00293                                   const VECTOR         &b,
00294                                   const PRECONDITIONER &precondition)
00295 {
00296   SolverControl::State conv=SolverControl::iterate;
00297 
00298                                    // Memory allocation
00299   Vr  = this->memory.alloc(); VECTOR& r  = *Vr; r.reinit(x);
00300   Vd  =this-> memory.alloc(); VECTOR& d  = *Vd; d.reinit(x);
00301 
00302   deallog.push("RichardsonT");
00303 
00304   try
00305     {
00306                                        // Main loop
00307       for(int iter=0; conv==SolverControl::iterate; iter++)
00308         {
00309                                            // Do not use Tresidual,
00310                                            // but do it in 2 steps
00311           A.Tvmult(r,x);
00312           r.sadd(-1.,1.,b);
00313           precondition.Tvmult(d,r);
00314 
00315           conv = this->control().check (iter, criterion());
00316           if (conv != SolverControl::iterate)
00317             break;
00318 
00319           x.add(additional_data.omega,d);
00320           print_vectors(iter,x,r,d);
00321         }
00322     }
00323   catch (...)
00324     {
00325       this->memory.free(Vr);
00326       this->memory.free(Vd);
00327       deallog.pop();
00328       throw;
00329     }
00330 
00331                                    // Deallocate Memory
00332   this->memory.free(Vr);
00333   this->memory.free(Vd);
00334   deallog.pop();
00335                                    // in case of failure: throw
00336                                    // exception
00337   if (this->control().last_check() != SolverControl::success)
00338     throw SolverControl::NoConvergence (this->control().last_step(),
00339                                         this->control().last_value());
00340                                    // otherwise exit as normal
00341 }
00342 
00343 
00344 template <class VECTOR>
00345 void
00346 SolverRichardson<VECTOR>::print_vectors(const unsigned int,
00347                                         const VECTOR&,
00348                                         const VECTOR&,
00349                                         const VECTOR&) const
00350 {}
00351 
00352 
00353 
00354 template <class VECTOR>
00355 inline typename VECTOR::value_type
00356 SolverRichardson<VECTOR>::criterion()
00357 {
00358   if (!additional_data.use_preconditioned_residual)
00359     res = Vr->l2_norm();
00360   else
00361     res = Vd->l2_norm();
00362   return res;
00363 }
00364 
00365 
00366 template <class VECTOR>
00367 inline void
00368 SolverRichardson<VECTOR>::set_omega (const double om)
00369 {
00370   additional_data.omega=om;
00371 }
00372 
00373 #endif // DOXYGEN
00374 
00375 DEAL_II_NAMESPACE_CLOSE
00376 
00377 #endif
00378 
 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