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> 55 template<
class MatrixType,
class ContainerType>
57 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
59 if (A.getRawPtr () != A_.getRawPtr ()) {
60 IsInitialized_ =
false;
62 Partitioner_ = Teuchos::null;
63 Importer_ = Teuchos::null;
65 hasBlockCrsMatrix_ =
false;
68 IsParallel_ = (A->getRowMap ()->getComm ()->getSize () > 1);
71 std::vector<Teuchos::RCP<ContainerType> > emptyVec;
72 std::swap (Containers_, emptyVec);
79 template<
class MatrixType,
class ContainerType>
83 Time_ (
Teuchos::rcp (new
Teuchos::Time (
"Ifpack2::BlockRelaxation"))),
85 PartitionerType_ (
"linear"),
89 DampingFactor_ (STS::one ()),
91 ZeroStartingSolution_ (true),
92 DoBackwardGS_ (false),
93 IsInitialized_ (false),
98 InitializeTime_ (0.0),
105 NumGlobalNonzeros_ (0),
106 hasBlockCrsMatrix_ (false)
109 template<
class MatrixType,
class ContainerType>
114 template<
class MatrixType,
class ContainerType>
119 Teuchos::ParameterList validparams;
121 List.validateParameters (validparams);
124 if (PrecType_ == Ifpack2::Details::JACOBI) {
126 }
else if (PrecType_ == Ifpack2::Details::GS) {
128 }
else if (PrecType_ == Ifpack2::Details::SGS) {
129 PT =
"Symmetric Gauss-Seidel";
134 if (PT ==
"Jacobi") {
135 PrecType_ = Ifpack2::Details::JACOBI;
137 else if (PT ==
"Gauss-Seidel") {
138 PrecType_ = Ifpack2::Details::GS;
140 else if (PT ==
"Symmetric Gauss-Seidel") {
141 PrecType_ = Ifpack2::Details::SGS;
144 TEUCHOS_TEST_FOR_EXCEPTION(
145 true, std::invalid_argument,
"Ifpack2::BlockRelaxation::setParameters: " 146 "Invalid parameter value \"" << PT <<
"\" for parameter \"relaxation: " 159 if (PrecType_ != Ifpack2::Details::JACOBI) {
162 if (NumLocalBlocks_ < 0) {
163 NumLocalBlocks_ = A_->getNodeNumRows() / (-NumLocalBlocks_);
168 TEUCHOS_TEST_FOR_EXCEPTION(
169 DoBackwardGS_, std::runtime_error,
170 "Ifpack2::BlockRelaxation:setParameters: Setting the \"relaxation: " 171 "backward mode\" parameter to true is not yet supported.");
178 template<
class MatrixType,
class ContainerType>
179 Teuchos::RCP<const Teuchos::Comm<int> >
182 TEUCHOS_TEST_FOR_EXCEPTION
183 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::getComm: " 184 "The matrix is null. You must call setMatrix() with a nonnull matrix " 185 "before you may call this method.");
186 return A_->getComm ();
189 template<
class MatrixType,
class ContainerType>
190 Teuchos::RCP<
const Tpetra::RowMatrix<
typename MatrixType::scalar_type,
191 typename MatrixType::local_ordinal_type,
192 typename MatrixType::global_ordinal_type,
193 typename MatrixType::node_type> >
198 template<
class MatrixType,
class ContainerType>
199 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
200 typename MatrixType::global_ordinal_type,
201 typename MatrixType::node_type> >
205 TEUCHOS_TEST_FOR_EXCEPTION
206 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::" 207 "getDomainMap: The matrix is null. You must call setMatrix() with a " 208 "nonnull matrix before you may call this method.");
209 return A_->getDomainMap ();
212 template<
class MatrixType,
class ContainerType>
213 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
214 typename MatrixType::global_ordinal_type,
215 typename MatrixType::node_type> >
219 TEUCHOS_TEST_FOR_EXCEPTION
220 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::" 221 "getRangeMap: The matrix is null. You must call setMatrix() with a " 222 "nonnull matrix before you may call this method.");
223 return A_->getRangeMap ();
226 template<
class MatrixType,
class ContainerType>
233 template<
class MatrixType,
class ContainerType>
237 return NumInitialize_;
240 template<
class MatrixType,
class ContainerType>
248 template<
class MatrixType,
class ContainerType>
256 template<
class MatrixType,
class ContainerType>
261 return InitializeTime_;
264 template<
class MatrixType,
class ContainerType>
272 template<
class MatrixType,
class ContainerType>
280 template<
class MatrixType,
class ContainerType>
285 return ComputeFlops_;
288 template<
class MatrixType,
class ContainerType>
296 template<
class MatrixType,
class ContainerType>
299 apply (
const Tpetra::MultiVector<
typename MatrixType::scalar_type,
300 typename MatrixType::local_ordinal_type,
301 typename MatrixType::global_ordinal_type,
302 typename MatrixType::node_type>& X,
303 Tpetra::MultiVector<
typename MatrixType::scalar_type,
304 typename MatrixType::local_ordinal_type,
305 typename MatrixType::global_ordinal_type,
306 typename MatrixType::node_type>& Y,
307 Teuchos::ETransp mode,
311 TEUCHOS_TEST_FOR_EXCEPTION
312 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::apply: " 313 "The matrix is null. You must call setMatrix() with a nonnull matrix, " 314 "then call initialize() and compute() (in that order), before you may " 315 "call this method.");
316 TEUCHOS_TEST_FOR_EXCEPTION(
317 !
isComputed (), std::runtime_error,
"Ifpack2::BlockRelaxation::apply: " 318 "isComputed() must be true prior to calling apply.");
319 TEUCHOS_TEST_FOR_EXCEPTION(
320 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
321 "Ifpack2::BlockRelaxation::apply: X.getNumVectors() = " 322 << X.getNumVectors () <<
" != Y.getNumVectors() = " 323 << Y.getNumVectors () <<
".");
324 TEUCHOS_TEST_FOR_EXCEPTION(
325 mode != Teuchos::NO_TRANS, std::logic_error,
"Ifpack2::BlockRelaxation::" 326 "apply: This method currently only implements the case mode == Teuchos::" 327 "NO_TRANS. Transposed modes are not currently supported.");
328 TEUCHOS_TEST_FOR_EXCEPTION(
329 alpha != Teuchos::ScalarTraits<scalar_type>::one (), std::logic_error,
330 "Ifpack2::BlockRelaxation::apply: This method currently only implements " 331 "the case alpha == 1. You specified alpha = " << alpha <<
".");
332 TEUCHOS_TEST_FOR_EXCEPTION(
333 beta != Teuchos::ScalarTraits<scalar_type>::zero (), std::logic_error,
334 "Ifpack2::BlockRelaxation::apply: This method currently only implements " 335 "the case beta == 0. You specified beta = " << beta <<
".");
341 Teuchos::RCP<const MV> X_copy;
343 auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
344 auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
345 if (X_lcl_host.ptr_on_device () == Y_lcl_host.ptr_on_device ()) {
346 X_copy = rcp (
new MV (X, Teuchos::Copy));
348 X_copy = rcpFromRef (X);
352 if (ZeroStartingSolution_) {
353 Y.putScalar (STS::zero ());
358 case Ifpack2::Details::JACOBI:
359 ApplyInverseJacobi(*X_copy,Y);
361 case Ifpack2::Details::GS:
362 ApplyInverseGS(*X_copy,Y);
364 case Ifpack2::Details::SGS:
365 ApplyInverseSGS(*X_copy,Y);
368 TEUCHOS_TEST_FOR_EXCEPTION
369 (
true, std::logic_error,
"Ifpack2::BlockRelaxation::apply: Invalid " 370 "PrecType_ enum value " << PrecType_ <<
". Valid values are Ifpack2::" 371 "Details::JACOBI = " << Ifpack2::Details::JACOBI <<
", Ifpack2::Details" 372 "::GS = " << Ifpack2::Details::GS <<
", and Ifpack2::Details::SGS = " 373 << Ifpack2::Details::SGS <<
". Please report this bug to the Ifpack2 " 379 ApplyTime_ += Time_->totalElapsedTime();
382 template<
class MatrixType,
class ContainerType>
385 applyMat (
const Tpetra::MultiVector<
typename MatrixType::scalar_type,
386 typename MatrixType::local_ordinal_type,
387 typename MatrixType::global_ordinal_type,
388 typename MatrixType::node_type>& X,
389 Tpetra::MultiVector<
typename MatrixType::scalar_type,
390 typename MatrixType::local_ordinal_type,
391 typename MatrixType::global_ordinal_type,
392 typename MatrixType::node_type>& Y,
393 Teuchos::ETransp mode)
const 395 A_->apply (X, Y, mode);
398 template<
class MatrixType,
class ContainerType>
404 typedef Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type>
407 TEUCHOS_TEST_FOR_EXCEPTION
408 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::initialize: " 409 "The matrix is null. You must call setMatrix() with a nonnull matrix " 410 "before you may call this method.");
413 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
414 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
415 if (A_bcrs.is_null ()) {
416 hasBlockCrsMatrix_ =
false;
419 hasBlockCrsMatrix_ =
true;
422 IsInitialized_ =
false;
425 NumMyRows_ = A_->getNodeNumRows ();
426 NumGlobalRows_ = A_->getGlobalNumRows ();
427 NumGlobalNonzeros_ = A_->getGlobalNumEntries ();
433 if (PartitionerType_ ==
"linear") {
436 }
else if (PartitionerType_ ==
"line") {
439 }
else if (PartitionerType_ ==
"user") {
445 TEUCHOS_TEST_FOR_EXCEPTION
446 (
true, std::logic_error,
"Ifpack2::BlockRelaxation::initialize: Unknown " 447 "partitioner type " << PartitionerType_ <<
". Valid values are " 448 "\"linear\", \"line\", and \"user\".");
452 Partitioner_->setParameters (List_);
453 Partitioner_->compute ();
456 NumLocalBlocks_ = Partitioner_->numLocalParts ();
461 if (A_->getComm()->getSize() != 1) {
469 InitializeTime_ += Time_->totalElapsedTime ();
470 IsInitialized_ =
true;
474 template<
class MatrixType,
class ContainerType>
479 typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
480 typedef Tpetra::Import<local_ordinal_type,global_ordinal_type, node_type> import_type;
483 using Teuchos::Array;
484 using Teuchos::ArrayView;
492 ExtractSubmatrices ();
496 TEUCHOS_TEST_FOR_EXCEPTION
497 (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0, std::runtime_error,
498 "Ifpack2::BlockRelaxation::computeBlockCrs: " 499 "We do not support overlapped Jacobi yet for Tpetra::BlockCrsMatrix. Sorry!");
507 if (IsParallel_ && (PrecType_ == Ifpack2::Details::GS ||
508 PrecType_ == Ifpack2::Details::SGS)) {
510 RCP<const block_crs_matrix_type> A_bcrs =
511 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
512 int bs = A_bcrs->getBlockSize();
516 RCP<const map_type> oldDomainMap = A_->getDomainMap();
517 RCP<const map_type> oldColMap = A_->getColMap();
521 global_size_t numGlobalElements = oldColMap->getGlobalNumElements()*bs;
523 RCP<const Teuchos::Comm<int> >comm = oldColMap->getComm();
524 ArrayView<const global_ordinal_type> oldColElements = oldColMap->getNodeElementList();
525 Array<global_ordinal_type> newColElements(bs*oldColElements.size());
527 for(
int i=0; i<oldColElements.size(); i++) {
528 for(
int j=0; j<bs; j++) {
529 newColElements[i*bs+j] = oldColElements[i]*bs+j;
532 RCP<map_type> colMap = rcp(
new map_type(numGlobalElements,newColElements,indexBase,comm));
535 Importer_ = rcp (
new import_type (oldDomainMap, colMap));
540 ComputeTime_ += Time_->totalElapsedTime();
545 template<
class MatrixType,
class ContainerType>
556 TEUCHOS_TEST_FOR_EXCEPTION
557 (A_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::compute: " 558 "The matrix is null. You must call setMatrix() with a nonnull matrix, " 559 "then call initialize(), before you may call this method.");
563 TEUCHOS_TEST_FOR_EXCEPTION
564 (NumSweeps_ < 0, std::logic_error,
"Ifpack2::BlockRelaxation::compute: " 565 "NumSweeps_ = " << NumSweeps_ <<
" < 0.");
571 if (hasBlockCrsMatrix_) {
582 ExtractSubmatrices ();
586 if (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0) {
588 W_ = rcp (
new vector_type (A_->getRowMap ()));
589 W_->putScalar (STS::zero ());
590 Teuchos::ArrayRCP<scalar_type > w_ptr = W_->getDataNonConst(0);
592 for (local_ordinal_type i = 0 ; i < NumLocalBlocks_ ; ++i) {
593 for (
size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
596 int LID = (*Partitioner_)(i,j);
597 w_ptr[LID]+= STS::one();
600 W_->reciprocal (*W_);
609 if (IsParallel_ && (PrecType_ == Ifpack2::Details::GS ||
610 PrecType_ == Ifpack2::Details::SGS)) {
611 Importer_ = A_->getGraph ()->getImporter ();
612 if (Importer_.is_null ()) {
613 Importer_ = rcp (
new import_type (A_->getDomainMap (), A_->getColMap ()));
619 ComputeTime_ += Time_->totalElapsedTime();
623 template<
class MatrixType,
class ContainerType>
628 TEUCHOS_TEST_FOR_EXCEPTION
629 (Partitioner_.is_null (), std::runtime_error,
"Ifpack2::BlockRelaxation::" 630 "ExtractSubmatrices: Partitioner object is null.");
632 NumLocalBlocks_ = Partitioner_->numLocalParts ();
633 Containers_.resize (NumLocalBlocks_);
636 const size_t numRows = Partitioner_->numRowsInPart (i);
639 Teuchos::Array<local_ordinal_type> localRows (numRows);
640 for (
size_t j = 0; j < numRows; ++j) {
641 localRows[j] = (*Partitioner_) (i,j);
643 if(numRows>1 || hasBlockCrsMatrix_) {
644 Containers_[i] = Teuchos::rcp (
new ContainerType (A_, localRows ()));
645 Containers_[i]->setParameters (List_);
646 Containers_[i]->initialize ();
647 Containers_[i]->compute ();
652 template<
class MatrixType,
class ContainerType>
657 const size_t NumVectors = X.getNumVectors ();
658 MV AY (Y.getMap (), NumVectors);
661 int starting_iteration = 0;
662 if (ZeroStartingSolution_) {
664 starting_iteration = 1;
668 for (
int j = starting_iteration; j < NumSweeps_; ++j) {
670 AY.update (ONE, X, -ONE);
674 ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
679 template<
class MatrixType,
class ContainerType>
684 const size_t NumVectors = X.getNumVectors ();
688 if (OverlapLevel_ == 0) {
692 if( Partitioner_->numRowsInPart (i) != 1 || hasBlockCrsMatrix_) {
693 if(Containers_[i]->getNumRows () == 0 )
continue;
694 Containers_[i]->apply (X, Y, Teuchos::NO_TRANS, DampingFactor_, one);
695 ApplyFlops_ += NumVectors * 2 * NumGlobalRows_;
700 Teuchos::ArrayRCP< const scalar_type > Diag = DiagRCP_->getData();
702 for(
unsigned int nv = 0;nv < NumVectors ; ++nv ) {
703 Teuchos::ArrayRCP< const scalar_type > xRCP = X.getData(nv);
705 Teuchos::ArrayRCP< scalar_type > yRCP = Y.getDataNonConst(nv);
717 if(Containers_[i]->getNumRows() == 0)
continue;
718 if ( Partitioner_->numRowsInPart (i) != 1 ) {
720 Containers_[i]->weightedApply(X,Y,*W_,Teuchos::NO_TRANS,DampingFactor_,one);
721 }
catch (std::exception& e) {
722 std::cerr <<
"BlockRelaxation::DoJacobi: Containers_[" << i
723 <<
"]->weightedApply() threw an exception: " << e.what ()
730 ApplyFlops_ += NumVectors * 4 * NumGlobalRows_;
735 template<
class MatrixType,
class ContainerType>
740 MV Xcopy (X, Teuchos::Copy);
741 for (
int j = 0; j < NumSweeps_; ++j) {
742 DoGaussSeidel (Xcopy, Y);
743 if (j != NumSweeps_ - 1) {
744 Tpetra::deep_copy (Xcopy, X);
749 template<
class MatrixType,
class ContainerType>
754 using Teuchos::Array;
755 using Teuchos::ArrayRCP;
756 using Teuchos::ArrayView;
759 using Teuchos::rcpFromRef;
764 const size_t Length = A_->getNodeMaxNumRowEntries();
765 const size_t NumVectors = X.getNumVectors();
766 Array<scalar_type> Values;
767 Array<local_ordinal_type> Indices;
768 Indices.resize (Length);
770 if(hasBlockCrsMatrix_)
772 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
773 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
774 int bs = A_bcrs->getBlockSize();
775 Values.resize (bs*bs*Length);
778 Values.resize (Length);
785 Y2 = rcp (
new MV (Importer_->getTargetMap (), NumVectors));
794 MV Residual (X.getMap (), NumVectors,
false);
796 ArrayRCP<ArrayRCP<scalar_type> > x_ptr = X.get2dViewNonConst();
797 ArrayRCP<ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();
798 ArrayRCP<ArrayRCP<scalar_type> > y2_ptr = Y2->get2dViewNonConst();
799 ArrayRCP<ArrayRCP<scalar_type> > residual_ptr = Residual.get2dViewNonConst();
802 if (IsParallel_) Y2->doImport(Y,*Importer_,Tpetra::INSERT);
805 if( Partitioner_->numRowsInPart (i) != 1 || hasBlockCrsMatrix_) {
806 if (Containers_[i]->getNumRows () == 0)
continue;
808 ArrayView<const local_ordinal_type> localRows =
809 Containers_[i]->getLocalRows ();
810 const size_t localNumRows = Containers_[i]->getNumRows ();
811 for (
size_t j = 0; j < localNumRows; ++j) {
814 A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
816 for (
size_t m = 0; m < NumVectors; ++m) {
817 ArrayView<const scalar_type> x_local = (x_ptr())[m]();
818 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
819 ArrayView<scalar_type> r_local = (residual_ptr())[m]();
821 if(hasBlockCrsMatrix_) {
822 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
823 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
824 int bs = A_bcrs->getBlockSize();
825 for (
int localR = 0; localR < bs; localR++)
826 r_local[LID*bs+localR] = x_local[LID*bs+localR];
827 for (
size_t k = 0; k < NumEntries; ++k) {
829 for (
int localR = 0; localR < bs; localR++) {
830 for(
int localC = 0; localC < bs; localC++) {
834 r_local[LID*bs+localR] -= Values[k*bs*bs+localR+localC*bs] * y2_local[col*bs+localC];
840 r_local[LID] = x_local[LID];
841 for (
size_t k = 0; k < NumEntries; ++k) {
843 r_local[LID] -= Values[k] * y2_local[col];
858 Containers_[i]->apply (Residual, *Y2, Teuchos::NO_TRANS,
864 ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
870 Teuchos::ArrayRCP< const scalar_type > Diag = DiagRCP_->getData();
872 ArrayRCP<ArrayRCP<scalar_type> > y2_ptr2 = Y2->get2dViewNonConst();
873 for(
unsigned int nv = 0;nv < NumVectors ; ++nv ) {
874 Teuchos::ArrayRCP< const scalar_type > xRCP = X.getData(nv);
876 ArrayView<scalar_type> y2_local = (y2_ptr2())[nv]();
878 y2_local[LRID]= newy;
885 if(hasBlockCrsMatrix_) {
886 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
887 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
888 int bs = A_bcrs->getBlockSize();
889 for (
size_t m = 0; m < NumVectors; ++m) {
890 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
891 ArrayView<scalar_type> y_local = (y_ptr())[m]();
892 for (
size_t i = 0; i < NumMyRows_*bs; ++i) {
893 y_local[i] = y2_local[i];
898 for (
size_t m = 0; m < NumVectors; ++m) {
899 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
900 ArrayView<scalar_type> y_local = (y_ptr())[m]();
901 for (
size_t i = 0; i < NumMyRows_; ++i) {
902 y_local[i] = y2_local[i];
909 template<
class MatrixType,
class ContainerType>
914 MV Xcopy (X, Teuchos::Copy);
915 for (
int j = 0; j < NumSweeps_; ++j) {
917 if (j != NumSweeps_ - 1) {
918 Tpetra::deep_copy (Xcopy, X);
923 template<
class MatrixType,
class ContainerType>
926 DoSGS (MV& X, MV& Y)
const 928 using Teuchos::Array;
929 using Teuchos::ArrayRCP;
930 using Teuchos::ArrayView;
933 using Teuchos::rcpFromRef;
936 const size_t Length = A_->getNodeMaxNumRowEntries();
937 const size_t NumVectors = X.getNumVectors();
938 Array<scalar_type> Values;
939 Array<local_ordinal_type> Indices;
940 Indices.resize(Length);
942 if(hasBlockCrsMatrix_)
944 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
945 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
946 int bs = A_bcrs->getBlockSize();
947 Values.resize (bs*bs*Length);
950 Values.resize (Length);
957 Y2 = rcp (
new MV (Importer_->getTargetMap (), NumVectors));
966 MV Residual (X.getMap (), NumVectors,
false);
968 ArrayRCP<ArrayRCP<scalar_type> > x_ptr = X.get2dViewNonConst();
969 ArrayRCP<ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();
970 ArrayRCP<ArrayRCP<scalar_type> > y2_ptr = Y2->get2dViewNonConst();
971 ArrayRCP<ArrayRCP<scalar_type> > residual_ptr = Residual.get2dViewNonConst();
975 Y2->doImport (Y, *Importer_, Tpetra::INSERT);
980 if( Partitioner_->numRowsInPart (i) != 1 || hasBlockCrsMatrix_) {
981 if (Containers_[i]->getNumRows () == 0) {
985 ArrayView<const local_ordinal_type> localRows =
986 Containers_[i]->getLocalRows ();
987 for (
size_t j = 0; j < Containers_[i]->getNumRows (); ++j) {
990 A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
993 for (
size_t m = 0; m < NumVectors; ++m) {
994 ArrayView<const scalar_type> x_local = (x_ptr())[m]();
995 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
996 ArrayView<scalar_type> r_local = (residual_ptr())[m]();
998 if(hasBlockCrsMatrix_) {
999 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
1000 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
1001 int bs = A_bcrs->getBlockSize();
1002 for (
int localR = 0; localR < bs; localR++)
1003 r_local[LID*bs+localR] = x_local[LID*bs+localR];
1004 for (
size_t k = 0; k < NumEntries; ++k) {
1006 for (
int localR = 0; localR < bs; localR++) {
1007 for(
int localC = 0; localC < bs; localC++) {
1008 r_local[LID*bs+localR] -= Values[k*bs*bs+localR+localC*bs] * y2_local[col*bs+localC];
1014 r_local[LID] = x_local[LID];
1015 for (
size_t k = 0 ; k < NumEntries ; k++) {
1017 r_local[LID] -= Values[k] * y2_local[col];
1028 Containers_[i]->apply (Residual, *Y2, Teuchos::NO_TRANS,
1029 DampingFactor_, one);
1032 ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
1038 Teuchos::ArrayRCP< const scalar_type > Diag = DiagRCP_->getData();
1040 for(
unsigned int nv = 0;nv < NumVectors ; ++nv ) {
1041 Teuchos::ArrayRCP< const scalar_type > xRCP = X.getData(nv);
1043 Teuchos::ArrayRCP< scalar_type > yRCP = Y.getDataNonConst(nv);
1059 if( hasBlockCrsMatrix_ || Partitioner_->numRowsInPart (i) != 1 ) {
1060 if (Containers_[i-1]->getNumRows () == 0)
continue;
1063 ArrayView<const local_ordinal_type> localRows =
1064 Containers_[i-1]->getLocalRows ();
1065 for (
size_t j = 0; j < Containers_[i-1]->getNumRows (); ++j) {
1068 A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
1071 for (
size_t m = 0; m < NumVectors; ++m) {
1072 ArrayView<const scalar_type> x_local = (x_ptr())[m]();
1073 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
1074 ArrayView<scalar_type> r_local = (residual_ptr())[m]();
1076 if(hasBlockCrsMatrix_) {
1077 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
1078 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
1079 int bs = A_bcrs->getBlockSize();
1080 for (
int localR = 0; localR < bs; localR++)
1081 r_local[LID*bs+localR] = x_local[LID*bs+localR];
1082 for (
size_t k = 0; k < NumEntries; ++k) {
1084 for (
int localR = 0; localR < bs; localR++) {
1085 for(
int localC = 0; localC < bs; localC++) {
1086 r_local[LID*bs+localR] -= Values[k*bs*bs+localR+localC*bs] * y2_local[col*bs+localC];
1092 r_local [LID] = x_local[LID];
1093 for (
size_t k = 0; k < NumEntries; ++k) {
1095 r_local[LID] -= Values[k] * y2_local[col];
1107 Containers_[i-1]->apply (Residual, *Y2, Teuchos::NO_TRANS,
1108 DampingFactor_, one);
1111 ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
1119 if(hasBlockCrsMatrix_) {
1120 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
1121 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
1122 int bs = A_bcrs->getBlockSize();
1123 for (
size_t m = 0; m < NumVectors; ++m) {
1124 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
1125 ArrayView<scalar_type> y_local = (y_ptr())[m]();
1126 for (
size_t i = 0; i < NumMyRows_*bs; ++i) {
1127 y_local[i] = y2_local[i];
1132 for (
size_t m = 0; m < NumVectors; ++m) {
1133 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
1134 ArrayView<scalar_type> y_local = (y_ptr())[m]();
1135 for (
size_t i = 0; i < NumMyRows_; ++i) {
1136 y_local[i] = y2_local[i];
1143 template<
class MatrixType,
class ContainerType>
1148 std::ostringstream out;
1153 out <<
"\"Ifpack2::BlockRelaxation\": {";
1154 if (this->getObjectLabel () !=
"") {
1155 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
1157 out <<
"Initialized: " << (
isInitialized () ?
"true" :
"false") <<
", ";
1158 out <<
"Computed: " << (
isComputed () ?
"true" :
"false") <<
", ";
1160 if (A_.is_null ()) {
1161 out <<
"Matrix: null, ";
1164 out <<
"Matrix: not null" 1165 <<
", Global matrix dimensions: [" 1166 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"], ";
1171 out <<
"\"relaxation: type\": ";
1172 if (PrecType_ == Ifpack2::Details::JACOBI) {
1173 out <<
"Block Jacobi";
1174 }
else if (PrecType_ == Ifpack2::Details::GS) {
1175 out <<
"Block Gauss-Seidel";
1176 }
else if (PrecType_ == Ifpack2::Details::SGS) {
1177 out <<
"Block Symmetric Gauss-Seidel";
1182 out <<
", " <<
"sweeps: " << NumSweeps_ <<
", " 1183 <<
"damping factor: " << DampingFactor_ <<
", ";
1189 template<
class MatrixType,
class ContainerType>
1193 const Teuchos::EVerbosityLevel verbLevel)
const 1197 using Teuchos::TypeNameTraits;
1198 using Teuchos::VERB_DEFAULT;
1199 using Teuchos::VERB_NONE;
1200 using Teuchos::VERB_LOW;
1201 using Teuchos::VERB_MEDIUM;
1202 using Teuchos::VERB_HIGH;
1203 using Teuchos::VERB_EXTREME;
1205 Teuchos::EVerbosityLevel vl = verbLevel;
1206 if (vl == VERB_DEFAULT) vl = VERB_LOW;
1207 const int myImageID = A_->getComm()->getRank();
1210 Teuchos::OSTab tab (out);
1217 if (vl != VERB_NONE && myImageID == 0) {
1218 out <<
"Ifpack2::BlockRelaxation<" 1219 << TypeNameTraits<MatrixType>::name () <<
", " 1220 << TypeNameTraits<ContainerType>::name () <<
" >:";
1221 Teuchos::OSTab tab1 (out);
1223 if (this->getObjectLabel () !=
"") {
1224 out <<
"label: \"" << this->getObjectLabel () <<
"\"" << endl;
1226 out <<
"initialized: " << (
isInitialized () ?
"true" :
"false") << endl
1227 <<
"computed: " << (
isComputed () ?
"true" :
"false") << endl;
1229 std::string precType;
1230 if (PrecType_ == Ifpack2::Details::JACOBI) {
1231 precType =
"Block Jacobi";
1232 }
else if (PrecType_ == Ifpack2::Details::GS) {
1233 precType =
"Block Gauss-Seidel";
1234 }
else if (PrecType_ == Ifpack2::Details::SGS) {
1235 precType =
"Block Symmetric Gauss-Seidel";
1237 out <<
"type: " << precType << endl
1238 <<
"global number of rows: " << A_->getGlobalNumRows () << endl
1239 <<
"global number of columns: " << A_->getGlobalNumCols () << endl
1240 <<
"number of sweeps: " << NumSweeps_ << endl
1241 <<
"damping factor: " << DampingFactor_ << endl
1242 <<
"backwards mode: " 1243 << ((PrecType_ == Ifpack2::Details::GS && DoBackwardGS_) ?
"true" :
"false")
1245 <<
"zero starting solution: " 1246 << (ZeroStartingSolution_ ?
"true" :
"false") << endl;
1248 out <<
"===============================================================================" << endl;
1249 out <<
"Phase # calls Total Time (s) Total MFlops MFlops/s " << endl;
1250 out <<
"------------ ------- --------------- --------------- ---------------" << endl;
1258 out <<
"===============================================================================" << endl;
1264 template<
class MatrixType,
class ContainerType>
1269 if(DiagRCP_ == Teuchos::null) {
1270 DiagRCP_ = Teuchos::rcp(
new vector_type(A_->getDomainMap ()));
1271 A_->getLocalDiagCopy (*DiagRCP_);
1278 #ifdef HAVE_IFPACK2_EXPLICIT_INSTANTIATION 1291 #define IFPACK2_BLOCKRELAXATION_INSTANT(S,LO,GO,N) \ 1293 class Ifpack2::BlockRelaxation< \ 1294 Tpetra::RowMatrix<S, LO, GO, N>, \ 1295 Ifpack2::SparseContainer< \ 1296 Tpetra::RowMatrix<S, LO, GO, N>, \ 1297 Ifpack2::ILUT< ::Tpetra::RowMatrix<S,LO,GO,N> > > >; \ 1299 class Ifpack2::BlockRelaxation< \ 1300 Tpetra::RowMatrix<S, LO, GO, N>, \ 1301 Ifpack2::DenseContainer< \ 1302 Tpetra::RowMatrix<S, LO, GO, N>, \ 1305 class Ifpack2::BlockRelaxation< \ 1306 Tpetra::RowMatrix<S, LO, GO, N>, \ 1307 Ifpack2::TriDiContainer< \ 1308 Tpetra::RowMatrix<S, LO, GO, N>, \ 1311 class Ifpack2::BlockRelaxation< \ 1312 Tpetra::RowMatrix<S, LO, GO, N>, \ 1313 Ifpack2::BandedContainer< \ 1314 Tpetra::RowMatrix<S, LO, GO, N>, \ 1317 #endif // HAVE_IFPACK2_EXPLICIT_INSTANTIATION 1319 #endif // IFPACK2_BLOCKRELAXATION_DEF_HPP double getApplyFlops() const
Returns the number of flops for the application of the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:291
void applyMat(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS) const
Applies the matrix to a Tpetra::MultiVector.
Definition: Ifpack2_BlockRelaxation_def.hpp:385
Ifpack2::BlockRelaxation class declaration.
Ifpack2::TriDiContainer class declaration.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the input matrix is distributed.
Definition: Ifpack2_BlockRelaxation_def.hpp:180
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:102
std::string description() const
A one-line description of this object.
Definition: Ifpack2_BlockRelaxation_def.hpp:1146
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:275
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:251
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 Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:217
int getNumCompute() const
Returns the number of calls to compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:243
void setParameters(const Teuchos::ParameterList ¶ms)
Sets all the parameters for the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:117
void compute()
compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation ...
Definition: Ifpack2_BlockRelaxation_def.hpp:548
Block relaxation preconditioners (or smoothers) for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse ma...
Definition: Ifpack2_BlockRelaxation_decl.hpp:81
Ifpack2 implementation details.
double getInitializeTime() const
Returns the time spent in initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:259
Declaration of a user-defined partitioner in which the user can define a nonoverlapping partition of ...
double getComputeFlops() const
Returns the number of flops in the computation phase.
Definition: Ifpack2_BlockRelaxation_def.hpp:283
double getComputeTime() const
Returns the time spent in compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:267
bool isComputed() const
Return true if compute() has been called.
Definition: Ifpack2_BlockRelaxation_decl.hpp:202
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_BlockRelaxation_def.hpp:57
virtual ~BlockRelaxation()
Destructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:111
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:203
Partition in which the user can define a nonoverlapping partition of the graph in any way they choose...
Definition: Ifpack2_Details_UserPartitioner_decl.hpp:66
bool isInitialized() const
Returns true if the preconditioner has been successfully initialized.
Definition: Ifpack2_BlockRelaxation_decl.hpp:192
int getNumInitialize() const
Returns the number of calls to initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:236
void initialize()
Initialize.
Definition: Ifpack2_BlockRelaxation_def.hpp:401
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:1192
Ifpack2::BandedContainer class declaration.
Ifpack2::DenseContainer class declaration.
MatrixType::node_type node_type
Node type of the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:105
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:99
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:63
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:96
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
Declaration of ILUT preconditioner.
Teuchos::RCP< const Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > getMatrix() const
The input matrix of this preconditioner's constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:194
void getValidParameters(Teuchos::ParameterList ¶ms)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:50
Ifpack2::SparseContainer class declaration.
Definition: Ifpack2_LinePartitioner_decl.hpp:75
A class to define linear partitions.
Definition: Ifpack2_LinearPartitioner_decl.hpp:54
void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &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:299
BlockRelaxation(const Teuchos::RCP< const row_matrix_type > &Matrix)
Constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:81