15#ifndef dealii_solver_idr_h
16#define dealii_solver_idr_h
45 namespace SolverIDRImplementation
51 template <
typename VectorType>
79 operator()(
const unsigned int i,
const VectorType &temp);
90 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
117template <
typename VectorType = Vector<
double>>
133 const unsigned int s;
158 template <
typename MatrixType,
typename PreconditionerType>
163 const PreconditionerType &preconditioner);
175 const VectorType &d)
const;
192 namespace SolverIDRImplementation
194 template <
typename VectorType>
203 template <
typename VectorType>
215 template <
typename VectorType>
218 const VectorType &temp)
221 if (data[i] ==
nullptr)
224 data[i]->reinit(temp,
true);
231 template <
typename VectorType,
232 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
235 n_blocks(
const VectorType &)
242 template <
typename VectorType,
243 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
246 n_blocks(
const VectorType &vector)
248 return vector.n_blocks();
253 template <
typename VectorType,
254 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
257 block(VectorType &vector,
const unsigned int b)
266 template <
typename VectorType,
267 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
269 typename VectorType::BlockType &
270 block(VectorType &vector,
const unsigned int b)
272 return vector.block(b);
280template <
typename VectorType>
283 const AdditionalData &data)
285 , additional_data(data)
290template <
typename VectorType>
293 , additional_data(data)
298template <
typename VectorType>
303 const VectorType &)
const
308template <
typename VectorType>
309template <
typename MatrixType,
typename PreconditionerType>
314 const PreconditionerType &preconditioner)
319 unsigned int step = 0;
321 const unsigned int s = additional_data.s;
328 VectorType &r = *r_pointer;
329 VectorType &v = *v_pointer;
330 VectorType &uhat = *uhat_pointer;
334 uhat.reinit(x,
true);
338 r.sadd(-1.0, 1.0, b);
340 using value_type =
typename VectorType::value_type;
344 real_type res = r.l2_norm();
345 iteration_state = this->iteration_status(step, res, x);
358 std::normal_distribution<> normal_distribution(0.0, 1.0);
359 for (
unsigned int i = 0; i < s; ++i)
370 VectorType &tmp_q = Q(i, x);
373 for (
unsigned int b = 0;
374 b < internal::SolverIDRImplementation::n_blocks(tmp_q);
376 for (
auto indx :
internal::SolverIDRImplementation::block(tmp_q,
b)
377 .locally_owned_elements())
378 internal::SolverIDRImplementation::block(tmp_q,
b)(indx) =
379 normal_distribution(rng);
385 for (
unsigned int j = 0; j < i; ++j)
388 v *= (v * tmp_q) / (tmp_q * tmp_q);
393 tmp_q *= 1.0 / tmp_q.l2_norm();
398 value_type omega = 1.;
400 bool early_exit =
false;
409 for (
unsigned int i = 0; i < s; ++i)
413 for (
unsigned int k = 0; k < s; ++k)
420 std::vector<unsigned int> indices;
422 for (
unsigned int i = k; i < s; ++i, ++j)
424 indices.push_back(i);
427 Mk.extract_submatrix_from(M, indices, indices);
431 Mk_inv.vmult(gamma, phik);
438 for (
unsigned int i = k, j = 0; i < s; ++i, ++j)
439 v.add(-
gamma(j), G[i]);
442 preconditioner.vmult(uhat, v);
446 uhat.sadd(omega,
gamma(0), U[k]);
447 for (
unsigned int i = k + 1, j = 1; i < s; ++i, ++j)
448 uhat.add(
gamma(j), U[i]);
459 value_type alpha = Q[0] * G[k] / M(0, 0);
460 for (
unsigned int i = 1; i < k; ++i)
462 const value_type alpha_old = alpha;
463 alpha = G[k].add_and_dot(-alpha, G[i - 1], Q[i]) / M(i, i);
467 uhat.add(-alpha_old, U[i - 1], -alpha, U[i]);
469 M(k, k) = G[k].add_and_dot(-alpha, G[k - 1], Q[k]);
471 uhat.add(-alpha, U[k - 1]);
474 M(k, k) = G[k] * Q[k];
479 for (
unsigned int i = k + 1; i < s; ++i)
480 M(i, k) = Q[i] * G[k];
484 const value_type beta = phi(k) / M(k, k);
488 print_vectors(step, x, r, U[k]);
494 iteration_state = this->iteration_status(step, res, x);
504 for (
unsigned int i = 0; i < k + 1; ++i)
506 for (
unsigned int i = k + 1; i < s; ++i)
507 phi(i) -= beta * M(i, k);
511 if (early_exit ==
true)
515 preconditioner.vmult(uhat, r);
518 omega = (v * r) / (v * v);
520 res =
std::sqrt(r.add_and_dot(-1.0 * omega, v, r));
523 print_vectors(step, x, r, uhat);
526 iteration_state = this->iteration_status(step, res, x);
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
SolverIDR(SolverControl &cn, const AdditionalData &data=AdditionalData())
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
SolverIDR(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
virtual ~SolverIDR() override=default
TmpVectors(const unsigned int s_param, VectorMemory< VectorType > &vmem)
VectorType & operator()(const unsigned int i, const VectorType &temp)
VectorType & operator[](const unsigned int i) const
VectorMemory< VectorType > & mem
std::vector< typename VectorMemory< VectorType >::Pointer > data
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
long double gamma(const unsigned int n)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int s=2)