17 #ifndef __deal2__solver_minres_h
18 #define __deal2__solver_minres_h
21 #include <deal.II/base/config.h>
22 #include <deal.II/lac/solver.h>
23 #include <deal.II/lac/solver_control.h>
24 #include <deal.II/base/logstream.h>
26 #include <deal.II/base/subscriptor.h>
62 template <
class VECTOR = Vector<
double> >
100 template<
class MATRIX,
class PRECONDITIONER>
105 const PRECONDITIONER &precondition);
134 const VECTOR &d)
const;
143 VECTOR *Vm0, *Vm1, *Vm2;
164 template<
class VECTOR>
167 const AdditionalData &)
174 template<
class VECTOR>
176 const AdditionalData &)
182 template<
class VECTOR>
188 template<
class VECTOR>
196 template<
class VECTOR>
201 const VECTOR &)
const
206 template<
class VECTOR>
207 template<
class MATRIX,
class PRECONDITIONER>
212 const PRECONDITIONER &precondition)
216 deallog.
push(
"minres");
219 Vu0 = this->memory.alloc();
220 Vu1 = this->memory.alloc();
221 Vu2 = this->memory.alloc();
222 Vv = this->memory.alloc();
223 Vm0 = this->memory.alloc();
224 Vm1 = this->memory.alloc();
225 Vm2 = this->memory.alloc();
227 typedef VECTOR *vecptr;
228 vecptr u[3] = {Vu0, Vu1, Vu2};
229 vecptr m[3] = {Vm0, Vm1, Vm2};
234 u[0]->reinit(b,
true);
235 u[1]->reinit(b,
true);
236 u[2]->reinit(b,
true);
237 m[0]->reinit(b,
true);
238 m[1]->reinit(b,
true);
239 m[2]->reinit(b,
true);
243 double delta[3] = { 0, 0, 0 };
244 double f[2] = { 0, 0 };
245 double e[2] = { 0, 0 };
269 precondition.vmult (v,*u[1]);
271 delta[1] = v * (*u[1]);
273 Assert (delta[1]>=0, ExcPreconditionerNotDefinite());
275 r0 = std::sqrt(delta[1]);
285 conv = this->control().check(0,r_l2);
290 v *= 1./std::sqrt(delta[1]);
295 u[2]->add (-std::sqrt(delta[1]/delta[0]), *u[0]);
298 u[2]->add (-gamma / std::sqrt(delta[1]), *u[1]);
304 precondition.vmult(v,*u[2]);
306 delta[2] = v * (*u[2]);
308 Assert (delta[2]>=0, ExcPreconditionerNotDefinite());
313 e[1] = std::sqrt(delta[2]);
317 d_ = s * e[0] - c * gamma;
318 e[0] = c * e[0] + s * gamma;
319 f[1] = s * std::sqrt(delta[2]);
320 e[1] = -c * std::sqrt(delta[2]);
323 d = std::sqrt (d_*d_ + delta[2]);
325 if (j>1) tau *= s / c;
329 s = std::sqrt(delta[2]) / d;
334 m[0]->add (-e[0], *m[1]);
336 m[0]->add (-f[0], *m[2]);
339 r_l2 *= std::fabs(s);
341 conv = this->control().check(j,r_l2);
373 this->memory.free(Vu0);
374 this->memory.free(Vu1);
375 this->memory.free(Vu2);
376 this->memory.free(Vv);
377 this->memory.free(Vm0);
378 this->memory.free(Vm1);
379 this->memory.free(Vm2);
386 this->control().last_value()));
392 DEAL_II_NAMESPACE_CLOSE
void vmult(VECTOR &u, const VECTOR &v) const
SolverMinRes(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
virtual void print_vectors(const unsigned int step, const VECTOR &x, const VECTOR &r, const VECTOR &d) const
#define AssertThrow(cond, exc)
DeclException0(ExcPreconditionerNotDefinite)
#define Assert(cond, exc)
void swap(VectorBase &u, VectorBase &v)
void push(const std::string &text)
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
virtual double criterion()
Stop iteration, goal reached.