#include <deal.II/base/config.h>
-#include <deal.II/lac/shifted_matrix.h>
-#include <deal.II/lac/solver.h>
-#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/vector_memory.h>
-#include <deal.II/lac/precondition.h>
#include <cmath>
SolverControl::State conv=SolverControl::iterate;
// Prepare matrix for solver
- ShiftedMatrix <MatrixType> A_s(A, -value);
+ auto A_op = linear_operator(A);
+ double current_shift = -value;
+ auto A_s = A_op + current_shift * identity_operator(A_op);
// Define solver
ReductionControl inner_control (5000, 1.e-16, 1.e-5, false, false);
// Compute unshifted eigenvalue
value = (entry * x (i) < 0.) ? -length : length;
value = 1./value;
- value -= A_s.shift ();
+ value -= current_shift;
if (iter==goal)
{
- const double new_shift = - additional_data.relaxation * value
- + (1.-additional_data.relaxation) * A_s.shift();
- A_s.shift(new_shift);
+ const auto &relaxation = additional_data.relaxation;
+ const double new_shift =
+ relaxation * (-value) + (1. - relaxation) * current_shift;
+
+ A_s = A_op + new_shift * identity_operator(A_op);
+ current_shift = new_shift;
+
++goal;
}