│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
214 template<
class M,
class RI>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
224 condition_estimate_(condition_estimate)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
227 condition_estimate_ =
false;
│ │ │ │ +
228 std::cerr <<
"WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
243 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
266 const std::vector<typename M::size_type>& rowsize_)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
279 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
285 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
288 if(!OwnerSet::contains(i->local().attribute())) {
│ │ │ │ -
│ │ │ │ -
290 std::cout<<rank<<
" Inserting diagonal for"<<i->local()<<std::endl;
│ │ │ │ -
│ │ │ │ -
292 sparsity[i->local()].insert(i->local());
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
240 scalar_real_type reduction,
int maxit,
int verbose,
bool condition_estimate) :
IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
│ │ │ │ +
241 condition_estimate_(condition_estimate)
│ │ │ │ +
│ │ │ │ +
243 if (condition_estimate && !(std::is_same<field_type,float>::value || std::is_same<field_type,double>::value)) {
│ │ │ │ +
244 condition_estimate_ =
false;
│ │ │ │ +
245 std::cerr <<
"WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
258 scalar_real_type reduction,
int maxit,
int verbose,
bool condition_estimate)
│ │ │ │ +
│ │ │ │ +
260 condition_estimate_(condition_estimate)
│ │ │ │ +
│ │ │ │ +
262 if (condition_estimate && !(std::is_same<field_type,float>::value || std::is_same<field_type,double>::value)) {
│ │ │ │ +
263 condition_estimate_ =
false;
│ │ │ │ +
264 std::cerr <<
"WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
284 _op->applyscaleadd(-1,x,b);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
287 if(iteration.step(0, def)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
296 std::vector<real_type> lambdas(0);
│ │ │ │ +
297 std::vector<real_type> betas(0);
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
301 m.setBuildMode(M::row_wise);
│ │ │ │ -
302 typename M::CreateIterator citer=m.createbegin();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
306 Dune::GlobalLookupIndexSet<I> global(
aggidxset);
│ │ │ │ -
│ │ │ │ -
308 typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
311 typedef typename std::set<size_type>::const_iterator SIter;
│ │ │ │ -
312 for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
315 if(i->find(idx)==i->end()) {
│ │ │ │ -
316 const typename I::IndexPair* gi=global.pair(idx);
│ │ │ │ -
│ │ │ │ -
318 std::cout<<rank<<
": row "<<idx<<
" is missing a diagonal entry! global="<<gi->global()<<
" attr="<<gi->local().attribute()<<
" "<<
│ │ │ │ -
319 OwnerSet::contains(gi->local().attribute())<<
│ │ │ │ -
320 " row size="<<i->size()<<std::endl;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
342 for (
unsigned int i = 0; i !=
sparsity.size(); ++i) {
│ │ │ │ -
343 if (add_sparsity[i].size() != 0) {
│ │ │ │ -
344 typedef std::set<size_type> Set;
│ │ │ │ -
│ │ │ │ -
346 std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin());
│ │ │ │ -
347 std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(),
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
356 const Dune::GlobalLookupIndexSet<I>&
idxset;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
362 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
376 static typename M::size_type
getSize(
const Type& t, std::size_t i)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
379 return t.
matrix[i].size();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
394 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
405 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
412 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_,
│ │ │ │ -
413 std::vector<size_t>& rowsize_)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
423 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
305 rholast =
_sp->dot(p,b);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
313 alpha =
_sp->dot(p,q);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
316 if (condition_estimate_)
│ │ │ │ +
317 lambdas.push_back(std::real(lambda));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
323 if(iteration.step(i, def))
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
332 if (condition_estimate_)
│ │ │ │ +
333 betas.push_back(std::real(beta));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
341 if (condition_estimate_) {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
354 row.insert(row.index()-1);
│ │ │ │ +
355 row.insert(row.index());
│ │ │ │ +
356 if (row.index() < T.
N() - 1)
│ │ │ │ +
357 row.insert(row.index()+1);
│ │ │ │ +
│ │ │ │ +
359 for (
int row = 0; row < i; ++row) {
│ │ │ │ +
│ │ │ │ +
361 T[row][row-1] = sqrt(betas[row-1]) / lambdas[row-1];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
364 T[row][row] = 1.0 / lambdas[row];
│ │ │ │ +
│ │ │ │ +
366 T[row][row] += betas[row-1] / lambdas[row-1];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
370 T[row][row+1] = sqrt(betas[row]) / lambdas[row];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
386 std::cout <<
"Min eigv estimate: " << Simd::io(min_eigv) <<
'\n';
│ │ │ │ +
387 std::cout <<
"Max eigv estimate: " << Simd::io(max_eigv) <<
'\n';
│ │ │ │ +
388 std::cout <<
"Condition estimate: "
│ │ │ │ +
389 << Simd::io(max_eigv / min_eigv) << std::endl;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
393 std::cerr <<
"WARNING: Condition estimate was requested. This requires ARPACK, but ARPACK was not found!" << std::endl;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
399 bool condition_estimate_ =
false;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
427 if(!OwnerSet::contains(i->local().attribute())) {
│ │ │ │ -
│ │ │ │ -
429 typedef typename M::ColIterator CIter;
│ │ │ │ -
430 for(CIter c=
matrix[i->local()].begin(), cend=
matrix[i->local()].end();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
434 if(c.index()==i->local()) {
│ │ │ │ -
435 auto setDiagonal = [](
auto&& scalarOrMatrix,
const auto& value) {
│ │ │ │ -
436 auto&& matrixView = Dune::Impl::asMatrix(scalarOrMatrix);
│ │ │ │ -
437 for (
auto rowIt = matrixView.begin(); rowIt != matrixView.end(); ++rowIt)
│ │ │ │ -
438 (*rowIt)[rowIt.index()] = value;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
448 const Dune::GlobalLookupIndexSet<I>&
idxset;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
455 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
464 typedef std::pair<typename I::GlobalIndex,typename M::block_type>
IndexedType;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
442 const Simd::Scalar<real_type> EPSILON=1e-80;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
445 field_type rho, rho_new, alpha, beta, h, omega;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
472 return t.
matrix[i].size();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
481 template<
class M,
class I,
class RI>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
466 _op->applyscaleadd(-1,x,r);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
471 if(iteration.step(0, norm)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
488 return cont.
matrix[i].size();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
493 cont.
rowsize.getRowSize(i)=rowsize;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
498 template<
class M,
class I,
class RI>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
505 return cont.
matrix[i].size();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
510 if (rowsize > cont.
rowsize.getCopyRowSize(i))
│ │ │ │ -
511 cont.
rowsize.getCopyRowSize(i)=rowsize;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
516 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
486 for (it = 0.5; it <
_maxit; it+=.5)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
493 rho_new =
_sp->dot(rt,r);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
496 if (Simd::allTrue(abs(rho) <= EPSILON))
│ │ │ │ +
497 DUNE_THROW(
SolverAbort,
"breakdown in BiCGSTAB - rho "
│ │ │ │ +
498 << Simd::io(rho) <<
" <= EPSILON " << EPSILON
│ │ │ │ +
499 <<
" after " << it <<
" iterations");
│ │ │ │ +
500 if (Simd::allTrue(abs(omega) <= EPSILON))
│ │ │ │ +
501 DUNE_THROW(
SolverAbort,
"breakdown in BiCGSTAB - omega "
│ │ │ │ +
502 << Simd::io(omega) <<
" <= EPSILON " << EPSILON
│ │ │ │ +
503 <<
" after " << it <<
" iterations");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
512 ( rho_new / rho ) * ( alpha / omega ));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
528 if ( Simd::allTrue(abs(h) < EPSILON) )
│ │ │ │ +
529 DUNE_THROW(
SolverAbort,
"abs(h) < EPSILON in BiCGSTAB - abs(h) "
│ │ │ │ +
530 << Simd::io(abs(h)) <<
" < EPSILON " << EPSILON
│ │ │ │ +
531 <<
" after " << it <<
" iterations");
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
538 numlimits = std::numeric_limits<GlobalIndex>::max();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
542 const typename I::IndexPair* index=cont.
idxset.pair(
col.index());
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
545 if ( index->local().attribute() != 2)
│ │ │ │ -
546 return index->global();
│ │ │ │ -
│ │ │ │ -
548 numlimits = std::numeric_limits<GlobalIndex>::max();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
556 if (gi != std::numeric_limits<GlobalIndex>::max()) {
│ │ │ │ -
557 const typename I::IndexPair& ip=cont.
aggidxset.at(gi);
│ │ │ │ -
558 assert(ip.global()==gi);
│ │ │ │ -
559 std::size_t column = ip.local();
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
549 if(iteration.step(it, norm)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
563 if(!OwnerSet::contains(ip.local().attribute()))
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
568 catch(
const Dune::RangeError&) {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
572 typedef typename GlobalLookup::IndexPair IndexPair;
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
576 const IndexPair* pi=lookup.pair(i);
│ │ │ │ -
│ │ │ │ -
578 if(OwnerSet::contains(pi->local().attribute())) {
│ │ │ │ -
│ │ │ │ -
580 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
│ │ │ │ -
581 std::cout<<rank<<cont.
aggidxset<<std::endl;
│ │ │ │ -
582 std::cout<<rank<<
": row "<<i<<
" (global="<<gi <<
") not in index set for owner index "<<pi->global()<<std::endl;
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
582 if(iteration.step(it, norm)){
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
590 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
593 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
597 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
603 typedef typename std::pair<GlobalIndex,typename M::block_type>
Data;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
619 numlimits = std::numeric_limits<GlobalIndex>::max();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
625 const typename I::IndexPair* index=cont.
idxset.pair(
col.index());
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
629 if ( index->local().attribute() != 2)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
632 numlimits = std::numeric_limits<GlobalIndex>::max();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
641 if (data.first != std::numeric_limits<GlobalIndex>::max()) {
│ │ │ │ -
642 typename M::size_type column=cont.
aggidxset.at(data.first).local();
│ │ │ │ -
643 cont.
matrix[i][column]=data.second;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
646 catch(
const Dune::RangeError&) {
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
597 template<
class CountType>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
636 _op->applyscaleadd(-1.0,x,b);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
645 if (iteration.step(0, def)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
653 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
656 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
659 template<
class M,
class I>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
662 template<
typename M,
typename C>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
666 typename C::CopySet copyflags;
│ │ │ │ -
667 typename C::OwnerSet ownerflags;
│ │ │ │ -
668 typedef typename C::ParallelIndexSet IndexSet;
│ │ │ │ -
│ │ │ │ -
670 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
│ │ │ │ -
671 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
│ │ │ │ -
672 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
653 std::array<real_type,2> c{{0.0,0.0}};
│ │ │ │ +
│ │ │ │ +
655 std::array<field_type,2> s{{0.0,0.0}};
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
658 std::array<field_type,3> T{{0.0,0.0,0.0}};
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
661 std::array<field_type,2> xi{{1.0,0.0}};
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
665 beta = sqrt(
_sp->dot(b,z));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
669 std::array<X,3> p{{b,b,b}};
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
676 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
│ │ │ │ -
│ │ │ │ -
678 origComm.buildGlobalLookup();
│ │ │ │ -
│ │ │ │ -
680 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
685 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
│ │ │ │ -
│ │ │ │ -
687 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
│ │ │ │ -
│ │ │ │ -
689 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
693 RemoteIndices<IndexSet> *ris =
new RemoteIndices<IndexSet>(origComm.indexSet(),
│ │ │ │ -
│ │ │ │ -
695 origComm.communicator());
│ │ │ │ -
696 ris->template rebuild<true>();
│ │ │ │ -
│ │ │ │ -
698 ri.getInterface().free();
│ │ │ │ -
699 ri.getInterface().build(*ris,copyflags,ownerflags);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
703 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
706 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
711 for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
717 origComm.globalLookup(),
│ │ │ │ -
│ │ │ │ -
719 backwardscopyrowsize);
│ │ │ │ -
│ │ │ │ -
721 newComm.indexSet(), copyrowsize);
│ │ │ │ -
722 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
731#ifdef DUNE_ISTL_WITH_CHECKING
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
734 typedef typename M::ConstRowIterator RIter;
│ │ │ │ -
735 for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row) {
│ │ │ │ -
736 typedef typename M::ConstColIterator CIter;
│ │ │ │ -
737 for(CIter
col=row->begin(), cend=row->end();
col!=cend; ++
col)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
740 newMatrix[
col.index()][row.index()];
│ │ │ │ -
│ │ │ │ -
742 std::cerr<<newComm.communicator().rank()<<
": entry ("
│ │ │ │ -
743 <<
col.index()<<
","<<row.index()<<
") missing! for symmetry!"<<std::endl;
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
675 std::array<X,3> q{{b,b,b}};
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
697 q[i2].axpy(-beta,q[i0]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
701 alpha =
_sp->dot(z,q[i2]);
│ │ │ │ +
702 q[i2].axpy(-alpha,q[i1]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
705 _prec->apply(z,q[i2]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
709 beta = sqrt(
_sp->dot(q[i2],z));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
726 T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
│ │ │ │ +
727 T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
733 generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
│ │ │ │ +
│ │ │ │ +
735 T[2] = c[i%2]*T[2] + s[i%2]*beta;
│ │ │ │ +
│ │ │ │ +
737 xi[i%2] = -s[i%2]*xi[(i+1)%2];
│ │ │ │ +
738 xi[(i+1)%2] *= c[i%2];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
742 p[i2].axpy(-T[1],p[i1]);
│ │ │ │ +
743 p[i2].axpy(-T[0],p[i0]);
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
752 DUNE_THROW(
ISTLError,
"Matrix not symmetric!");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
756 template<
typename M,
typename C>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
760 typedef typename C::ParallelIndexSet IndexSet;
│ │ │ │ -
761 typename C::OwnerSet ownerflags;
│ │ │ │ -
762 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
│ │ │ │ -
763 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
│ │ │ │ -
764 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
│ │ │ │ +
│ │ │ │ +
747 x.axpy(beta0*xi[(i+1)%2],p[i2]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
754 def = abs(beta0*xi[i%2]);
│ │ │ │ +
755 if(iteration.step(i, def)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
766 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
773 for (std::size_t i=0; i < origComm.indexSet().size(); i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
781 newComm.indexSet(), backwardscopyrowsize);
│ │ │ │ -
│ │ │ │ -
783 newComm.indexSet(),copyrowsize);
│ │ │ │ -
784 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy,
│ │ │ │ -
│ │ │ │ -
786 ri.getInterface().free();
│ │ │ │ -
787 RemoteIndices<IndexSet> *ris =
new RemoteIndices<IndexSet>(origComm.indexSet(),
│ │ │ │ -
│ │ │ │ -
789 origComm.communicator());
│ │ │ │ -
790 ris->template rebuild<true>();
│ │ │ │ -
791 ri.getInterface().build(*ris,ownerflags,ownerflags);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
795 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
│ │ │ │ -
│ │ │ │ -
797 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
│ │ │ │ -
798 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
819 template<
typename M,
typename C>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
837 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
845 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
Functionality for redistributing a parallel index set using graph partitioning.
│ │ │ │ -
Classes providing communication interfaces for overlapping Schwarz methods.
│ │ │ │ -
Col col
Definition matrixmatrix.hh:351
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
775 real_type norm_max = max(norm_dx, norm_dy);
│ │ │ │ +
776 real_type norm_min = min(norm_dx, norm_dy);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
779 cs = Simd::cond(norm_dy < eps,
│ │ │ │ +
│ │ │ │ +
781 Simd::cond(norm_dx < eps,
│ │ │ │ +
│ │ │ │ +
783 Simd::cond(norm_dy > norm_dx,
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
787 sn = Simd::cond(norm_dy < eps,
│ │ │ │ +
│ │ │ │ +
789 Simd::cond(norm_dx < eps,
│ │ │ │ +
│ │ │ │ +
791 Simd::cond(norm_dy > norm_dx,
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
802 using IterativeSolver<X,X>
::_op;
│ │ │ │ +
803 using IterativeSolver<X,X>
::_prec;
│ │ │ │ +
804 using IterativeSolver<X,X>
::_sp;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
825 template<
class X,
class Y=X,
class F = Y>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
926 const Simd::Scalar<real_type> EPSILON = 1e-80;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
930 std::vector<field_type,fAlloc> s(m+1), sn(m);
│ │ │ │ +
931 std::vector<real_type,rAlloc> cs(m);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
936 std::vector< std::vector<field_type,fAlloc> > H(m+1,s);
│ │ │ │ +
937 std::vector<F> v(m+1,b);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
945 _op->applyscaleadd(-1.0,x,b);
│ │ │ │ +
│ │ │ │ +
947 v[0] = 0.0;
_prec->apply(v[0],b);
│ │ │ │ +
948 norm =
_sp->norm(v[0]);
│ │ │ │ +
949 if(iteration.step(0, norm)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
969 _op->apply(v[i],v[i+1]);
│ │ │ │ +
970 _prec->apply(w,v[i+1]);
│ │ │ │ +
971 for(
int k=0; k<i+1; k++) {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
976 H[k][i] =
_sp->dot(v[k],w);
│ │ │ │ +
│ │ │ │ +
978 w.axpy(-H[k][i],v[k]);
│ │ │ │ +
│ │ │ │ +
980 H[i+1][i] =
_sp->norm(w);
│ │ │ │ +
981 if(Simd::allTrue(abs(H[i+1][i]) < EPSILON))
│ │ │ │ +
│ │ │ │ +
983 "breakdown in GMRes - |w| == 0.0 after " << j <<
" iterations");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
987 v[i+1] *= Simd::cond(norm==
real_type(0.),
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
992 for(
int k=0; k<i; k++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1004 iteration.step(j, norm);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1020 std::cout <<
"=== GMRes::restart" << std::endl;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1024 _op->applyscaleadd(-1.0,x,b);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1027 _prec->apply(v[0],b);
│ │ │ │ +
1028 norm =
_sp->norm(v[0]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1040 const std::vector<std::vector<field_type,fAlloc> >& H,
│ │ │ │ +
1041 const std::vector<field_type,fAlloc>& s,
│ │ │ │ +
1042 const std::vector<X>& v) {
│ │ │ │ +
│ │ │ │ +
1044 std::vector<field_type,fAlloc> y(s);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1047 for(
int a=i-1; a>=0; a--) {
│ │ │ │ +
│ │ │ │ +
1049 for(
int b=a+1; b<i; b++)
│ │ │ │ +
1050 rhs -= H[a][b]*y[b];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1061 template<
typename T>
│ │ │ │ +
│ │ │ │ +
1062 typename std::enable_if<std::is_same<field_type,real_type>::value,T>::type
conjugate(
const T& t) {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1066 template<
typename T>
│ │ │ │ +
│ │ │ │ +
1067 typename std::enable_if<!std::is_same<field_type,real_type>::value,T>::type
conjugate(
const T& t) {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1082 real_type norm_max = max(norm_dx, norm_dy);
│ │ │ │ +
1083 real_type norm_min = min(norm_dx, norm_dy);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1086 cs = Simd::cond(norm_dy < eps,
│ │ │ │ +
│ │ │ │ +
1088 Simd::cond(norm_dx < eps,
│ │ │ │ +
│ │ │ │ +
1090 Simd::cond(norm_dy > norm_dx,
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1094 sn = Simd::cond(norm_dy < eps,
│ │ │ │ +
│ │ │ │ +
1096 Simd::cond(norm_dx < eps,
│ │ │ │ +
│ │ │ │ +
1098 Simd::cond(norm_dy > norm_dx,
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1137 template<
class X,
class Y=X,
class F = Y>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1172 const Simd::Scalar<real_type> EPSILON = 1e-80;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1176 std::vector<field_type,fAlloc> s(m+1), sn(m);
│ │ │ │ +
1177 std::vector<real_type,rAlloc> cs(m);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1180 std::vector< std::vector<field_type,fAlloc> > H(m+1,s);
│ │ │ │ +
1181 std::vector<F> v(m+1,b);
│ │ │ │ +
1182 std::vector<X> w(m+1,b);
│ │ │ │ +
│ │ │ │ +
1184 Iteration iteration(*
this,res);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1190 _op->applyscaleadd(-1.0, x, v[0]);
│ │ │ │ +
│ │ │ │ +
1192 norm =
_sp->norm(v[0]);
│ │ │ │ +
1193 if(iteration.step(0, norm)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1202 v[0] *= (1.0 / norm);
│ │ │ │ +
│ │ │ │ +
1204 for(i=1; i<m+1; ++i)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1212 _prec->apply(w[i], v[i]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1215 _op->apply(w[i], v[i+1]);
│ │ │ │ +
│ │ │ │ +
1217 for(
int kk=0; kk<i+1; kk++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1223 H[kk][i] =
_sp->dot(v[kk],v[i+1]);
│ │ │ │ +
│ │ │ │ +
1225 v[i+1].axpy(-H[kk][i], v[kk]);
│ │ │ │ +
│ │ │ │ +
1227 H[i+1][i] =
_sp->norm(v[i+1]);
│ │ │ │ +
1228 if(Simd::allTrue(abs(H[i+1][i]) < EPSILON))
│ │ │ │ +
1229 DUNE_THROW(
SolverAbort,
"breakdown in fGMRes - |w| (-> "
│ │ │ │ +
1230 << w[i] <<
") == 0.0 after "
│ │ │ │ +
1231 << j <<
" iterations");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1250 iteration.step(j, norm);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1255 this->
update(tmp, i, H, s, w);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1265 std::cout <<
"=== fGMRes::restart" << std::endl;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1269 _op->applyscaleadd(-1.0, x,v[0]);
│ │ │ │ +
│ │ │ │ +
1271 norm =
_sp->norm(v[0]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1363 _restart(configuration.
get<int>(
"restart"))
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1368 _restart(configuration.
get<int>(
"restart"))
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1393 Iteration iteration(*
this, res);
│ │ │ │ +
│ │ │ │ +
1395 _op->applyscaleadd(-1,x,b);
│ │ │ │ +
│ │ │ │ +
1397 std::vector<std::shared_ptr<X> > p(_restart);
│ │ │ │ +
1398 std::vector<field_type,fAlloc> pp(_restart);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1402 p[0].reset(
new X(x));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1405 if(iteration.step(0, def)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1415 _prec->apply(*(p[0]),b);
│ │ │ │ +
1416 rho =
_sp->dot(*(p[0]),b);
│ │ │ │ +
1417 _op->apply(*(p[0]),q);
│ │ │ │ +
1418 pp[0] =
_sp->dot(*(p[0]),q);
│ │ │ │ +
│ │ │ │ +
1420 x.axpy(lambda,*(p[0]));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1426 if(iteration.step(i, def)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1433 int end=std::min(_restart,
_maxit-i+1);
│ │ │ │ +
1434 for (
int ii = 1; ii < end; ++ii)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1439 _prec->apply(prec_res,b);
│ │ │ │ +
│ │ │ │ +
1441 p[ii].reset(
new X(prec_res));
│ │ │ │ +
1442 _op->apply(prec_res, q);
│ │ │ │ +
│ │ │ │ +
1444 for(
int j=0; j<ii; ++j) {
│ │ │ │ +
1445 rho =
_sp->dot(q,*(p[j]))/pp[j];
│ │ │ │ +
1446 p[ii]->axpy(-rho, *(p[j]));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1450 _op->apply(*(p[ii]),q);
│ │ │ │ +
1451 pp[ii] =
_sp->dot(*(p[ii]),q);
│ │ │ │ +
1452 rho =
_sp->dot(*(p[ii]),b);
│ │ │ │ +
1453 lambda = rho/pp[ii];
│ │ │ │ +
1454 x.axpy(lambda,*(p[ii]));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1461 iteration.step(i, def);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1466 *(p[0])=*(p[_restart-1]);
│ │ │ │ +
1467 pp[0]=pp[_restart-1];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1529 scalar_real_type reduction,
int maxit,
int verbose,
int mmax = 10) :
IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
_mmax(mmax)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1560 const ParameterTree& config)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1567 const ParameterTree& config)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1589 _op->applyscaleadd(-1,x,b);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1592 std::vector<X> d(
_mmax+1, x);
│ │ │ │ +
1593 std::vector<X> Ad(
_mmax+1, x);
│ │ │ │ +
1594 std::vector<field_type,rAlloc> ddotAd(
_mmax+1,0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1598 if(iteration.step(0, def)){
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1610 for (; i_bounded <=
_mmax && i<=
_maxit; i_bounded++) {
│ │ │ │ +
│ │ │ │ +
1612 _prec->apply(d[i_bounded], b);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1615 orthogonalizations(i_bounded,Ad,w,ddotAd,d);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1618 _op->apply(d[i_bounded], Ad[i_bounded]);
│ │ │ │ +
1619 ddotAd[i_bounded]=
_sp->dot(d[i_bounded],Ad[i_bounded]);
│ │ │ │ +
1620 alpha =
_sp->dot(d[i_bounded], b)/ddotAd[i_bounded];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1623 x.axpy(alpha, d[i_bounded]);
│ │ │ │ +
1624 b.axpy(-alpha, Ad[i_bounded]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1629 iteration.step(i, def);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1633 cycle(Ad,d,ddotAd,i_bounded);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1646 for (
int k = 0; k < i_bounded; k++) {
│ │ │ │ +
1647 d[i_bounded].axpy(-
_sp->dot(Ad[k], w) / ddotAd[k], d[k]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1652 virtual void cycle(std::vector<X>& Ad,std::vector<X>& d,std::vector<
field_type,ReboundAllocatorType<X,field_type> >& ddotAd,
int& i_bounded) {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1655 std::swap(Ad[0], Ad[
_mmax]);
│ │ │ │ +
1656 std::swap(d[0], d[
_mmax]);
│ │ │ │ +
1657 std::swap(ddotAd[0], ddotAd[
_mmax]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1703 for (
int k = 0; k < _k_limit; k++) {
│ │ │ │ +
│ │ │ │ +
1705 d[i_bounded].axpy(-
_sp->dot(Ad[k], w) / ddotAd[k], d[k]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1708 if(_k_limit<=i_bounded)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1714 virtual void cycle(std::vector<X>& Ad, [[maybe_unused]] std::vector<X>& d, [[maybe_unused]] std::vector<
field_type,ReboundAllocatorType<X,field_type> >& ddotAd,
int& i_bounded)
override {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1718 _k_limit = Ad.size();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
#define DUNE_REGISTER_ITERATIVE_SOLVER(name,...)
Definition solverregistry.hh:19
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
Define general, extensible interface for operators. The available implementation wraps a matrix.
│ │ │ │ +
Define base class for scalar product and norm.
│ │ │ │ +
│ │ │ │ +
Define general, extensible interface for inverse operators.
│ │ │ │ +
Implementation of the BCRSMatrix class.
│ │ │ │ +
│ │ │ │
Definition allocator.hh:11
│ │ │ │ -
void redistributeMatrixEntries(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Definition matrixredistribute.hh:757
│ │ │ │ -
void redistributeSparsityPattern(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Definition matrixredistribute.hh:663
│ │ │ │ -
void redistributeMatrix(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Redistribute a matrix according to given domain decompositions.
Definition matrixredistribute.hh:820
│ │ │ │ -
derive error class from the base class in common
Definition istlexception.hh:19
│ │ │ │ -
Definition matrixredistribute.hh:22
│ │ │ │ -
void setNoBackwardsCopyRows(std::size_t size)
Definition matrixredistribute.hh:44
│ │ │ │ -
void redistribute(const D &from, D &to) const
Definition matrixredistribute.hh:28
│ │ │ │ -
void resetSetup()
Definition matrixredistribute.hh:35
│ │ │ │ -
void setNoCopyRows(std::size_t size)
Definition matrixredistribute.hh:41
│ │ │ │ -
bool isSetup() const
Definition matrixredistribute.hh:23
│ │ │ │ -
void setNoRows(std::size_t size)
Definition matrixredistribute.hh:38
│ │ │ │ -
void redistributeBackward(D &from, const D &to) const
Definition matrixredistribute.hh:32
│ │ │ │ -
std::size_t getBackwardsCopyRowSize(std::size_t index) const
Definition matrixredistribute.hh:57
│ │ │ │ -
std::size_t getRowSize(std::size_t index) const
Definition matrixredistribute.hh:47
│ │ │ │ -
std::size_t getCopyRowSize(std::size_t index) const
Definition matrixredistribute.hh:52
│ │ │ │ -
std::size_t getRowSize(std::size_t index) const
Definition matrixredistribute.hh:158
│ │ │ │ -
std::size_t & getBackwardsCopyRowSize(std::size_t index)
Definition matrixredistribute.hh:173
│ │ │ │ -
RedistributeInterface & getInterface()
Definition matrixredistribute.hh:75
│ │ │ │ -
void redistribute(const D &from, D &to) const
Definition matrixredistribute.hh:136
│ │ │ │ -
void setNoBackwardsCopyRows(std::size_t rows)
Definition matrixredistribute.hh:193
│ │ │ │ -
std::size_t & getCopyRowSize(std::size_t index)
Definition matrixredistribute.hh:163
│ │ │ │ -
RedistributeInformation()
Definition matrixredistribute.hh:71
│ │ │ │ -
std::size_t getCopyRowSize(std::size_t index) const
Definition matrixredistribute.hh:168
│ │ │ │ -
void setNoRows(std::size_t rows)
Definition matrixredistribute.hh:183
│ │ │ │ -
void reserve(std::size_t size)
Definition matrixredistribute.hh:150
│ │ │ │ -
OwnerOverlapCopyCommunication< T, T1 > Comm
Definition matrixredistribute.hh:69
│ │ │ │ -
void setNoCopyRows(std::size_t rows)
Definition matrixredistribute.hh:188
│ │ │ │ -
void redistributeBackward(D &from, const D &to) const
Definition matrixredistribute.hh:141
│ │ │ │ -
void setSetup()
Definition matrixredistribute.hh:106
│ │ │ │ -
std::size_t getBackwardsCopyRowSize(std::size_t index) const
Definition matrixredistribute.hh:178
│ │ │ │ -
void redistribute(const D &from, D &to) const
Definition matrixredistribute.hh:118
│ │ │ │ -
void resetSetup()
Definition matrixredistribute.hh:112
│ │ │ │ -
void redistributeBackward(D &from, const D &to) const
Definition matrixredistribute.hh:126
│ │ │ │ -
void checkInterface(const IS &source, const IS &target, MPI_Comm comm)
Definition matrixredistribute.hh:80
│ │ │ │ -
bool isSetup() const
Definition matrixredistribute.hh:145
│ │ │ │ -
std::size_t & getRowSize(std::size_t index)
Definition matrixredistribute.hh:153
│ │ │ │ -
Utility class to communicate and set the row sizes of a redistributed matrix.
Definition matrixredistribute.hh:216
│ │ │ │ -
M::size_type size_type
Definition matrixredistribute.hh:219
│ │ │ │ -
M::size_type value_type
Definition matrixredistribute.hh:218
│ │ │ │ -
RI & rowsize
Definition matrixredistribute.hh:230
│ │ │ │ -
const M & matrix
Definition matrixredistribute.hh:229
│ │ │ │ -
CommMatrixRowSize(const M &m_, RI &rowsize_)
Constructor.
Definition matrixredistribute.hh:226
│ │ │ │ -
Utility class to communicate and build the sparsity pattern of a redistributed matrix.
Definition matrixredistribute.hh:245
│ │ │ │ -
M::size_type size_type
Definition matrixredistribute.hh:246
│ │ │ │ -
const Dune::GlobalLookupIndexSet< I > & idxset
Definition matrixredistribute.hh:356
│ │ │ │ -
void storeSparsityPattern(M &m)
Creates and stores the sparsity pattern of the redistributed matrix.
Definition matrixredistribute.hh:276
│ │ │ │ -
const I & aggidxset
Definition matrixredistribute.hh:357
│ │ │ │ -
const std::vector< size_type > * rowsize
Definition matrixredistribute.hh:359
│ │ │ │ -
void completeSparsityPattern(std::vector< std::set< size_type > > add_sparsity)
Completes the sparsity pattern of the redistributed matrix with data from copy rows for the novlp cas...
Definition matrixredistribute.hh:340
│ │ │ │ -
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor for the original side.
Definition matrixredistribute.hh:254
│ │ │ │ -
const M & matrix
Definition matrixredistribute.hh:354
│ │ │ │ -
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, const std::vector< typename M::size_type > &rowsize_)
Constructor for the redistruted side.
Definition matrixredistribute.hh:265
│ │ │ │ -
std::vector< std::set< size_type > > sparsity
Definition matrixredistribute.hh:358
│ │ │ │ -
Dune::GlobalLookupIndexSet< I > LookupIndexSet
Definition matrixredistribute.hh:355
│ │ │ │ -
static M::size_type getSize(const Type &t, std::size_t i)
Definition matrixredistribute.hh:376
│ │ │ │ -
CommMatrixSparsityPattern< M, I > Type
Definition matrixredistribute.hh:365
│ │ │ │ -
I::GlobalIndex IndexedType
The indexed type we send. This is the global index indentitfying the column.
Definition matrixredistribute.hh:371
│ │ │ │ -
VariableSize IndexedTypeFlag
Each row varies in size.
Definition matrixredistribute.hh:374
│ │ │ │ -
Utility class for comunicating the matrix entries.
Definition matrixredistribute.hh:396
│ │ │ │ -
std::vector< size_t > * rowsize
row size information for the receiving side.
Definition matrixredistribute.hh:452
│ │ │ │ -
M & matrix
The matrix to communicate the values of.
Definition matrixredistribute.hh:446
│ │ │ │ -
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, std::vector< size_t > &rowsize_)
Constructor.
Definition matrixredistribute.hh:412
│ │ │ │ -
const Dune::GlobalLookupIndexSet< I > & idxset
Index set for the original matrix.
Definition matrixredistribute.hh:448
│ │ │ │ -
void setOverlapRowsToDirichlet()
Sets the non-owner rows correctly as Dirichlet boundaries.
Definition matrixredistribute.hh:421
│ │ │ │ -
const I & aggidxset
Index set for the redistributed matrix.
Definition matrixredistribute.hh:450
│ │ │ │ -
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor.
Definition matrixredistribute.hh:405
│ │ │ │ -
std::pair< typename I::GlobalIndex, typename M::block_type > IndexedType
The indexed type we send. This is the pair of global index indentitfying the column and the value its...
Definition matrixredistribute.hh:464
│ │ │ │ -
CommMatrixRow< M, I > Type
Definition matrixredistribute.hh:458
│ │ │ │ -
static std::size_t getSize(const Type &t, std::size_t i)
Definition matrixredistribute.hh:469
│ │ │ │ -
VariableSize IndexedTypeFlag
Each row varies in size.
Definition matrixredistribute.hh:467
│ │ │ │ -
Definition matrixredistribute.hh:483
│ │ │ │ -
static void scatter(Container &cont, const typename M::size_type &rowsize, std::size_t i)
Definition matrixredistribute.hh:490
│ │ │ │ -
static const M::size_type gather(const Container &cont, std::size_t i)
Definition matrixredistribute.hh:486
│ │ │ │ -
CommMatrixRowSize< M, RI > Container
Definition matrixredistribute.hh:484
│ │ │ │ -
Definition matrixredistribute.hh:500
│ │ │ │ -
static const M::size_type gather(const Container &cont, std::size_t i)
Definition matrixredistribute.hh:503
│ │ │ │ -
static void scatter(Container &cont, const typename M::size_type &rowsize, std::size_t i)
Definition matrixredistribute.hh:507
│ │ │ │ -
CommMatrixRowSize< M, RI > Container
Definition matrixredistribute.hh:501
│ │ │ │ -
Definition matrixredistribute.hh:518
│ │ │ │ -
M::ConstColIterator ColIter
Definition matrixredistribute.hh:521
│ │ │ │ -
static void scatter(Container &cont, const GlobalIndex &gi, std::size_t i, std::size_t j)
Definition matrixredistribute.hh:553
│ │ │ │ -
CommMatrixSparsityPattern< M, I > Container
Definition matrixredistribute.hh:520
│ │ │ │ -
static GlobalIndex numlimits
Definition matrixredistribute.hh:524
│ │ │ │ -
static ColIter col
Definition matrixredistribute.hh:523
│ │ │ │ -
I::GlobalIndex GlobalIndex
Definition matrixredistribute.hh:519
│ │ │ │ -
static const GlobalIndex & gather(const Container &cont, std::size_t i, std::size_t j)
Definition matrixredistribute.hh:526
│ │ │ │ -
Definition matrixredistribute.hh:599
│ │ │ │ -
I::GlobalIndex GlobalIndex
Definition matrixredistribute.hh:600
│ │ │ │ -
static Data datastore
Definition matrixredistribute.hh:605
│ │ │ │ -
static GlobalIndex numlimits
Definition matrixredistribute.hh:606
│ │ │ │ -
M::ConstColIterator ColIter
Definition matrixredistribute.hh:602
│ │ │ │ -
static const Data & gather(const Container &cont, std::size_t i, std::size_t j)
Definition matrixredistribute.hh:608
│ │ │ │ -
std::pair< GlobalIndex, typename M::block_type > Data
Definition matrixredistribute.hh:603
│ │ │ │ -
static void scatter(Container &cont, const Data &data, std::size_t i, std::size_t j)
Definition matrixredistribute.hh:638
│ │ │ │ -
static ColIter col
Definition matrixredistribute.hh:604
│ │ │ │ -
CommMatrixRow< M, I > Container
Definition matrixredistribute.hh:601
│ │ │ │ -
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition owneroverlapcopy.hh:174
│ │ │ │ -
EnumItem< AttributeSet, OwnerOverlapCopyAttributeSet::owner > OwnerSet
Definition owneroverlapcopy.hh:194
│ │ │ │ -
│ │ │ │ -
Definition repartition.hh:260
│ │ │ │ -
@ nonoverlapping
Category for non-overlapping solvers.
Definition solvercategory.hh:27
│ │ │ │ -
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition solvercategory.hh:34
│ │ │ │ +
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition dependency.hh:293
│ │ │ │ +
typename std::allocator_traits< typename AllocatorTraits< T >::type >::template rebind_alloc< X > ReboundAllocatorType
Definition allocator.hh:37
│ │ │ │ +
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:466
│ │ │ │ +
@ row_wise
Build in a row-wise manner.
Definition bcrsmatrix.hh:517
│ │ │ │ +
CreateIterator createend()
get create iterator pointing to one after the last block
Definition bcrsmatrix.hh:1100
│ │ │ │ +
CreateIterator createbegin()
get initial create iterator
Definition bcrsmatrix.hh:1094
│ │ │ │ +
size_type N() const
number of rows (counted in blocks)
Definition bcrsmatrix.hh:2001
│ │ │ │ +
A vector of blocks with memory management.
Definition bvector.hh:392
│ │ │ │ +
Wrapper to use a range of ARPACK++ eigenvalue solvers.
Definition arpackpp.hh:245
│ │ │ │ +
void computeSymMaxMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition arpackpp.hh:289
│ │ │ │ +
void computeSymMinMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition arpackpp.hh:391
│ │ │ │ +
Thrown when a solver aborts due to some problem.
Definition istlexception.hh:46
│ │ │ │ +
A linear operator.
Definition operators.hh:69
│ │ │ │ +
Base class for matrix free definition of preconditioners.
Definition preconditioner.hh:33
│ │ │ │ +
Base class for scalar product and norm computation.
Definition scalarproducts.hh:52
│ │ │ │ +
Statistics about the application of an inverse operator.
Definition solver.hh:50
│ │ │ │ +
double condition_estimate
Estimate of condition number.
Definition solver.hh:81
│ │ │ │ +
void clear()
Resets all data.
Definition solver.hh:58
│ │ │ │ +
bool converged
True if convergence criterion has been met.
Definition solver.hh:75
│ │ │ │ +
Simd::Scalar< real_type > scalar_real_type
scalar type underlying the field_type
Definition solver.hh:116
│ │ │ │ +
Y range_type
Type of the range of the operator to be inverted.
Definition solver.hh:107
│ │ │ │ +
X domain_type
Type of the domain of the operator to be inverted.
Definition solver.hh:104
│ │ │ │ +
X::field_type field_type
The field type of the operator.
Definition solver.hh:110
│ │ │ │ +
FieldTraits< field_type >::real_type real_type
The real type of the field type (is the same if using real numbers, but differs for std::complex)
Definition solver.hh:113
│ │ │ │ +
Base class for all implementations of iterative solvers.
Definition solver.hh:205
│ │ │ │ +
std::shared_ptr< const ScalarProduct< X > > _sp
Definition solver.hh:508
│ │ │ │ +
std::shared_ptr< const LinearOperator< X, X > > _op
Definition solver.hh:506
│ │ │ │ +
int _maxit
Definition solver.hh:510
│ │ │ │ +
int _verbose
Definition solver.hh:511
│ │ │ │ +
scalar_real_type _reduction
Definition solver.hh:509
│ │ │ │ +
std::shared_ptr< Preconditioner< X, X > > _prec
Definition solver.hh:507
│ │ │ │ +
Preconditioned loop solver.
Definition solvers.hh:59
│ │ │ │ +
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition solvers.hh:73
│ │ │ │ +
typename IterativeSolver< X, X >::template Iteration< unsigned int > Iteration
Definition solvers.hh:116
│ │ │ │ +
gradient method
Definition solvers.hh:124
│ │ │ │ +
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition solvers.hh:142
│ │ │ │ +
typename IterativeSolver< X, X >::template Iteration< unsigned int > Iteration
Definition solvers.hh:187
│ │ │ │ +
conjugate gradient method
Definition solvers.hh:193
│ │ │ │ +
CGSolver(std::shared_ptr< const LinearOperator< X, X > > op, std::shared_ptr< ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, X > > prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition solvers.hh:256
│ │ │ │ +
static constexpr bool enableConditionEstimate
Definition solvers.hh:208
│ │ │ │ +
CGSolver(const LinearOperator< X, X > &op, const ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition solvers.hh:239
│ │ │ │ +
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition solvers.hh:279
│ │ │ │ +
CGSolver(const LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition solvers.hh:222
│ │ │ │ +
typename IterativeSolver< X, X >::template Iteration< unsigned int > Iteration
Definition solvers.hh:412
│ │ │ │ +
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition solvers.hh:419
│ │ │ │ +
typename IterativeSolver< X, X >::template Iteration< CountType > Iteration
Definition solvers.hh:598
│ │ │ │ +
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition solvers.hh:439
│ │ │ │ +
Minimal Residual Method (MINRES)
Definition solvers.hh:609
│ │ │ │ +
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition solvers.hh:627
│ │ │ │ +
typename IterativeSolver< X, X >::template Iteration< unsigned int > Iteration
Definition solvers.hh:808
│ │ │ │ +
implements the Generalized Minimal Residual (GMRes) method
Definition solvers.hh:827
│ │ │ │ +
RestartedGMResSolver(std::shared_ptr< const LinearOperator< X, Y > > op, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Constructor.
Definition solvers.hh:878
│ │ │ │ +
std::enable_if<!std::is_same< field_type, real_type >::value, T >::type conjugate(const T &t)
Definition solvers.hh:1067
│ │ │ │ +
RestartedGMResSolver(std::shared_ptr< const LinearOperator< X, Y > > op, std::shared_ptr< const ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Definition solvers.hh:883
│ │ │ │ +
void update(X &w, int i, const std::vector< std::vector< field_type, fAlloc > > &H, const std::vector< field_type, fAlloc > &s, const std::vector< X > &v)
Definition solvers.hh:1039
│ │ │ │ +
std::enable_if< std::is_same< field_type, real_type >::value, T >::type conjugate(const T &t)
Definition solvers.hh:1062
│ │ │ │ +
ReboundAllocatorType< X, field_type > fAlloc
field_type Allocator retrieved from domain type
Definition solvers.hh:838
│ │ │ │ +
int _restart
Definition solvers.hh:1120
│ │ │ │ +
ReboundAllocatorType< X, real_type > rAlloc
real_type Allocator retrieved from domain type
Definition solvers.hh:840
│ │ │ │ +
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition solvers.hh:923
│ │ │ │ +
RestartedGMResSolver(const LinearOperator< X, Y > &op, Preconditioner< X, Y > &prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition solvers.hh:850
│ │ │ │ +
RestartedGMResSolver(const LinearOperator< X, Y > &op, const ScalarProduct< X > &sp, Preconditioner< X, Y > &prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition solvers.hh:861
│ │ │ │ +
void generatePlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
Definition solvers.hh:1073
│ │ │ │ +
RestartedGMResSolver(std::shared_ptr< const LinearOperator< X, Y > > op, std::shared_ptr< const ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, Y > > prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition solvers.hh:894
│ │ │ │ +
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator.
Definition solvers.hh:910
│ │ │ │ +
typename IterativeSolver< X, X >::template Iteration< unsigned int > Iteration
Definition solvers.hh:1119
│ │ │ │ +
void applyPlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
Definition solvers.hh:1106
│ │ │ │ +
implements the Flexible Generalized Minimal Residual (FGMRes) method (right preconditioned)
Definition solvers.hh:1139
│ │ │ │ +
void apply(X &x, Y &b, double reduction, InverseOperatorResult &res) override
Apply inverse operator.
Definition solvers.hh:1169
│ │ │ │ +
Generalized preconditioned conjugate gradient solver.
Definition solvers.hh:1307
│ │ │ │ +
GeneralizedPCGSolver(const LinearOperator< X, X > &op, const ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition solvers.hh:1343
│ │ │ │ +
GeneralizedPCGSolver(const LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition solvers.hh:1331
│ │ │ │ +
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition solvers.hh:1391
│ │ │ │ +
GeneralizedPCGSolver(std::shared_ptr< const LinearOperator< X, X > > op, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Constructor.
Definition solvers.hh:1361
│ │ │ │ +
GeneralizedPCGSolver(std::shared_ptr< const LinearOperator< X, X > > op, std::shared_ptr< const ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, X > > prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition solvers.hh:1377
│ │ │ │ +
GeneralizedPCGSolver(std::shared_ptr< const LinearOperator< X, X > > op, std::shared_ptr< const ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Definition solvers.hh:1366
│ │ │ │ +
Accelerated flexible conjugate gradient method.
Definition solvers.hh:1500
│ │ │ │ +
RestartedFCGSolver(const LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int mmax=10)
Constructor to initialize a RestartedFCG solver.
Definition solvers.hh:1518
│ │ │ │ +
RestartedFCGSolver(std::shared_ptr< const LinearOperator< X, X > > op, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &config)
Constructor.
Definition solvers.hh:1558
│ │ │ │ +
RestartedFCGSolver(std::shared_ptr< const LinearOperator< X, X > > op, std::shared_ptr< const ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, X > > prec, scalar_real_type reduction, int maxit, int verbose, int mmax=10)
Constructor to initialize a RestartedFCG solver.
Definition solvers.hh:1538
│ │ │ │ +
int _mmax
Definition solvers.hh:1661
│ │ │ │ +
RestartedFCGSolver(std::shared_ptr< const LinearOperator< X, X > > op, std::shared_ptr< const ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &config)
Definition solvers.hh:1564
│ │ │ │ +
typename IterativeSolver< X, X >::template Iteration< unsigned int > Iteration
Definition solvers.hh:1668
│ │ │ │ +
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition solvers.hh:1583
│ │ │ │ +
RestartedFCGSolver(const LinearOperator< X, X > &op, const ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int mmax=10)
Constructor to initialize a RestartedFCG solver.
Definition solvers.hh:1528
│ │ │ │ +
Complete flexible conjugate gradient method.
Definition solvers.hh:1679
│ │ │ │ +
virtual void apply(X &x, X &b, InverseOperatorResult &res) override
Apply inverse operator.
Definition solvers.hh:1693
│ │ │ │
