15#ifndef dealii_solver_idr_h
16#define dealii_solver_idr_h
46 namespace SolverIDRImplementation
52 template <
typename VectorType>
91 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
118template <
typename VectorType = Vector<
double>>
135 const unsigned int s;
160 template <
typename MatrixType,
typename PreconditionerType>
167 const PreconditionerType &preconditioner);
176 print_vectors(
const unsigned int step,
196 namespace SolverIDRImplementation
198 template <
typename VectorType>
207 template <
typename VectorType>
219 template <
typename VectorType>
222 const VectorType &
temp)
225 if (
data[i] ==
nullptr)
235 template <
typename VectorType,
236 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
239 n_blocks(
const VectorType &)
246 template <
typename VectorType,
247 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
250 n_blocks(
const VectorType &vector)
252 return vector.n_blocks();
257 template <
typename VectorType,
258 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
261 block(VectorType &vector,
const unsigned int b)
270 template <
typename VectorType,
271 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
273 typename VectorType::BlockType &
274 block(VectorType &vector,
const unsigned int b)
276 return vector.block(b);
284template <
typename VectorType>
288 const AdditionalData &
data)
290 , additional_data(
data)
295template <
typename VectorType>
299 , additional_data(
data)
304template <
typename VectorType>
309 const VectorType &)
const
314template <
typename VectorType>
316template <
typename MatrixType,
typename PreconditionerType>
323 const PreconditionerType &preconditioner)
328 unsigned int step = 0;
330 const unsigned int s = additional_data.
s;
343 uhat.reinit(
x,
true);
347 r.sadd(-1.0, 1.0, b);
349 using value_type =
typename VectorType::value_type;
353 real_type
res =
r.l2_norm();
368 for (
unsigned int i = 0; i < s; ++i)
382 for (
unsigned int b = 0;
383 b < internal::SolverIDRImplementation::n_blocks(
tmp_q);
386 .locally_owned_elements())
394 for (
unsigned int j = 0;
j < i; ++
j)
418 for (
unsigned int i = 0; i < s; ++i)
422 for (
unsigned int k = 0;
k < s; ++
k)
429 std::vector<unsigned int> indices;
431 for (
unsigned int i =
k; i < s; ++i, ++
j)
433 indices.push_back(i);
436 Mk.extract_submatrix_from(
M, indices, indices);
447 for (
unsigned int i =
k,
j = 0; i < s; ++i, ++
j)
451 preconditioner.vmult(
uhat, v);
456 for (
unsigned int i =
k + 1,
j = 1; i < s; ++i, ++
j)
469 for (
unsigned int i = 1; i <
k; ++i)
488 for (
unsigned int i =
k + 1; i < s; ++i)
489 M(i,
k) =
Q[i] * G[
k];
497 print_vectors(step,
x,
r, U[
k]);
513 for (
unsigned int i = 0; i <
k + 1; ++i)
515 for (
unsigned int i =
k + 1; i < s; ++i)
524 preconditioner.vmult(
uhat,
r);
527 omega = (v *
r) / (v * v);
532 print_vectors(step,
x,
r,
uhat);
@ 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())
virtual ~SolverIDR() override=default
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
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_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
std::vector< index_type > data
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
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)