00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__solver_relaxation_h
00013 #define __deal2__solver_relaxation_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
00045 template <class VECTOR = Vector<double> >
00046 class SolverRelaxation : public Solver<VECTOR>
00047 {
00048 public:
00055 struct AdditionalData {};
00056
00060 SolverRelaxation (SolverControl &cn,
00061 const AdditionalData &data=AdditionalData());
00062
00066 virtual ~SolverRelaxation ();
00067
00075 template<class MATRIX, class RELAXATION>
00076 void
00077 solve (const MATRIX& A,
00078 VECTOR& x,
00079 const VECTOR& b,
00080 const RELAXATION& R);
00081 };
00082
00083
00084
00085 template <class VECTOR>
00086 SolverRelaxation<VECTOR>::SolverRelaxation(SolverControl &cn,
00087 const AdditionalData &)
00088 :
00089 Solver<VECTOR> (cn)
00090 {}
00091
00092
00093
00094 template <class VECTOR>
00095 SolverRelaxation<VECTOR>::~SolverRelaxation()
00096 {}
00097
00098
00099 template <class VECTOR>
00100 template <class MATRIX, class RELAXATION>
00101 void
00102 SolverRelaxation<VECTOR>::solve (
00103 const MATRIX& A,
00104 VECTOR& x,
00105 const VECTOR& b,
00106 const RELAXATION &R)
00107 {
00108 GrowingVectorMemory<VECTOR> mem;
00109 SolverControl::State conv=SolverControl::iterate;
00110
00111
00112 typename VectorMemory<VECTOR>::Pointer Vr(mem); VECTOR& r = *Vr; r.reinit(x);
00113 typename VectorMemory<VECTOR>::Pointer Vd(mem); VECTOR& d = *Vd; d.reinit(x);
00114
00115 deallog.push("Relaxation");
00116
00117 try
00118 {
00119
00120 for(int iter=0; conv==SolverControl::iterate; iter++)
00121 {
00122
00123 A.vmult(r,x);
00124 r.sadd(-1.,1.,b);
00125
00126
00127
00128
00129
00130
00131 conv = this->control().check (iter, r.l2_norm());
00132 if (conv != SolverControl::iterate)
00133 break;
00134 R.step(x,b);
00135 }
00136 }
00137 catch (...)
00138 {
00139 deallog.pop();
00140 throw;
00141 }
00142 deallog.pop();
00143
00144
00145
00146 if (this->control().last_check() != SolverControl::success)
00147 throw SolverControl::NoConvergence (this->control().last_step(),
00148 this->control().last_value());
00149
00150 }
00151
00152
00153 DEAL_II_NAMESPACE_CLOSE
00154
00155 #endif
00156