17 #ifndef __deal2__eigen_h
18 #define __deal2__eigen_h
21 #include <deal.II/base/config.h>
22 #include <deal.II/lac/shifted_matrix.h>
23 #include <deal.II/lac/solver.h>
24 #include <deal.II/lac/solver_control.h>
25 #include <deal.II/lac/solver_cg.h>
26 #include <deal.II/lac/solver_gmres.h>
27 #include <deal.II/lac/solver_minres.h>
28 #include <deal.II/lac/vector_memory.h>
29 #include <deal.II/lac/precondition.h>
54 template <
class VECTOR = Vector<
double> >
92 const AdditionalData &data=AdditionalData());
110 template <
class MATRIX>
112 solve (
double &value,
148 template <
class VECTOR = Vector<
double> >
182 unsigned int start_adaption = 6,
183 bool use_residual =
true)
185 relaxation(relaxation),
186 start_adaption(start_adaption),
187 use_residual(use_residual)
197 const AdditionalData &data=AdditionalData());
218 template <
class MATRIX>
220 solve (
double &value,
235 template <
class VECTOR>
241 additional_data(data)
246 template <
class VECTOR>
252 template <
class VECTOR>
253 template <
class MATRIX>
261 deallog.
push(
"Power method");
263 VECTOR *Vy = this->memory.alloc ();
266 VECTOR *Vr = this->memory.alloc ();
270 double length = x.l2_norm ();
271 double old_length = 0.;
279 y.add(additional_data.shift, x);
283 length = y.l2_norm ();
289 double thresh = length/x.size();
295 while (std::fabs(entry) < thresh);
300 value = (entry * x (i) < 0.) ? -length : length;
301 value -= additional_data.shift;
311 conv = this->control().check (iter, std::fabs(1./length-1./old_length));
314 this->memory.free(Vy);
315 this->memory.free(Vr);
322 this->control().last_value()));
328 template <
class VECTOR>
334 additional_data(data)
339 template <
class VECTOR>
345 template <
class VECTOR>
346 template <
class MATRIX>
352 deallog.
push(
"Wielandt");
363 solver(inner_control, this->memory);
366 unsigned int goal = additional_data.start_adaption;
369 VECTOR *Vy = this->memory.alloc ();
372 VECTOR *Vr = this->memory.alloc ();
376 double length = x.l2_norm ();
377 double old_value = value;
384 solver.
solve (A_s, y, x, prec);
387 length = y.l2_norm ();
393 double thresh = length/x.size();
399 while (std::fabs(entry) < thresh);
404 value = (entry * x (i) < 0.) ? -length : length;
406 value -= A_s.
shift ();
410 const double new_shift = - additional_data.relaxation * value
411 + (1.-additional_data.relaxation) * A_s.
shift();
412 A_s.
shift(new_shift);
417 x.equ (1./length, y);
419 if (additional_data.use_residual)
423 r.sadd(-1., value, x);
424 double res = r.l2_norm();
426 conv = this->control().check (iter, res);
430 conv = this->control().check (iter, std::fabs(1./value-1./old_value));
435 this->memory.free(Vy);
436 this->memory.free(Vr);
444 this->control().last_value());
448 DEAL_II_NAMESPACE_CLOSE
EigenInverse(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
void solve(double &value, const MATRIX &A, VECTOR &x)
void vmult(VECTOR &u, const VECTOR &v) const
types::global_dof_index size_type
void shift(const double sigma)
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
#define AssertThrow(cond, exc)
AdditionalData additional_data
unsigned int global_dof_index
#define Assert(cond, exc)
AdditionalData additional_data
AdditionalData(const double shift=0.)
void solve(double &value, const MATRIX &A, VECTOR &x)
unsigned int start_adaption
void push(const std::string &text)
types::global_dof_index size_type
AdditionalData(double relaxation=1., unsigned int start_adaption=6, bool use_residual=true)
::ExceptionBase & ExcInternalError()
EigenPower(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
Stop iteration, goal reached.