43 #ifndef IFPACK2_BLOCKRELAXATION_DEF_HPP 44 #define IFPACK2_BLOCKRELAXATION_DEF_HPP 47 #include "Ifpack2_LinearPartitioner.hpp" 48 #include "Ifpack2_LinePartitioner.hpp" 50 #include "Ifpack2_Details_UserPartitioner_def.hpp" 51 #include <Ifpack2_Parameters.hpp> 52 #include "Teuchos_TimeMonitor.hpp" 56 template<
class MatrixType,
class ContainerType>
58 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
60 if (A.getRawPtr () != A_.getRawPtr ()) {
61 IsInitialized_ =
false;
63 Partitioner_ = Teuchos::null;
67 IsParallel_ = (A->getRowMap ()->getComm ()->getSize () != 1);
68 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
69 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A);
70 hasBlockCrsMatrix_ = !A_bcrs.is_null();
72 if (! Container_.is_null ()) {
75 Container_->clearBlocks();
84 template<
class MatrixType,
class ContainerType>
90 PartitionerType_ (
"linear"),
93 containerType_ (
"TriDi"),
95 ZeroStartingSolution_ (true),
96 hasBlockCrsMatrix_ (false),
97 DoBackwardGS_ (false),
99 DampingFactor_ (STS::one ()),
100 IsInitialized_ (false),
105 InitializeTime_ (0.0),
109 NumGlobalNonzeros_ (0),
116 template<
class MatrixType,
class ContainerType>
121 template<
class MatrixType,
class ContainerType>
122 Teuchos::RCP<const Teuchos::ParameterList>
126 Teuchos::RCP<Teuchos::ParameterList> validParams = Teuchos::parameterList (
"Ifpack2::BlockRelaxation");
128 validParams->set(
"relaxation: container",
"TriDi");
129 validParams->set(
"relaxation: backward mode",
false);
130 validParams->set(
"relaxation: type",
"Jacobi");
131 validParams->set(
"relaxation: sweeps", 1);
132 validParams->set(
"relaxation: damping factor", STS::one());
133 validParams->set(
"relaxation: zero starting solution",
true);
134 validParams->set(
"block relaxation: decouple dofs",
false);
135 validParams->set(
"schwarz: compute condest",
false);
136 validParams->set(
"schwarz: combine mode",
"ZERO");
137 validParams->set(
"schwarz: use reordering",
true);
138 validParams->set(
"schwarz: filter singletons",
false);
139 validParams->set(
"schwarz: overlap level", 0);
140 validParams->set(
"partitioner: type",
"greedy");
141 validParams->set(
"partitioner: local parts", 1);
142 validParams->set(
"partitioner: overlap", 0);
143 Teuchos::Array<Teuchos::ArrayRCP<int>> tmp0;
144 validParams->set(
"partitioner: parts", tmp0);
145 validParams->set(
"partitioner: maintain sparsity",
false);
146 validParams->set(
"fact: ilut level-of-fill", 1.0);
147 validParams->set(
"fact: absolute threshold", 0.0);
148 validParams->set(
"fact: relative threshold", 1.0);
149 validParams->set(
"fact: relax value", 0.0);
151 Teuchos::ParameterList dummyList;
152 validParams->set(
"Amesos2",dummyList);
153 validParams->sublist(
"Amesos2").disableRecursiveValidation();
154 validParams->set(
"Amesos2 solver name",
"KLU2");
156 Teuchos::ArrayRCP<int> tmp;
157 validParams->set(
"partitioner: map", tmp);
159 validParams->set(
"partitioner: line detection threshold", 0.0);
160 validParams->set(
"partitioner: PDE equations", 1);
161 Teuchos::RCP<Tpetra::MultiVector<double,
162 typename MatrixType::local_ordinal_type,
163 typename MatrixType::global_ordinal_type,
164 typename MatrixType::node_type> > dummy;
165 validParams->set(
"partitioner: coordinates",dummy);
170 template<
class MatrixType,
class ContainerType>
176 Teuchos::RCP<const Teuchos::ParameterList> validparams;
178 List.validateParameters (*validparams);
180 if (List.isParameter (
"relaxation: container")) {
186 containerType_ = List.get<std::string> (
"relaxation: container");
188 if(containerType_ ==
"Tridiagonal") {
189 containerType_ =
"TriDi";
191 if(containerType_ ==
"Block Tridiagonal") {
192 containerType_ =
"BlockTriDi";
196 if (List.isParameter (
"relaxation: type")) {
197 const std::string relaxType = List.get<std::string> (
"relaxation: type");
199 if (relaxType ==
"Jacobi") {
200 PrecType_ = Ifpack2::Details::JACOBI;
202 else if (relaxType ==
"MT Split Jacobi") {
203 PrecType_ = Ifpack2::Details::MTSPLITJACOBI;
205 else if (relaxType ==
"Gauss-Seidel") {
206 PrecType_ = Ifpack2::Details::GS;
208 else if (relaxType ==
"Symmetric Gauss-Seidel") {
209 PrecType_ = Ifpack2::Details::SGS;
212 TEUCHOS_TEST_FOR_EXCEPTION
213 (
true, std::invalid_argument,
"Ifpack2::BlockRelaxation::" 214 "setParameters: Invalid parameter value \"" << relaxType
215 <<
"\" for parameter \"relaxation: type\".");
219 if (List.isParameter (
"relaxation: sweeps")) {
220 NumSweeps_ = List.get<
int> (
"relaxation: sweeps");
226 if (List.isParameter (
"relaxation: damping factor")) {
227 if (List.isType<
double> (
"relaxation: damping factor")) {
228 const double dampingFactor =
229 List.get<
double> (
"relaxation: damping factor");
230 DampingFactor_ =
static_cast<scalar_type> (dampingFactor);
232 else if (List.isType<
scalar_type> (
"relaxation: damping factor")) {
233 DampingFactor_ = List.get<
scalar_type> (
"relaxation: damping factor");
235 else if (List.isType<
magnitude_type> (
"relaxation: damping factor")) {
238 DampingFactor_ =
static_cast<scalar_type> (dampingFactor);
241 TEUCHOS_TEST_FOR_EXCEPTION
242 (
true, std::invalid_argument,
"Ifpack2::BlockRelaxation::" 243 "setParameters: Parameter \"relaxation: damping factor\" " 244 "has the wrong type.");
248 if (List.isParameter (
"relaxation: zero starting solution")) {
249 ZeroStartingSolution_ = List.get<
bool> (
"relaxation: zero starting solution");
252 if (List.isParameter (
"relaxation: backward mode")) {
253 DoBackwardGS_ = List.get<
bool> (
"relaxation: backward mode");
256 if (List.isParameter (
"partitioner: type")) {
257 PartitionerType_ = List.get<std::string> (
"partitioner: type");
262 if (List.isParameter (
"partitioner: local parts")) {
266 else if (! std::is_same<int, local_ordinal_type>::value &&
267 List.isType<
int> (
"partitioner: local parts")) {
268 NumLocalBlocks_ = List.get<
int> (
"partitioner: local parts");
271 TEUCHOS_TEST_FOR_EXCEPTION
272 (
true, std::invalid_argument,
"Ifpack2::BlockRelaxation::" 273 "setParameters: Parameter \"partitioner: local parts\" " 274 "has the wrong type.");
278 if (List.isParameter (
"partitioner: overlap level")) {
279 if (List.isType<
int> (
"partitioner: overlap level")) {
280 OverlapLevel_ = List.get<
int> (
"partitioner: overlap level");
282 else if (! std::is_same<int, local_ordinal_type>::value &&
287 TEUCHOS_TEST_FOR_EXCEPTION
288 (
true, std::invalid_argument,
"Ifpack2::BlockRelaxation::" 289 "setParameters: Parameter \"partitioner: overlap level\" " 290 "has the wrong type.");
294 std::string defaultContainer =
"TriDi";
301 if (PrecType_ != Ifpack2::Details::JACOBI) {
304 if (NumLocalBlocks_ < static_cast<local_ordinal_type> (0)) {
305 NumLocalBlocks_ = A_->getNodeNumRows() / (-NumLocalBlocks_);
309 if(List.isParameter(
"block relaxation: decouple dofs"))
310 decouple_ = List.get<
bool>(
"block relaxation: decouple dofs");
315 TEUCHOS_TEST_FOR_EXCEPTION(
316 DoBackwardGS_, std::runtime_error,
317 "Ifpack2::BlockRelaxation:setParameters: Setting the \"relaxation: " 318 "backward mode\" parameter to true is not yet supported.");
325 template<
class MatrixType,
class ContainerType>
326 Teuchos::RCP<const Teuchos::Comm<int> >
329 TEUCHOS_TEST_FOR_EXCEPTION
330 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::getComm: " 331 "The matrix is null. You must call setMatrix() with a nonnull matrix " 332 "before you may call this method.");
333 return A_->getComm ();
336 template<
class MatrixType,
class ContainerType>
337 Teuchos::RCP<
const Tpetra::RowMatrix<
typename MatrixType::scalar_type,
338 typename MatrixType::local_ordinal_type,
339 typename MatrixType::global_ordinal_type,
340 typename MatrixType::node_type> >
345 template<
class MatrixType,
class ContainerType>
346 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
347 typename MatrixType::global_ordinal_type,
348 typename MatrixType::node_type> >
352 TEUCHOS_TEST_FOR_EXCEPTION
353 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::" 354 "getDomainMap: The matrix is null. You must call setMatrix() with a " 355 "nonnull matrix before you may call this method.");
356 return A_->getDomainMap ();
359 template<
class MatrixType,
class ContainerType>
360 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
361 typename MatrixType::global_ordinal_type,
362 typename MatrixType::node_type> >
366 TEUCHOS_TEST_FOR_EXCEPTION
367 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::" 368 "getRangeMap: The matrix is null. You must call setMatrix() with a " 369 "nonnull matrix before you may call this method.");
370 return A_->getRangeMap ();
373 template<
class MatrixType,
class ContainerType>
380 template<
class MatrixType,
class ContainerType>
384 return NumInitialize_;
387 template<
class MatrixType,
class ContainerType>
395 template<
class MatrixType,
class ContainerType>
403 template<
class MatrixType,
class ContainerType>
408 return InitializeTime_;
411 template<
class MatrixType,
class ContainerType>
419 template<
class MatrixType,
class ContainerType>
428 template<
class MatrixType,
class ContainerType>
430 TEUCHOS_TEST_FOR_EXCEPTION(
431 A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::getNodeSmootherComplexity: " 432 "The input matrix A is null. Please call setMatrix() with a nonnull " 433 "input matrix, then call compute(), before calling this method.");
436 size_t block_nnz = 0;
438 block_nnz += Partitioner_->numRowsInPart(i) *Partitioner_->numRowsInPart(i);
440 return block_nnz + A_->getNodeNumEntries();
443 template<
class MatrixType,
class ContainerType>
446 apply (
const Tpetra::MultiVector<
typename MatrixType::scalar_type,
447 typename MatrixType::local_ordinal_type,
448 typename MatrixType::global_ordinal_type,
449 typename MatrixType::node_type>& X,
450 Tpetra::MultiVector<
typename MatrixType::scalar_type,
451 typename MatrixType::local_ordinal_type,
452 typename MatrixType::global_ordinal_type,
453 typename MatrixType::node_type>& Y,
454 Teuchos::ETransp mode,
458 TEUCHOS_TEST_FOR_EXCEPTION
459 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::apply: " 460 "The matrix is null. You must call setMatrix() with a nonnull matrix, " 461 "then call initialize() and compute() (in that order), before you may " 462 "call this method.");
463 TEUCHOS_TEST_FOR_EXCEPTION(
464 ! isComputed (), std::runtime_error,
"Ifpack2::BlockRelaxation::apply: " 465 "isComputed() must be true prior to calling apply.");
466 TEUCHOS_TEST_FOR_EXCEPTION(
467 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
468 "Ifpack2::BlockRelaxation::apply: X.getNumVectors() = " 469 << X.getNumVectors () <<
" != Y.getNumVectors() = " 470 << Y.getNumVectors () <<
".");
471 TEUCHOS_TEST_FOR_EXCEPTION(
472 mode != Teuchos::NO_TRANS, std::logic_error,
"Ifpack2::BlockRelaxation::" 473 "apply: This method currently only implements the case mode == Teuchos::" 474 "NO_TRANS. Transposed modes are not currently supported.");
475 TEUCHOS_TEST_FOR_EXCEPTION(
476 alpha != Teuchos::ScalarTraits<scalar_type>::one (), std::logic_error,
477 "Ifpack2::BlockRelaxation::apply: This method currently only implements " 478 "the case alpha == 1. You specified alpha = " << alpha <<
".");
479 TEUCHOS_TEST_FOR_EXCEPTION(
480 beta != Teuchos::ScalarTraits<scalar_type>::zero (), std::logic_error,
481 "Ifpack2::BlockRelaxation::apply: This method currently only implements " 482 "the case beta == 0. You specified beta = " << beta <<
".");
484 const std::string timerName (
"Ifpack2::BlockRelaxation::apply");
485 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
486 if (timer.is_null ()) {
487 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
490 double startTime = timer->wallTime();
493 Teuchos::TimeMonitor timeMon (*timer);
497 Teuchos::RCP<const MV> X_copy;
500 X_copy = rcp (
new MV (X, Teuchos::Copy));
502 X_copy = rcpFromRef (X);
506 if (ZeroStartingSolution_) {
507 Y.putScalar (STS::zero ());
512 case Ifpack2::Details::JACOBI:
513 ApplyInverseJacobi(*X_copy,Y);
515 case Ifpack2::Details::GS:
516 ApplyInverseGS(*X_copy,Y);
518 case Ifpack2::Details::SGS:
519 ApplyInverseSGS(*X_copy,Y);
521 case Ifpack2::Details::MTSPLITJACOBI:
523 Container_->applyInverseJacobi(*X_copy, Y, DampingFactor_, ZeroStartingSolution_, NumSweeps_);
526 TEUCHOS_TEST_FOR_EXCEPTION
527 (
true, std::logic_error,
"Ifpack2::BlockRelaxation::apply: Invalid " 528 "PrecType_ enum value " << PrecType_ <<
". Valid values are Ifpack2::" 529 "Details::JACOBI = " << Ifpack2::Details::JACOBI <<
", Ifpack2::Details" 530 "::GS = " << Ifpack2::Details::GS <<
", and Ifpack2::Details::SGS = " 531 << Ifpack2::Details::SGS <<
". Please report this bug to the Ifpack2 " 536 ApplyTime_ += (timer->wallTime() - startTime);
540 template<
class MatrixType,
class ContainerType>
543 applyMat (
const Tpetra::MultiVector<
typename MatrixType::scalar_type,
544 typename MatrixType::local_ordinal_type,
545 typename MatrixType::global_ordinal_type,
546 typename MatrixType::node_type>& X,
547 Tpetra::MultiVector<
typename MatrixType::scalar_type,
548 typename MatrixType::local_ordinal_type,
549 typename MatrixType::global_ordinal_type,
550 typename MatrixType::node_type>& Y,
551 Teuchos::ETransp mode)
const 553 A_->apply (X, Y, mode);
556 template<
class MatrixType,
class ContainerType>
562 typedef Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type>
565 TEUCHOS_TEST_FOR_EXCEPTION
566 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::initialize: " 567 "The matrix is null. You must call setMatrix() with a nonnull matrix " 568 "before you may call this method.");
570 Teuchos::RCP<Teuchos::Time> timer =
571 Teuchos::TimeMonitor::getNewCounter (
"Ifpack2::BlockRelaxation::initialize");
572 double startTime = timer->wallTime();
575 Teuchos::TimeMonitor timeMon (*timer);
576 IsInitialized_ =
false;
579 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
580 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
581 hasBlockCrsMatrix_ = !A_bcrs.is_null();
582 if (A_bcrs.is_null ()) {
583 hasBlockCrsMatrix_ =
false;
586 hasBlockCrsMatrix_ =
true;
589 NumLocalRows_ = A_->getNodeNumRows ();
590 NumGlobalRows_ = A_->getGlobalNumRows ();
591 NumGlobalNonzeros_ = A_->getGlobalNumEntries ();
596 Partitioner_ = Teuchos::null;
598 if (PartitionerType_ ==
"linear") {
601 }
else if (PartitionerType_ ==
"line") {
604 }
else if (PartitionerType_ ==
"user") {
610 TEUCHOS_TEST_FOR_EXCEPTION
611 (
true, std::logic_error,
"Ifpack2::BlockRelaxation::initialize: Unknown " 612 "partitioner type " << PartitionerType_ <<
". Valid values are " 613 "\"linear\", \"line\", and \"user\".");
617 Partitioner_->setParameters (List_);
618 Partitioner_->compute ();
621 NumLocalBlocks_ = Partitioner_->numLocalParts ();
626 if (A_->getComm()->getSize() != 1) {
634 TEUCHOS_TEST_FOR_EXCEPTION
635 (NumSweeps_ < 0, std::logic_error,
"Ifpack2::BlockRelaxation::initialize: " 636 "NumSweeps_ = " << NumSweeps_ <<
" < 0.");
639 ExtractSubmatricesStructure ();
643 if (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0) {
644 TEUCHOS_TEST_FOR_EXCEPTION
645 (hasBlockCrsMatrix_, std::runtime_error,
646 "Ifpack2::BlockRelaxation::initialize: " 647 "We do not support overlapped Jacobi yet for Tpetra::BlockCrsMatrix. Sorry!");
650 W_ = rcp (
new vector_type (A_->getRowMap ()));
651 W_->putScalar (STS::zero ());
652 Teuchos::ArrayRCP<scalar_type > w_ptr = W_->getDataNonConst(0);
655 for (
size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
657 w_ptr[LID] += STS::one();
660 W_->reciprocal (*W_);
664 InitializeTime_ += (timer->wallTime() - startTime);
666 IsInitialized_ =
true;
670 template<
class MatrixType,
class ContainerType>
677 TEUCHOS_TEST_FOR_EXCEPTION
678 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::compute: " 679 "The matrix is null. You must call setMatrix() with a nonnull matrix, " 680 "then call initialize(), before you may call this method.");
682 if (! isInitialized ()) {
686 Teuchos::RCP<Teuchos::Time> timer =
687 Teuchos::TimeMonitor::getNewCounter (
"Ifpack2::BlockRelaxation::compute");
689 double startTime = timer->wallTime();
692 Teuchos::TimeMonitor timeMon (*timer);
697 Container_->compute();
700 ComputeTime_ += (timer->wallTime() - startTime);
705 template<
class MatrixType,
class ContainerType>
710 TEUCHOS_TEST_FOR_EXCEPTION
711 (Partitioner_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::" 712 "ExtractSubmatricesStructure: Partitioner object is null.");
714 std::string containerType = ContainerType::getName ();
715 if (containerType ==
"Generic") {
719 containerType = containerType_;
723 bool pointIndexed = decouple_ && hasBlockCrsMatrix_;
724 Teuchos::Array<Teuchos::Array<local_ordinal_type> > blockRows;
728 auto A_bcrs = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(A_);
729 local_ordinal_type dofs = hasBlockCrsMatrix_ ?
730 A_bcrs->getBlockSize() : List_.get<
int>(
"partitioner: PDE equations");
731 blockRows.resize(NumLocalBlocks_ * dofs);
732 if(hasBlockCrsMatrix_)
734 for(local_ordinal_type i = 0; i < NumLocalBlocks_; i++)
736 size_t blockSize = Partitioner_->numRowsInPart(i);
739 for(local_ordinal_type j = 0; j < dofs; j++)
741 local_ordinal_type blockIndex = i * dofs + j;
742 blockRows[blockIndex].resize(blockSize);
743 for(
size_t k = 0; k < blockSize; k++)
747 blockRows[blockIndex][k] = (*Partitioner_)(i, k) * dofs + j;
756 for(local_ordinal_type i = 0; i < NumLocalBlocks_; i++)
759 TEUCHOS_TEST_FOR_EXCEPTION(Partitioner_->numRowsInPart(i) % dofs != 0, std::logic_error,
760 "Expected size of all blocks (partitions) to be divisible by MueLu dofs/node.");
761 size_t blockSize = Partitioner_->numRowsInPart(i) / dofs;
765 for(local_ordinal_type j = 0; j < dofs; j++)
767 local_ordinal_type blockIndex = i * dofs + j;
768 blockRows[blockIndex].resize(blockSize);
769 for(
size_t k = 0; k < blockSize; k++)
771 blockRows[blockIndex][k] = (*Partitioner_)(i, k * dofs + j);
780 blockRows.resize(NumLocalBlocks_);
781 for(local_ordinal_type i = 0; i < NumLocalBlocks_; ++i)
783 const size_t numRows = Partitioner_->numRowsInPart (i);
784 blockRows[i].resize(numRows);
786 for (
size_t j = 0; j < numRows; ++j)
788 blockRows[i][j] = (*Partitioner_) (i,j);
794 Container_->setParameters(List_);
795 Container_->initialize();
798 template<
class MatrixType,
class ContainerType>
800 BlockRelaxation<MatrixType,ContainerType>::
801 ApplyInverseJacobi (
const MV& X, MV& Y)
const 803 const size_t NumVectors = X.getNumVectors ();
805 MV AY (Y.getMap (), NumVectors);
808 int starting_iteration = 0;
809 if (OverlapLevel_ > 0)
812 auto WView = W_->getLocalViewHost (Tpetra::Access::ReadOnly);
813 if(ZeroStartingSolution_) {
814 auto XView = X.getLocalViewHost (Tpetra::Access::ReadOnly);
815 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
816 Container_->DoOverlappingJacobi(XView, YView, WView, DampingFactor_);
817 starting_iteration = 1;
819 const scalar_type ONE = STS::one();
820 for(
int j = starting_iteration; j < NumSweeps_; j++)
823 AY.update (ONE, X, -ONE);
825 auto AYView = AY.getLocalViewHost (Tpetra::Access::ReadOnly);
826 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
827 Container_->DoOverlappingJacobi (AYView, YView, WView, DampingFactor_);
834 if(ZeroStartingSolution_)
836 auto XView = X.getLocalViewHost (Tpetra::Access::ReadOnly);
837 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
838 Container_->DoJacobi (XView, YView, DampingFactor_);
839 starting_iteration = 1;
841 const scalar_type ONE = STS::one();
842 for(
int j = starting_iteration; j < NumSweeps_; j++)
845 AY.update (ONE, X, -ONE);
847 auto AYView = AY.getLocalViewHost (Tpetra::Access::ReadOnly);
848 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
849 Container_->DoJacobi (AYView, YView, DampingFactor_);
855 template<
class MatrixType,
class ContainerType>
857 BlockRelaxation<MatrixType,ContainerType>::
858 ApplyInverseGS (
const MV& X, MV& Y)
const 862 size_t numVecs = X.getNumVectors();
864 auto XView = X.getLocalViewHost (Tpetra::Access::ReadOnly);
865 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
868 bool deleteY2 =
false;
871 Y2 = ptr(
new MV(Importer_->getTargetMap(), numVecs));
878 for(
int j = 0; j < NumSweeps_; ++j)
881 Y2->doImport(Y, *Importer_, Tpetra::INSERT);
882 auto Y2View = Y2->getLocalViewHost (Tpetra::Access::ReadWrite);
883 Container_->DoGaussSeidel(XView, YView, Y2View, DampingFactor_);
888 auto Y2View = Y2->getLocalViewHost (Tpetra::Access::ReadWrite);
889 for(
int j = 0; j < NumSweeps_; ++j)
891 Container_->DoGaussSeidel(XView, YView, Y2View, DampingFactor_);
898 template<
class MatrixType,
class ContainerType>
900 BlockRelaxation<MatrixType,ContainerType>::
901 ApplyInverseSGS (
const MV& X, MV& Y)
const 906 auto XView = X.getLocalViewHost (Tpetra::Access::ReadOnly);
907 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
910 bool deleteY2 =
false;
913 Y2 = ptr(
new MV(Importer_->getTargetMap(), X.getNumVectors()));
920 for(
int j = 0; j < NumSweeps_; ++j)
923 Y2->doImport(Y, *Importer_, Tpetra::INSERT);
924 auto Y2View = Y2->getLocalViewHost (Tpetra::Access::ReadWrite);
925 Container_->DoSGS(XView, YView, Y2View, DampingFactor_);
930 auto Y2View = Y2->getLocalViewHost (Tpetra::Access::ReadWrite);
931 for(
int j = 0; j < NumSweeps_; ++j)
933 Container_->DoSGS(XView, YView, Y2View, DampingFactor_);
940 template<
class MatrixType,
class ContainerType>
941 void BlockRelaxation<MatrixType,ContainerType>::computeImporter ()
const 946 using Teuchos::ArrayView;
947 using Teuchos::Array;
949 using Teuchos::rcp_dynamic_cast;
952 if(hasBlockCrsMatrix_)
954 const RCP<const block_crs_matrix_type> bcrs =
955 rcp_dynamic_cast<
const block_crs_matrix_type>(A_);
956 int bs_ = bcrs->getBlockSize();
957 RCP<const map_type> oldDomainMap = A_->getDomainMap();
958 RCP<const map_type> oldColMap = A_->getColMap();
961 global_size_t numGlobalElements = oldColMap->getGlobalNumElements() * bs_;
962 global_ordinal_type indexBase = oldColMap->getIndexBase();
963 RCP<const Comm<int>> comm = oldColMap->getComm();
964 ArrayView<const global_ordinal_type> oldColElements = oldColMap->getNodeElementList();
965 Array<global_ordinal_type> newColElements(bs_ * oldColElements.size());
966 for(
int i = 0; i < oldColElements.size(); i++)
968 for(
int j = 0; j < bs_; j++)
969 newColElements[i * bs_ + j] = oldColElements[i] * bs_ + j;
971 RCP<map_type> colMap(
new map_type(numGlobalElements, newColElements, indexBase, comm));
973 Importer_ = rcp(
new import_type(oldDomainMap, colMap));
975 else if(!A_.is_null())
977 Importer_ = A_->getGraph()->getImporter();
978 if(Importer_.is_null())
979 Importer_ = rcp(
new import_type(A_->getDomainMap(), A_->getColMap()));
985 template<
class MatrixType,
class ContainerType>
990 std::ostringstream out;
995 out <<
"\"Ifpack2::BlockRelaxation\": {";
996 if (this->getObjectLabel () !=
"") {
997 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
1001 if (A_.is_null ()) {
1002 out <<
"Matrix: null, ";
1007 out <<
"Global matrix dimensions: [" 1008 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"], ";
1013 out <<
"\"relaxation: type\": ";
1014 if (PrecType_ == Ifpack2::Details::JACOBI) {
1015 out <<
"Block Jacobi";
1016 }
else if (PrecType_ == Ifpack2::Details::GS) {
1017 out <<
"Block Gauss-Seidel";
1018 }
else if (PrecType_ == Ifpack2::Details::SGS) {
1019 out <<
"Block Symmetric Gauss-Seidel";
1020 }
else if (PrecType_ == Ifpack2::Details::MTSPLITJACOBI) {
1021 out <<
"MT Split Jacobi";
1026 int approx_rows_per_part = A_->getNodeNumRows()/Partitioner_->numLocalParts();
1027 out <<
", blocksize: "<<approx_rows_per_part;
1029 out <<
", overlap: " << OverlapLevel_;
1031 out <<
", " <<
"sweeps: " << NumSweeps_ <<
", " 1032 <<
"damping factor: " << DampingFactor_ <<
", ";
1034 std::string containerType = ContainerType::getName();
1035 out <<
"block container: " << ((containerType ==
"Generic") ? containerType_ : containerType);
1036 if(List_.isParameter(
"partitioner: PDE equations"))
1037 out <<
", dofs/node: "<<List_.get<
int>(
"partitioner: PDE equations");
1044 template<
class MatrixType,
class ContainerType>
1048 const Teuchos::EVerbosityLevel verbLevel)
const 1052 using Teuchos::TypeNameTraits;
1053 using Teuchos::VERB_DEFAULT;
1054 using Teuchos::VERB_NONE;
1055 using Teuchos::VERB_LOW;
1056 using Teuchos::VERB_MEDIUM;
1057 using Teuchos::VERB_HIGH;
1058 using Teuchos::VERB_EXTREME;
1060 Teuchos::EVerbosityLevel vl = verbLevel;
1061 if (vl == VERB_DEFAULT) vl = VERB_LOW;
1062 const int myImageID = A_->getComm()->getRank();
1065 Teuchos::OSTab tab (out);
1072 if (vl != VERB_NONE && myImageID == 0) {
1073 out <<
"Ifpack2::BlockRelaxation<" 1074 << TypeNameTraits<MatrixType>::name () <<
", " 1075 << TypeNameTraits<ContainerType>::name () <<
" >:";
1076 Teuchos::OSTab tab1 (out);
1078 if (this->getObjectLabel () !=
"") {
1079 out <<
"label: \"" << this->getObjectLabel () <<
"\"" << endl;
1081 out <<
"initialized: " << (isInitialized () ?
"true" :
"false") << endl
1082 <<
"computed: " << (isComputed () ?
"true" :
"false") << endl;
1084 std::string precType;
1085 if (PrecType_ == Ifpack2::Details::JACOBI) {
1086 precType =
"Block Jacobi";
1087 }
else if (PrecType_ == Ifpack2::Details::GS) {
1088 precType =
"Block Gauss-Seidel";
1089 }
else if (PrecType_ == Ifpack2::Details::SGS) {
1090 precType =
"Block Symmetric Gauss-Seidel";
1092 out <<
"type: " << precType << endl
1093 <<
"global number of rows: " << A_->getGlobalNumRows () << endl
1094 <<
"global number of columns: " << A_->getGlobalNumCols () << endl
1095 <<
"number of sweeps: " << NumSweeps_ << endl
1096 <<
"damping factor: " << DampingFactor_ << endl
1097 <<
"backwards mode: " 1098 << ((PrecType_ == Ifpack2::Details::GS && DoBackwardGS_) ?
"true" :
"false")
1100 <<
"zero starting solution: " 1101 << (ZeroStartingSolution_ ?
"true" :
"false") << endl;
1119 #ifdef HAVE_IFPACK2_EXPLICIT_INSTANTIATION 1121 #define IFPACK2_BLOCKRELAXATION_INSTANT(S,LO,GO,N) \ 1123 class Ifpack2::BlockRelaxation< \ 1124 Tpetra::RowMatrix<S, LO, GO, N>, \ 1125 Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N> > >; 1127 #endif // HAVE_IFPACK2_EXPLICIT_INSTANTIATION 1129 #endif // IFPACK2_BLOCKRELAXATION_DEF_HPP Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the input matrix is distributed.
Definition: Ifpack2_BlockRelaxation_def.hpp:327
Ifpack2::BlockRelaxation class declaration.
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:422
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_BlockRelaxation_def.hpp:429
double getComputeTime() const
Returns the time spent in compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:414
int getNumInitialize() const
Returns the number of calls to initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:383
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Definition: Ifpack2_BlockRelaxation_def.hpp:1047
double getInitializeTime() const
Returns the time spent in initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:406
Teuchos::RCP< const row_matrix_type > getMatrix() const
The input matrix of this preconditioner's constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:341
Teuchos::RCP< const map_type > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:350
void getParameter(const Teuchos::ParameterList ¶ms, const std::string &name, T &value)
Set a value from a ParameterList if a parameter with the specified name exists.
Definition: Ifpack2_Parameters.hpp:59
void apply(const MV &X, MV &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Applies the preconditioner to X, returns the result in Y.
Definition: Ifpack2_BlockRelaxation_def.hpp:446
std::string description() const
A one-line description of this object.
Definition: Ifpack2_BlockRelaxation_def.hpp:988
Teuchos::RCP< const map_type > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:364
void setParameters(const Teuchos::ParameterList ¶ms)
Sets all the parameters for the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:173
void compute()
compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation ...
Definition: Ifpack2_BlockRelaxation_def.hpp:673
Block relaxation preconditioners (or smoothers) for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse ma...
Definition: Ifpack2_BlockRelaxation_decl.hpp:82
Ifpack2 implementation details.
Declaration of a user-defined partitioner in which the user can define a partition of the graph...
void applyMat(const MV &X, MV &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS) const
Applies the matrix to a Tpetra::MultiVector.
Definition: Ifpack2_BlockRelaxation_def.hpp:543
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_BlockRelaxation_decl.hpp:108
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_BlockRelaxation_def.hpp:58
virtual ~BlockRelaxation()
Destructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:118
Partition in which the user can define a nonoverlapping partition of the graph in any way they choose...
Definition: Ifpack2_Details_UserPartitioner_decl.hpp:69
void initialize()
Initialize.
Definition: Ifpack2_BlockRelaxation_def.hpp:559
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:100
int getNumCompute() const
Returns the number of calls to compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:390
Definition: Ifpack2_Container_decl.hpp:576
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:398
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:97
Interface for creating and solving a set of local linear problems.
Definition: Ifpack2_Container_decl.hpp:112
static Teuchos::RCP< BaseContainer > build(std::string containerType, const Teuchos::RCP< const MatrixType > &A, const Teuchos::Array< Teuchos::Array< local_ordinal_type >> &partitions, const Teuchos::RCP< const import_type > importer, bool pointIndexed)
Build a specialization of Ifpack2::Container given a key that has been registered.
Definition: Ifpack2_ContainerFactory_def.hpp:89
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
void getValidParameters(Teuchos::ParameterList ¶ms)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:51
Ifpack2::LinePartitioner: A class to define partitions into a set of lines.
Definition: Ifpack2_LinePartitioner_decl.hpp:77
A class to define linear partitions.
Definition: Ifpack2_LinearPartitioner_decl.hpp:60
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a list of all the parameters that this class accepts.
Definition: Ifpack2_BlockRelaxation_def.hpp:124
BlockRelaxation(const Teuchos::RCP< const row_matrix_type > &Matrix)
Constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:86