16 #ifndef dealii_eigen_h 17 #define dealii_eigen_h 20 #include <deal.II/base/config.h> 22 #include <deal.II/lac/linear_operator.h> 23 #include <deal.II/lac/precondition.h> 24 #include <deal.II/lac/solver.h> 25 #include <deal.II/lac/solver_cg.h> 26 #include <deal.II/lac/solver_control.h> 27 #include <deal.II/lac/solver_gmres.h> 28 #include <deal.II/lac/solver_minres.h> 29 #include <deal.II/lac/vector_memory.h> 33 DEAL_II_NAMESPACE_OPEN
54 template <
typename VectorType = Vector<
double>>
99 template <
typename MatrixType>
101 solve(
double &value,
const MatrixType &A, VectorType &x);
133 template <
typename VectorType = Vector<
double>>
164 unsigned int start_adaption = 6,
165 bool use_residual =
true)
166 : relaxation(relaxation)
167 , start_adaption(start_adaption)
168 , use_residual(use_residual)
192 template <
typename MatrixType>
194 solve(
double &value,
const MatrixType &A, VectorType &x);
207 template <
class VectorType>
211 :
Solver<VectorType>(cn, mem)
217 template <
class VectorType>
223 template <
class VectorType>
224 template <
typename MatrixType>
233 VectorType & y = *Vy;
236 VectorType & r = *Vr;
239 double length = x.l2_norm();
240 double old_length = 0.;
253 length = y.l2_norm();
259 double thresh = length / x.size();
265 while (std::fabs(entry) < thresh);
270 value = (entry * x(i) < 0.) ? -length : length;
274 x.equ(1 / length, y);
282 std::fabs(1. / length - 1. / old_length),
289 iter, std::fabs(1. / length - 1. / old_length)));
296 template <
class VectorType>
300 :
Solver<VectorType>(cn, mem)
306 template <
class VectorType>
312 template <
class VectorType>
313 template <
typename MatrixType>
325 double current_shift = -value;
338 VectorType & y = *Vy;
341 VectorType & r = *Vr;
344 double length = x.l2_norm();
345 double old_value = value;
350 double res = -std::numeric_limits<double>::max();
354 solver.
solve(A_s, y, x, prec);
357 length = y.l2_norm();
363 double thresh = length / x.size();
369 while (std::fabs(entry) < thresh);
374 value = (entry * x(i) < 0. ? -1. : 1.) / length - current_shift;
379 const double new_shift =
380 relaxation * (-value) + (1. - relaxation) * current_shift;
383 current_shift = new_shift;
389 x.equ(1. / length, y);
395 r.sadd(-1., value, x);
402 res = std::fabs(1. / value - 1. / old_value);
415 DEAL_II_NAMESPACE_CLOSE
types::global_dof_index size_type
AdditionalData additional_data
void solve(double &value, const MatrixType &A, VectorType &x)
#define AssertThrow(cond, exc)
unsigned long long int global_dof_index
void solve(double &value, const MatrixType &A, VectorType &x)
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
Stop iteration, goal reached.
#define Assert(cond, exc)
AdditionalData additional_data
EigenPower(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType ¤t_iterate), StateCombiner > iteration_status
AdditionalData(const double shift=0.)
EigenInverse(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
VectorMemory< VectorType > & memory
unsigned int start_adaption
LinearOperator< Range, Range, Payload > identity_operator(const std::function< void(Range &, bool)> &reinit_vector)
types::global_dof_index size_type
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
AdditionalData(double relaxation=1., unsigned int start_adaption=6, bool use_residual=true)
static::ExceptionBase & ExcInternalError()