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>
174 template<
class VECTOR>
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");
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]);
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);
386 this->
control().last_value()));
392 DEAL_II_NAMESPACE_CLOSE
VectorMemory< VECTOR > & memory
virtual State check(const unsigned int step, const double check_value)
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)
SolverControl & control() const
Stop iteration, goal reached.
#define Assert(cond, exc)
void push(const std::string &text)
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
virtual double criterion()