43 #ifndef IFPACK2_BLOCKRELAXATION_DEF_HPP
44 #define IFPACK2_BLOCKRELAXATION_DEF_HPP
47 #include "Ifpack2_LinearPartitioner.hpp"
48 #include "Ifpack2_LinePartitioner.hpp"
49 #include "Ifpack2_Zoltan2Partitioner_decl.hpp"
50 #include "Ifpack2_Zoltan2Partitioner_def.hpp"
52 #include "Ifpack2_Details_UserPartitioner_def.hpp"
53 #include "Ifpack2_LocalFilter.hpp"
54 #include "Ifpack2_Parameters.hpp"
56 #include "Tpetra_BlockCrsMatrix_Helpers_decl.hpp"
60 template<
class MatrixType,
class ContainerType>
65 IsInitialized_ =
false;
67 Partitioner_ = Teuchos::null;
71 IsParallel_ = (A->getRowMap ()->getComm ()->getSize () != 1);
73 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A);
74 hasBlockCrsMatrix_ = !A_bcrs.
is_null();
76 if (! Container_.is_null ()) {
79 Container_->clearBlocks();
88 template<
class MatrixType,
class ContainerType>
92 Container_ (Teuchos::null),
93 Partitioner_ (Teuchos::null),
94 PartitionerType_ (
"linear"),
97 containerType_ (
"TriDi"),
98 PrecType_ (Ifpack2::Details::JACOBI),
99 ZeroStartingSolution_ (true),
100 hasBlockCrsMatrix_ (false),
101 DoBackwardGS_ (false),
104 schwarzCombineMode_(
"ZERO"),
105 DampingFactor_ (
STS::one ()),
106 IsInitialized_ (false),
110 TimerForApply_(true),
112 InitializeTime_ (0.0),
116 NumGlobalNonzeros_ (0),
118 Importer_ (Teuchos::null)
123 template<
class MatrixType,
class ContainerType>
128 template<
class MatrixType,
class ContainerType>
135 validParams->
set(
"relaxation: container",
"TriDi");
136 validParams->
set(
"relaxation: backward mode",
false);
137 validParams->
set(
"relaxation: type",
"Jacobi");
138 validParams->
set(
"relaxation: sweeps", 1);
139 validParams->
set(
"relaxation: damping factor", STS::one());
140 validParams->
set(
"relaxation: zero starting solution",
true);
141 validParams->
set(
"block relaxation: decouple dofs",
false);
142 validParams->
set(
"schwarz: compute condest",
false);
143 validParams->
set(
"schwarz: combine mode",
"ZERO");
144 validParams->
set(
"schwarz: use reordering",
true);
145 validParams->
set(
"schwarz: filter singletons",
false);
146 validParams->
set(
"schwarz: overlap level", 0);
147 validParams->
set(
"partitioner: type",
"greedy");
148 validParams->
set(
"zoltan2: algorithm",
"phg");
149 validParams->
set(
"partitioner: local parts", 1);
150 validParams->
set(
"partitioner: overlap", 0);
151 validParams->
set(
"partitioner: combine mode",
"ZERO");
153 validParams->
set(
"partitioner: parts", tmp0);
155 validParams->
set(
"partitioner: global ID parts", tmp1);
156 validParams->
set(
"partitioner: nonsymmetric overlap combine",
false);
157 validParams->
set(
"partitioner: maintain sparsity",
false);
158 validParams->
set(
"fact: ilut level-of-fill", 1.0);
159 validParams->
set(
"fact: absolute threshold", 0.0);
160 validParams->
set(
"fact: relative threshold", 1.0);
161 validParams->
set(
"fact: relax value", 0.0);
164 validParams->
set(
"Amesos2",dummyList);
166 validParams->
set(
"Amesos2 solver name",
"KLU2");
169 validParams->
set(
"partitioner: map", tmp);
171 validParams->
set(
"partitioner: line detection threshold", 0.0);
172 validParams->
set(
"partitioner: PDE equations", 1);
174 typename MatrixType::local_ordinal_type,
175 typename MatrixType::global_ordinal_type,
176 typename MatrixType::node_type> > dummy;
177 validParams->
set(
"partitioner: coordinates",dummy);
178 validParams->
set(
"timer for apply",
true);
179 validParams->
set(
"partitioner: subparts per part", 1);
180 validParams->
set(
"partitioner: block size", -1);
181 validParams->
set(
"partitioner: print level",
false);
186 template<
class MatrixType,
class ContainerType>
194 this->setParametersImpl(const_cast<Teuchos::ParameterList&>(pl));
197 template<
class MatrixType,
class ContainerType>
202 if (List.
isType<
double>(
"relaxation: damping factor")) {
205 scalar_type df = List.
get<
double>(
"relaxation: damping factor");
206 List.
remove(
"relaxation: damping factor");
207 List.
set(
"relaxation: damping factor",df);
221 containerType_ = List.
get<std::string> (
"relaxation: container");
223 if(containerType_ ==
"Tridiagonal") {
224 containerType_ =
"TriDi";
226 if(containerType_ ==
"Block Tridiagonal") {
227 containerType_ =
"BlockTriDi";
232 const std::string relaxType = List.
get<std::string> (
"relaxation: type");
234 if (relaxType ==
"Jacobi") {
235 PrecType_ = Ifpack2::Details::JACOBI;
237 else if (relaxType ==
"MT Split Jacobi") {
238 PrecType_ = Ifpack2::Details::MTSPLITJACOBI;
240 else if (relaxType ==
"Gauss-Seidel") {
241 PrecType_ = Ifpack2::Details::GS;
243 else if (relaxType ==
"Symmetric Gauss-Seidel") {
244 PrecType_ = Ifpack2::Details::SGS;
248 (
true, std::invalid_argument,
"Ifpack2::BlockRelaxation::"
249 "setParameters: Invalid parameter value \"" << relaxType
250 <<
"\" for parameter \"relaxation: type\".");
255 NumSweeps_ = List.
get<
int> (
"relaxation: sweeps");
261 if (List.
isParameter (
"relaxation: damping factor")) {
262 if (List.
isType<
double> (
"relaxation: damping factor")) {
263 const double dampingFactor =
264 List.
get<
double> (
"relaxation: damping factor");
265 DampingFactor_ =
static_cast<scalar_type
> (dampingFactor);
267 else if (List.
isType<scalar_type> (
"relaxation: damping factor")) {
268 DampingFactor_ = List.
get<scalar_type> (
"relaxation: damping factor");
270 else if (List.
isType<magnitude_type> (
"relaxation: damping factor")) {
271 const magnitude_type dampingFactor =
272 List.
get<magnitude_type> (
"relaxation: damping factor");
273 DampingFactor_ =
static_cast<scalar_type
> (dampingFactor);
277 (
true, std::invalid_argument,
"Ifpack2::BlockRelaxation::"
278 "setParameters: Parameter \"relaxation: damping factor\" "
279 "has the wrong type.");
283 if (List.
isParameter (
"relaxation: zero starting solution")) {
284 ZeroStartingSolution_ = List.
get<
bool> (
"relaxation: zero starting solution");
287 if (List.
isParameter (
"relaxation: backward mode")) {
288 DoBackwardGS_ = List.
get<
bool> (
"relaxation: backward mode");
292 PartitionerType_ = List.
get<std::string> (
"partitioner: type");
297 if (List.
isParameter (
"partitioner: local parts")) {
298 if (List.
isType<local_ordinal_type> (
"partitioner: local parts")) {
299 NumLocalBlocks_ = List.
get<local_ordinal_type> (
"partitioner: local parts");
301 else if (! std::is_same<int, local_ordinal_type>::value &&
302 List.
isType<
int> (
"partitioner: local parts")) {
303 NumLocalBlocks_ = List.
get<
int> (
"partitioner: local parts");
307 (
true, std::invalid_argument,
"Ifpack2::BlockRelaxation::"
308 "setParameters: Parameter \"partitioner: local parts\" "
309 "has the wrong type.");
313 if (List.
isParameter (
"partitioner: overlap level")) {
314 if (List.
isType<
int> (
"partitioner: overlap level")) {
315 OverlapLevel_ = List.
get<
int> (
"partitioner: overlap level");
317 else if (! std::is_same<int, local_ordinal_type>::value &&
318 List.
isType<local_ordinal_type> (
"partitioner: overlap level")) {
319 OverlapLevel_ = List.
get<local_ordinal_type> (
"partitioner: overlap level");
323 (
true, std::invalid_argument,
"Ifpack2::BlockRelaxation::"
324 "setParameters: Parameter \"partitioner: overlap level\" "
325 "has the wrong type.");
330 if ( ( List.
isParameter(
"partitioner: global ID parts")) && (OverlapLevel_ < 1)) OverlapLevel_ = 1;
332 if (List.
isParameter (
"partitioner: nonsymmetric overlap combine"))
333 nonsymCombine_ = List.
get<
bool> (
"partitioner: nonsymmetric overlap combine");
335 if (List.
isParameter (
"partitioner: combine mode"))
336 schwarzCombineMode_ = List.
get<std::string> (
"partitioner: combine mode");
338 std::string defaultContainer =
"TriDi";
339 if(std::is_same<ContainerType, Container<MatrixType> >::value)
345 if (PrecType_ != Ifpack2::Details::JACOBI) {
348 if (NumLocalBlocks_ < static_cast<local_ordinal_type> (0)) {
349 NumLocalBlocks_ = A_->getLocalNumRows() / (-NumLocalBlocks_);
353 if(List.
isParameter(
"block relaxation: decouple dofs"))
354 decouple_ = List.
get<
bool>(
"block relaxation: decouple dofs");
360 DoBackwardGS_, std::runtime_error,
361 "Ifpack2::BlockRelaxation:setParameters: Setting the \"relaxation: "
362 "backward mode\" parameter to true is not yet supported.");
365 TimerForApply_ = List.
get<
bool>(
"timer for apply");
372 template<
class MatrixType,
class ContainerType>
377 (A_.
is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::getComm: "
378 "The matrix is null. You must call setMatrix() with a nonnull matrix "
379 "before you may call this method.");
380 return A_->getComm ();
383 template<
class MatrixType,
class ContainerType>
384 Teuchos::RCP<
const Tpetra::RowMatrix<
typename MatrixType::scalar_type,
385 typename MatrixType::local_ordinal_type,
386 typename MatrixType::global_ordinal_type,
387 typename MatrixType::node_type> >
392 template<
class MatrixType,
class ContainerType>
393 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
394 typename MatrixType::global_ordinal_type,
395 typename MatrixType::node_type> >
400 (A_.
is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::"
401 "getDomainMap: The matrix is null. You must call setMatrix() with a "
402 "nonnull matrix before you may call this method.");
403 return A_->getDomainMap ();
406 template<
class MatrixType,
class ContainerType>
407 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
408 typename MatrixType::global_ordinal_type,
409 typename MatrixType::node_type> >
414 (A_.
is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::"
415 "getRangeMap: The matrix is null. You must call setMatrix() with a "
416 "nonnull matrix before you may call this method.");
417 return A_->getRangeMap ();
420 template<
class MatrixType,
class ContainerType>
427 template<
class MatrixType,
class ContainerType>
431 return NumInitialize_;
434 template<
class MatrixType,
class ContainerType>
442 template<
class MatrixType,
class ContainerType>
450 template<
class MatrixType,
class ContainerType>
455 return InitializeTime_;
458 template<
class MatrixType,
class ContainerType>
466 template<
class MatrixType,
class ContainerType>
475 template<
class MatrixType,
class ContainerType>
478 A_.
is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::getNodeSmootherComplexity: "
479 "The input matrix A is null. Please call setMatrix() with a nonnull "
480 "input matrix, then call compute(), before calling this method.");
483 size_t block_nnz = 0;
485 block_nnz += Partitioner_->numRowsInPart(i) *Partitioner_->numRowsInPart(i);
487 return block_nnz + A_->getLocalNumEntries();
490 template<
class MatrixType,
class ContainerType>
493 apply (
const Tpetra::MultiVector<
typename MatrixType::scalar_type,
494 typename MatrixType::local_ordinal_type,
495 typename MatrixType::global_ordinal_type,
496 typename MatrixType::node_type>& X,
497 Tpetra::MultiVector<
typename MatrixType::scalar_type,
498 typename MatrixType::local_ordinal_type,
499 typename MatrixType::global_ordinal_type,
500 typename MatrixType::node_type>& Y,
506 (A_.
is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::apply: "
507 "The matrix is null. You must call setMatrix() with a nonnull matrix, "
508 "then call initialize() and compute() (in that order), before you may "
509 "call this method.");
511 ! isComputed (), std::runtime_error,
"Ifpack2::BlockRelaxation::apply: "
512 "isComputed() must be true prior to calling apply.");
513 TEUCHOS_TEST_FOR_EXCEPTION(
514 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
515 "Ifpack2::BlockRelaxation::apply: X.getNumVectors() = "
516 << X.getNumVectors () <<
" != Y.getNumVectors() = "
517 << Y.getNumVectors () <<
".");
518 TEUCHOS_TEST_FOR_EXCEPTION(
520 "apply: This method currently only implements the case mode == Teuchos::"
521 "NO_TRANS. Transposed modes are not currently supported.");
522 TEUCHOS_TEST_FOR_EXCEPTION(
524 "Ifpack2::BlockRelaxation::apply: This method currently only implements "
525 "the case alpha == 1. You specified alpha = " << alpha <<
".");
526 TEUCHOS_TEST_FOR_EXCEPTION(
528 "Ifpack2::BlockRelaxation::apply: This method currently only implements "
529 "the case beta == 0. You specified beta = " << beta <<
".");
531 const std::string timerName (
"Ifpack2::BlockRelaxation::apply");
533 if (TimerForApply_) {
555 X_copy = rcpFromRef (X);
559 if (ZeroStartingSolution_) {
560 Y.putScalar (STS::zero ());
565 case Ifpack2::Details::JACOBI:
566 ApplyInverseJacobi(*X_copy,Y);
568 case Ifpack2::Details::GS:
569 ApplyInverseGS(*X_copy,Y);
571 case Ifpack2::Details::SGS:
572 ApplyInverseSGS(*X_copy,Y);
574 case Ifpack2::Details::MTSPLITJACOBI:
576 Container_->applyInverseJacobi(*X_copy, Y, DampingFactor_, ZeroStartingSolution_, NumSweeps_);
579 TEUCHOS_TEST_FOR_EXCEPTION
580 (
true, std::logic_error,
"Ifpack2::BlockRelaxation::apply: Invalid "
581 "PrecType_ enum value " << PrecType_ <<
". Valid values are Ifpack2::"
582 "Details::JACOBI = " << Ifpack2::Details::JACOBI <<
", Ifpack2::Details"
583 "::GS = " << Ifpack2::Details::GS <<
", and Ifpack2::Details::SGS = "
584 << Ifpack2::Details::SGS <<
". Please report this bug to the Ifpack2 "
589 ApplyTime_ += (time.
wallTime() - startTime);
593 template<
class MatrixType,
class ContainerType>
596 applyMat (
const Tpetra::MultiVector<
typename MatrixType::scalar_type,
597 typename MatrixType::local_ordinal_type,
598 typename MatrixType::global_ordinal_type,
599 typename MatrixType::node_type>& X,
600 Tpetra::MultiVector<
typename MatrixType::scalar_type,
601 typename MatrixType::local_ordinal_type,
602 typename MatrixType::global_ordinal_type,
603 typename MatrixType::node_type>& Y,
606 A_->apply (X, Y, mode);
609 template<
class MatrixType,
class ContainerType>
615 typedef Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type>
619 (A_.
is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::initialize: "
620 "The matrix is null. You must call setMatrix() with a nonnull matrix "
621 "before you may call this method.");
625 double startTime = timer->wallTime();
629 IsInitialized_ =
false;
633 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
634 hasBlockCrsMatrix_ = !A_bcrs.
is_null();
636 if(!hasBlockCrsMatrix_ && List_.isParameter(
"relaxation: container") && List_.get<std::string>(
"relaxation: container") ==
"BlockTriDi" ) {
638 int block_size = List_.get<
int>(
"partitioner: block size");
640 (block_size == -1,
"A pointwise matrix and block_size = -1 were given as inputs.");
641 A_bcrs = Tpetra::convertToBlockCrsMatrix(*Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_), block_size);
643 hasBlockCrsMatrix_ =
true;
644 Kokkos::DefaultExecutionSpace().fence();
647 NumLocalRows_ = A_->getLocalNumRows ();
648 NumGlobalRows_ = A_->getGlobalNumRows ();
649 NumGlobalNonzeros_ = A_->getGlobalNumEntries ();
654 Partitioner_ = Teuchos::null;
656 if (PartitionerType_ ==
"linear") {
659 }
else if (PartitionerType_ ==
"line") {
662 }
else if (PartitionerType_ ==
"user") {
665 }
else if (PartitionerType_ ==
"zoltan2") {
666 #if defined(HAVE_IFPACK2_ZOLTAN2)
667 if (A_->getGraph ()->getComm ()->getSize () == 1) {
670 rcp (
new Ifpack2::Zoltan2Partitioner<row_graph_type> (A_->getGraph ()) );
675 rcp (
new Ifpack2::Zoltan2Partitioner<row_graph_type> (A_local->getGraph ()) );
679 (
true, std::logic_error,
"Ifpack2::BlockRelaxation::initialize: Zoltan2 not enabled.");
685 (
true, std::logic_error,
"Ifpack2::BlockRelaxation::initialize: Unknown "
686 "partitioner type " << PartitionerType_ <<
". Valid values are "
687 "\"linear\", \"line\", and \"user\".");
691 Partitioner_->setParameters (List_);
692 Partitioner_->compute ();
695 NumLocalBlocks_ = Partitioner_->numLocalParts ();
700 if (A_->getComm()->getSize() != 1) {
709 (NumSweeps_ < 0, std::logic_error,
"Ifpack2::BlockRelaxation::initialize: "
710 "NumSweeps_ = " << NumSweeps_ <<
" < 0.");
713 ExtractSubmatricesStructure ();
717 if (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0) {
719 (hasBlockCrsMatrix_, std::runtime_error,
720 "Ifpack2::BlockRelaxation::initialize: "
721 "We do not support overlapped Jacobi yet for Tpetra::BlockCrsMatrix. Sorry!");
724 W_ =
rcp (
new vector_type (A_->getRowMap ()));
725 W_->putScalar (STS::zero ());
730 for (
size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
732 w_ptr[LID] += STS::one();
740 if (schwarzCombineMode_ ==
"ADD") {
741 typedef Tpetra::MultiVector< typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type,typename MatrixType::node_type> scMV;
744 scMV nonOverLapW(theImport->getSourceMap(), 1,
false);
747 nonOverLapW.putScalar(STS::zero ());
748 for (
int ii = 0; ii < (int) theImport->getSourceMap()->getLocalNumElements(); ii++) nonOverLapWArray[ii] = w_ptr[ii];
749 nonOverLapWArray = Teuchos::null;
750 w_ptr = Teuchos::null;
751 nonOverLapW.doExport (*W_, *theImport, Tpetra::ADD);
752 W_->doImport( nonOverLapW, *theImport, Tpetra::INSERT);
756 W_->reciprocal (*W_);
760 InitializeTime_ += (timer->wallTime() - startTime);
762 IsInitialized_ =
true;
766 template<
class MatrixType,
class ContainerType>
774 (A_.
is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::compute: "
775 "The matrix is null. You must call setMatrix() with a nonnull matrix, "
776 "then call initialize(), before you may call this method.");
778 if (! isInitialized ()) {
785 double startTime = timer->
wallTime();
793 Container_->compute();
796 ComputeTime_ += (timer->
wallTime() - startTime);
801 template<
class MatrixType,
class ContainerType>
807 (Partitioner_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::"
808 "ExtractSubmatricesStructure: Partitioner object is null.");
810 std::string containerType = ContainerType::getName ();
811 if (containerType ==
"Generic") {
815 containerType = containerType_;
819 bool pointIndexed = decouple_ && hasBlockCrsMatrix_;
824 auto A_bcrs = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(A_);
825 local_ordinal_type dofs = hasBlockCrsMatrix_ ?
826 A_bcrs->getBlockSize() : List_.get<
int>(
"partitioner: PDE equations");
827 blockRows.
resize(NumLocalBlocks_ * dofs);
828 if(hasBlockCrsMatrix_)
830 for(local_ordinal_type i = 0; i < NumLocalBlocks_; i++)
832 size_t blockSize = Partitioner_->numRowsInPart(i);
835 for(local_ordinal_type j = 0; j < dofs; j++)
837 local_ordinal_type blockIndex = i * dofs + j;
838 blockRows[blockIndex].
resize(blockSize);
839 for(
size_t k = 0; k < blockSize; k++)
843 blockRows[blockIndex][k] = (*Partitioner_)(i, k) * dofs + j;
852 for(local_ordinal_type i = 0; i < NumLocalBlocks_; i++)
856 "Expected size of all blocks (partitions) to be divisible by MueLu dofs/node.");
857 size_t blockSize = Partitioner_->numRowsInPart(i) / dofs;
861 for(local_ordinal_type j = 0; j < dofs; j++)
863 local_ordinal_type blockIndex = i * dofs + j;
864 blockRows[blockIndex].
resize(blockSize);
865 for(
size_t k = 0; k < blockSize; k++)
867 blockRows[blockIndex][k] = (*Partitioner_)(i, k * dofs + j);
876 blockRows.
resize(NumLocalBlocks_);
877 for(local_ordinal_type i = 0; i < NumLocalBlocks_; ++i)
879 const size_t numRows = Partitioner_->numRowsInPart (i);
880 blockRows[i].
resize(numRows);
882 for (
size_t j = 0; j < numRows; ++j)
884 blockRows[i][j] = (*Partitioner_) (i,j);
890 Container_->setParameters(List_);
891 Container_->initialize();
894 template<
class MatrixType,
class ContainerType>
896 BlockRelaxation<MatrixType,ContainerType>::
897 ApplyInverseJacobi (
const MV& X, MV& Y)
const
899 const size_t NumVectors = X.getNumVectors ();
901 MV AY (Y.getMap (), NumVectors);
904 int starting_iteration = 0;
905 if (OverlapLevel_ > 0)
908 auto WView = W_->getLocalViewHost (Tpetra::Access::ReadOnly);
909 if(ZeroStartingSolution_) {
910 auto XView = X.getLocalViewHost (Tpetra::Access::ReadOnly);
911 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
912 Container_->DoOverlappingJacobi(XView, YView, WView, DampingFactor_, nonsymCombine_);
913 starting_iteration = 1;
915 const scalar_type ONE = STS::one();
916 for(
int j = starting_iteration; j < NumSweeps_; j++)
919 AY.update (ONE, X, -ONE);
921 auto AYView = AY.getLocalViewHost (Tpetra::Access::ReadOnly);
922 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
923 Container_->DoOverlappingJacobi (AYView, YView, WView, DampingFactor_, nonsymCombine_);
930 if(ZeroStartingSolution_)
932 auto XView = X.getLocalViewHost (Tpetra::Access::ReadOnly);
933 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
934 Container_->DoJacobi (XView, YView, DampingFactor_);
935 starting_iteration = 1;
937 const scalar_type ONE = STS::one();
938 for(
int j = starting_iteration; j < NumSweeps_; j++)
941 AY.update (ONE, X, -ONE);
943 auto AYView = AY.getLocalViewHost (Tpetra::Access::ReadOnly);
944 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
945 Container_->DoJacobi (AYView, YView, DampingFactor_);
951 template<
class MatrixType,
class ContainerType>
953 BlockRelaxation<MatrixType,ContainerType>::
954 ApplyInverseGS (
const MV& X, MV& Y)
const
958 size_t numVecs = X.getNumVectors();
960 auto XView = X.getLocalViewHost (Tpetra::Access::ReadOnly);
961 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
964 bool deleteY2 =
false;
967 Y2 = ptr(
new MV(Importer_->getTargetMap(), numVecs));
974 for(
int j = 0; j < NumSweeps_; ++j)
977 Y2->doImport(Y, *Importer_, Tpetra::INSERT);
978 auto Y2View = Y2->getLocalViewHost (Tpetra::Access::ReadWrite);
979 Container_->DoGaussSeidel(XView, YView, Y2View, DampingFactor_);
984 auto Y2View = Y2->getLocalViewHost (Tpetra::Access::ReadWrite);
985 for(
int j = 0; j < NumSweeps_; ++j)
987 Container_->DoGaussSeidel(XView, YView, Y2View, DampingFactor_);
994 template<
class MatrixType,
class ContainerType>
996 BlockRelaxation<MatrixType,ContainerType>::
997 ApplyInverseSGS (
const MV& X, MV& Y)
const
1002 auto XView = X.getLocalViewHost (Tpetra::Access::ReadOnly);
1003 auto YView = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
1006 bool deleteY2 =
false;
1009 Y2 = ptr(
new MV(Importer_->getTargetMap(), X.getNumVectors()));
1016 for(
int j = 0; j < NumSweeps_; ++j)
1019 Y2->doImport(Y, *Importer_, Tpetra::INSERT);
1020 auto Y2View = Y2->getLocalViewHost (Tpetra::Access::ReadWrite);
1021 Container_->DoSGS(XView, YView, Y2View, DampingFactor_);
1026 auto Y2View = Y2->getLocalViewHost (Tpetra::Access::ReadWrite);
1027 for(
int j = 0; j < NumSweeps_; ++j)
1029 Container_->DoSGS(XView, YView, Y2View, DampingFactor_);
1036 template<
class MatrixType,
class ContainerType>
1037 void BlockRelaxation<MatrixType,ContainerType>::computeImporter ()
const
1045 using Teuchos::rcp_dynamic_cast;
1048 if(hasBlockCrsMatrix_)
1050 const RCP<const block_crs_matrix_type> bcrs =
1051 rcp_dynamic_cast<
const block_crs_matrix_type>(A_);
1052 int bs_ = bcrs->getBlockSize();
1053 RCP<const map_type> oldDomainMap = A_->getDomainMap();
1054 RCP<const map_type> oldColMap = A_->getColMap();
1057 global_size_t numGlobalElements = oldColMap->getGlobalNumElements() * bs_;
1058 global_ordinal_type indexBase = oldColMap->getIndexBase();
1059 RCP<const Comm<int>> comm = oldColMap->getComm();
1060 ArrayView<const global_ordinal_type> oldColElements = oldColMap->getLocalElementList();
1061 Array<global_ordinal_type> newColElements(bs_ * oldColElements.size());
1062 for(
int i = 0; i < oldColElements.size(); i++)
1064 for(
int j = 0; j < bs_; j++)
1065 newColElements[i * bs_ + j] = oldColElements[i] * bs_ + j;
1067 RCP<map_type> colMap(
new map_type(numGlobalElements, newColElements, indexBase, comm));
1069 Importer_ =
rcp(
new import_type(oldDomainMap, colMap));
1073 Importer_ = A_->getGraph()->getImporter();
1074 if(Importer_.is_null())
1075 Importer_ =
rcp(
new import_type(A_->getDomainMap(), A_->getColMap()));
1081 template<
class MatrixType,
class ContainerType>
1086 std::ostringstream out;
1091 out <<
"\"Ifpack2::BlockRelaxation\": {";
1092 if (this->getObjectLabel () !=
"") {
1093 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
1098 out <<
"Matrix: null, ";
1103 out <<
"Global matrix dimensions: ["
1104 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"], ";
1109 out <<
"\"relaxation: type\": ";
1110 if (PrecType_ == Ifpack2::Details::JACOBI) {
1111 out <<
"Block Jacobi";
1112 }
else if (PrecType_ == Ifpack2::Details::GS) {
1113 out <<
"Block Gauss-Seidel";
1114 }
else if (PrecType_ == Ifpack2::Details::SGS) {
1115 out <<
"Block Symmetric Gauss-Seidel";
1116 }
else if (PrecType_ == Ifpack2::Details::MTSPLITJACOBI) {
1117 out <<
"MT Split Jacobi";
1123 if(hasBlockCrsMatrix_)
1127 int approx_rows_per_part = A_->getLocalNumRows()/Partitioner_->numLocalParts();
1128 out <<
", blocksize: "<<approx_rows_per_part;
1130 out <<
", overlap: " << OverlapLevel_;
1132 out <<
", " <<
"sweeps: " << NumSweeps_ <<
", "
1133 <<
"damping factor: " << DampingFactor_ <<
", ";
1135 std::string containerType = ContainerType::getName();
1136 out <<
"block container: " << ((containerType ==
"Generic") ? containerType_ : containerType);
1137 if(List_.isParameter(
"partitioner: PDE equations"))
1138 out <<
", dofs/node: "<<List_.get<
int>(
"partitioner: PDE equations");
1145 template<
class MatrixType,
class ContainerType>
1163 const int myImageID = A_->getComm()->getRank();
1173 if (vl !=
VERB_NONE && myImageID == 0) {
1174 out <<
"Ifpack2::BlockRelaxation<"
1175 << TypeNameTraits<MatrixType>::name () <<
", "
1176 << TypeNameTraits<ContainerType>::name () <<
" >:";
1179 if (this->getObjectLabel () !=
"") {
1180 out <<
"label: \"" << this->getObjectLabel () <<
"\"" << endl;
1182 out <<
"initialized: " << (isInitialized () ?
"true" :
"false") << endl
1183 <<
"computed: " << (isComputed () ?
"true" :
"false") << endl;
1185 std::string precType;
1186 if (PrecType_ == Ifpack2::Details::JACOBI) {
1187 precType =
"Block Jacobi";
1188 }
else if (PrecType_ == Ifpack2::Details::GS) {
1189 precType =
"Block Gauss-Seidel";
1190 }
else if (PrecType_ == Ifpack2::Details::SGS) {
1191 precType =
"Block Symmetric Gauss-Seidel";
1193 out <<
"type: " << precType << endl
1194 <<
"global number of rows: " << A_->getGlobalNumRows () << endl
1195 <<
"global number of columns: " << A_->getGlobalNumCols () << endl
1196 <<
"number of sweeps: " << NumSweeps_ << endl
1197 <<
"damping factor: " << DampingFactor_ << endl
1198 <<
"nonsymmetric overlap combine" << nonsymCombine_ << endl
1199 <<
"backwards mode: "
1200 << ((PrecType_ == Ifpack2::Details::GS && DoBackwardGS_) ?
"true" :
"false")
1202 <<
"zero starting solution: "
1203 << (ZeroStartingSolution_ ?
"true" :
"false") << endl;
1221 #ifdef HAVE_IFPACK2_EXPLICIT_INSTANTIATION
1223 #define IFPACK2_BLOCKRELAXATION_INSTANT(S,LO,GO,N) \
1225 class Ifpack2::BlockRelaxation< \
1226 Tpetra::RowMatrix<S, LO, GO, N>, \
1227 Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N> > >;
1229 #endif // HAVE_IFPACK2_EXPLICIT_INSTANTIATION
1231 #endif // IFPACK2_BLOCKRELAXATION_DEF_HPP
Ifpack2::BlockRelaxation class declaration.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the input matrix is distributed.
Definition: Ifpack2_BlockRelaxation_def.hpp:374
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
ParameterList & disableRecursiveValidation()
std::string description() const
A one-line description of this object.
Definition: Ifpack2_BlockRelaxation_def.hpp:1084
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:469
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:445
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
Teuchos::RCP< const map_type > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:411
int getNumCompute() const
Returns the number of calls to compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:437
void setParameters(const Teuchos::ParameterList ¶ms)
Sets all the parameters for the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:189
void compute()
compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation ...
Definition: Ifpack2_BlockRelaxation_def.hpp:769
Block relaxation preconditioners (or smoothers) for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse ma...
Definition: Ifpack2_BlockRelaxation_decl.hpp:82
double getInitializeTime() const
Returns the time spent in initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:453
Declaration of a user-defined partitioner in which the user can define a partition of the graph...
bool isParameter(const std::string &name) const
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
double getComputeTime() const
Returns the time spent in compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:461
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_BlockRelaxation_def.hpp:476
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_BlockRelaxation_def.hpp:62
virtual ~BlockRelaxation()
Destructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:125
Teuchos::RCP< const map_type > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:397
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 resize(size_type new_size, const value_type &x=value_type())
int getNumInitialize() const
Returns the number of calls to initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:430
void initialize()
Initialize.
Definition: Ifpack2_BlockRelaxation_def.hpp:612
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
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:1148
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:100
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:596
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:97
bool isType(const std::string &name) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Access only local rows and columns of a sparse matrix.
Definition: Ifpack2_LocalFilter_decl.hpp:161
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
Teuchos::RCP< const row_matrix_type > getMatrix() const
The input matrix of this preconditioner's constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:388
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:78
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:131
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:493
BlockRelaxation(const Teuchos::RCP< const row_matrix_type > &Matrix)
Constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:90