16 #ifndef dealii__solver_minres_h 17 #define dealii__solver_minres_h 20 #include <deal.II/base/config.h> 21 #include <deal.II/lac/solver.h> 22 #include <deal.II/lac/solver_control.h> 23 #include <deal.II/base/logstream.h> 25 #include <deal.II/base/subscriptor.h> 26 #include <deal.II/base/signaling_nan.h> 28 DEAL_II_NAMESPACE_OPEN
68 template <
class VectorType = Vector<
double> >
102 template<
typename MatrixType,
typename PreconditionerType>
104 solve (
const MatrixType &A,
107 const PreconditionerType &precondition);
133 const VectorType &d)
const;
139 VectorType *
Vu0, *Vu1, *Vu2;
140 VectorType *Vm0, *Vm1, *Vm2;
157 template<
class VectorType>
162 Solver<VectorType>(cn,mem)
167 template<
class VectorType>
179 res2(numbers::signaling_nan<double>())
183 template<
class VectorType>
189 template<
class VectorType>
197 template<
class VectorType>
202 const VectorType &)
const 207 template<
class VectorType>
208 template<
typename MatrixType,
typename PreconditionerType>
213 const PreconditionerType &precondition)
215 deallog.
push(
"minres");
219 Vu1 = this->
memory.alloc();
220 Vu2 = this->
memory.alloc();
221 Vv = this->
memory.alloc();
222 Vm0 = this->
memory.alloc();
223 Vm1 = this->
memory.alloc();
224 Vm2 = this->
memory.alloc();
226 typedef VectorType *vecptr;
227 vecptr u[3] = {
Vu0, Vu1, Vu2};
228 vecptr m[3] = {Vm0, Vm1, Vm2};
233 u[0]->reinit(b,
true);
234 u[1]->reinit(b,
true);
235 u[2]->reinit(b,
true);
236 m[0]->reinit(b,
true);
237 m[1]->reinit(b,
true);
238 m[2]->reinit(b,
true);
242 double delta[3] = { 0, 0, 0 };
243 double f[2] = { 0, 0 };
244 double e[2] = { 0, 0 };
266 precondition.vmult (v,*u[1]);
268 delta[1] = v * (*u[1]);
272 r0 = std::sqrt(delta[1]);
286 v *= 1./std::sqrt(delta[1]);
291 u[2]->add (-std::sqrt(delta[1]/delta[0]), *u[0]);
293 const double gamma = *u[2] * v;
294 u[2]->add (-gamma / std::sqrt(delta[1]), *u[1]);
300 precondition.vmult(v,*u[2]);
302 delta[2] = v * (*u[2]);
309 e[1] = std::sqrt(delta[2]);
313 d_ = s * e[0] - c * gamma;
314 e[0] = c * e[0] + s * gamma;
315 f[1] = s * std::sqrt(delta[2]);
316 e[1] = -c * std::sqrt(delta[2]);
319 const double d = std::sqrt (d_*d_ + delta[2]);
326 s = std::sqrt(delta[2]) / d;
331 m[0]->add (-e[0], *m[1]);
333 m[0]->add (-f[0], *m[2]);
336 r_l2 *= std::fabs(s);
389 DEAL_II_NAMESPACE_CLOSE
SolverMinRes(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
#define AssertThrow(cond, exc)
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
Stop iteration, goal reached.
#define Assert(cond, exc)
#define DeclException0(Exception0)
VectorMemory< VectorType > & memory
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType ¤t_iterate), StateCombiner > iteration_status
void push(const std::string &text)
virtual double criterion()
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &precondition)
static::ExceptionBase & ExcPreconditionerNotDefinite()