16 #ifndef dealii_solver_idr_h
17 #define dealii_solver_idr_h
46 namespace SolverIDRImplementation
52 template <
typename VectorType>
80 operator()(
const unsigned int i,
const VectorType &temp);
91 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
118 template <
typename VectorType = Vector<
double>>
134 const unsigned int s;
159 template <
typename MatrixType,
typename PreconditionerType>
164 const PreconditionerType &preconditioner);
176 const VectorType &
d)
const;
193 namespace SolverIDRImplementation
195 template <
typename VectorType>
204 template <
typename VectorType>
216 template <
typename VectorType>
219 const VectorType &temp)
222 if (data[i] ==
nullptr)
225 data[i]->reinit(temp,
true);
232 template <
typename VectorType,
233 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
243 template <
typename VectorType,
244 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
249 return vector.n_blocks();
254 template <
typename VectorType,
255 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
258 block(VectorType &vector,
const unsigned int b)
267 template <
typename VectorType,
268 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
270 typename VectorType::BlockType &
271 block(VectorType &vector,
const unsigned int b)
273 return vector.block(
b);
281 template <
typename VectorType>
284 const AdditionalData &data)
286 , additional_data(data)
291 template <
typename VectorType>
294 , additional_data(data)
299 template <
typename VectorType>
304 const VectorType &)
const
309 template <
typename VectorType>
310 template <
typename MatrixType,
typename PreconditionerType>
315 const PreconditionerType &preconditioner)
320 unsigned int step = 0;
322 const unsigned int s = additional_data.s;
329 VectorType &r = *r_pointer;
330 VectorType &v = *v_pointer;
331 VectorType &uhat = *uhat_pointer;
335 uhat.reinit(x,
true);
339 r.sadd(-1.0, 1.0,
b);
341 using value_type =
typename VectorType::value_type;
345 real_type res = r.l2_norm();
346 iteration_state = this->iteration_status(step, res, x);
359 std::normal_distribution<> normal_distribution(0.0, 1.0);
360 for (
unsigned int i = 0; i < s; ++i)
371 VectorType &tmp_q = Q(i, x);
374 for (
unsigned int b = 0;
377 for (
auto indx : internal::SolverIDRImplementation::block(tmp_q,
b)
378 .locally_owned_elements())
379 internal::SolverIDRImplementation::block(tmp_q,
b)(indx) =
380 normal_distribution(rng);
386 for (
unsigned int j = 0; j < i; ++j)
389 v *= (v * tmp_q) / (tmp_q * tmp_q);
394 tmp_q *= 1.0 / tmp_q.l2_norm();
399 value_type omega = 1.;
401 bool early_exit =
false;
410 for (
unsigned int i = 0; i < s; ++i)
414 for (
unsigned int k = 0; k < s; ++k)
421 std::vector<unsigned int> indices;
423 for (
unsigned int i = k; i < s; ++i, ++j)
425 indices.push_back(i);
428 Mk.extract_submatrix_from(M, indices, indices);
432 Mk_inv.vmult(
gamma, phik);
439 for (
unsigned int i = k, j = 0; i < s; ++i, ++j)
440 v.add(-
gamma(j), G[i]);
443 preconditioner.vmult(uhat, v);
447 uhat.sadd(omega,
gamma(0),
U[k]);
448 for (
unsigned int i = k + 1, j = 1; i < s; ++i, ++j)
460 value_type alpha = Q[0] * G[k] / M(0, 0);
461 for (
unsigned int i = 1; i < k; ++i)
463 const value_type alpha_old = alpha;
464 alpha = G[k].add_and_dot(-alpha, G[i - 1], Q[i]) / M(i, i);
468 uhat.add(-alpha_old,
U[i - 1], -alpha,
U[i]);
470 M(k, k) = G[k].add_and_dot(-alpha, G[k - 1], Q[k]);
472 uhat.add(-alpha,
U[k - 1]);
475 M(k, k) = G[k] * Q[k];
480 for (
unsigned int i = k + 1; i < s; ++i)
481 M(i, k) = Q[i] * G[k];
485 const value_type beta = phi(k) / M(k, k);
489 print_vectors(step, x, r,
U[k]);
495 iteration_state = this->iteration_status(step, res, x);
505 for (
unsigned int i = 0; i < k + 1; ++i)
507 for (
unsigned int i = k + 1; i < s; ++i)
508 phi(i) -= beta * M(i, k);
512 if (early_exit ==
true)
516 preconditioner.vmult(uhat, r);
519 omega = (v * r) / (v * v);
521 res =
std::sqrt(r.add_and_dot(-1.0 * omega, v, r));
524 print_vectors(step, x, r, uhat);
527 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)
VectorMemory< VectorType > & mem
VectorType & operator[](const unsigned int i) const
std::vector< typename VectorMemory< VectorType >::Pointer > data
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)