16 #ifndef dealii_solver_idr_h
17 #define dealii_solver_idr_h
44 namespace SolverIDRImplementation
50 template <
typename VectorType>
88 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
117 template <
class VectorType = Vector<
double>>
133 const unsigned int s;
158 template <
typename MatrixType,
typename PreconditionerType>
160 solve(
const MatrixType &
A,
163 const PreconditionerType &preconditioner);
192 namespace SolverIDRImplementation
194 template <
class VectorType>
203 template <
class VectorType>
215 template <
class VectorType>
221 if (data[i] ==
nullptr)
224 data[i]->reinit(temp);
233 template <
class VectorType>
236 const AdditionalData & data)
238 , additional_data(data)
243 template <
class VectorType>
246 , additional_data(data)
251 template <
class VectorType>
261 template <
class VectorType>
262 template <
typename MatrixType,
typename PreconditionerType>
267 const PreconditionerType &preconditioner)
272 unsigned int step = 0;
274 const unsigned int s = additional_data.s;
291 vhat.reinit(x,
true);
292 uhat.reinit(x,
true);
293 ghat.reinit(x,
true);
297 r.sadd(-1.0, 1.0,
b);
300 double res = r.l2_norm();
301 iteration_state = this->iteration_status(step, res, x);
314 std::normal_distribution<> normal_distribution(0.0, 1.0);
315 for (
unsigned int i = 0; i < s; ++i)
330 for (
auto indx : tmp_q.locally_owned_elements())
331 tmp_q(indx) = normal_distribution(rng);
337 for (
unsigned int j = 0; j < i; ++j)
340 v *= (v * tmp_q) / (tmp_q * tmp_q);
345 tmp_q *= 1.0 / tmp_q.l2_norm();
352 bool early_exit =
false;
361 for (
unsigned int i = 0; i < s; ++i)
365 for (
unsigned int k = 0; k < s; ++k)
372 std::vector<unsigned int> indices;
374 for (
unsigned int i = k; i < s; ++i, ++j)
376 indices.push_back(i);
379 Mk.extract_submatrix_from(M, indices, indices);
383 Mk_inv.vmult(
gamma, phik);
390 for (
unsigned int i = k; i < s; ++i, ++j)
391 v.add(-1.0 *
gamma(j), G[i]);
392 preconditioner.vmult(vhat, v);
397 for (
unsigned int i = k; i < s; ++i, ++j)
405 for (
unsigned int i = 0; i < k; ++i)
407 double alpha = (Q[i] * ghat) / M(i, i);
408 ghat.add(-alpha, G[i]);
409 uhat.add(-alpha,
U[i]);
415 for (
unsigned int i = k; i < s; ++i)
416 M(i, k) = Q[i] * G[k];
421 double beta = phi(k) / M(k, k);
422 r.add(-1.0 * beta, G[k]);
425 print_vectors(step, x, r,
U[k]);
431 iteration_state = this->iteration_status(step, res, x);
441 for (
unsigned int i = 0; i < k + 1; ++i)
443 for (
unsigned int i = k + 1; i < s; ++i)
444 phi(i) -= beta * M(i, k);
448 if (early_exit ==
true)
452 preconditioner.vmult(vhat, r);
455 omega = (v * r) / (v * v);
457 r.add(-1.0 * omega, v);
460 print_vectors(step, x, r, vhat);
464 iteration_state = this->iteration_status(step, res, x);