00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
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
00245 for(int iter=0; conv==SolverControl::iterate; iter++)
00246 {
00247
00248
00249 A.vmult(r,x);
00250 r.sadd(-1.,1.,b);
00251 precondition.vmult(d,r);
00252
00253
00254
00255
00256
00257
00258 conv = this->control().check (iter, criterion());
00259
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
00275 this->memory.free(Vr);
00276 this->memory.free(Vd);
00277 deallog.pop();
00278
00279
00280
00281 if (this->control().last_check() != SolverControl::success)
00282 throw SolverControl::NoConvergence (this->control().last_step(),
00283 this->control().last_value());
00284
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
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
00307 for(int iter=0; conv==SolverControl::iterate; iter++)
00308 {
00309
00310
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
00332 this->memory.free(Vr);
00333 this->memory.free(Vd);
00334 deallog.pop();
00335
00336
00337 if (this->control().last_check() != SolverControl::success)
00338 throw SolverControl::NoConvergence (this->control().last_step(),
00339 this->control().last_value());
00340
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