43 #ifndef IFPACK2_RELAXATION_DEF_HPP 44 #define IFPACK2_RELAXATION_DEF_HPP 46 #include "Teuchos_StandardParameterEntryValidators.hpp" 48 #include "Tpetra_CrsMatrix.hpp" 49 #include "Tpetra_Experimental_BlockCrsMatrix.hpp" 51 #include "Ifpack2_Relaxation_decl.hpp" 53 #ifdef HAVE_IFPACK2_EXPERIMENTAL_KOKKOSKERNELS_FEATURES 54 # include "KokkosKernels_GaussSeidel.hpp" 55 #endif // HAVE_IFPACK2_EXPERIMENTAL_KOKKOSKERNELS_FEATURES 57 #ifdef HAVE_IFPACK2_DUMP_MTX_MATRIX 58 # include "MatrixMarket_Tpetra.hpp" 72 NonnegativeIntValidator () {}
82 const std::string& paramName,
83 const std::string& sublistName)
const 90 anyVal.
type () !=
typeid (int),
92 "Parameter \"" << paramName <<
"\" in sublist \"" << sublistName
93 <<
"\" has the wrong type." << endl <<
"Parameter: " << paramName
94 << endl <<
"Type specified: " << entryName << endl
95 <<
"Type required: int" << endl);
97 const int val = Teuchos::any_cast<
int> (anyVal);
100 "Parameter \"" << paramName <<
"\" in sublist \"" << sublistName
101 <<
"\" is negative." << endl <<
"Parameter: " << paramName
102 << endl <<
"Value specified: " << val << endl
103 <<
"Required range: [0, INT_MAX]" << endl);
108 return "NonnegativeIntValidator";
113 printDoc (
const std::string& docString,
114 std::ostream &out)
const 117 out <<
"#\tValidator Used: " << std::endl;
118 out <<
"#\t\tNonnegativeIntValidator" << std::endl;
125 template<class Scalar, const bool isOrdinal=Teuchos::ScalarTraits<Scalar>::isOrdinal>
129 static const Scalar eps ();
133 template<
class Scalar>
134 class SmallTraits<Scalar, true> {
136 static const Scalar eps () {
142 template<
class Scalar>
143 class SmallTraits<Scalar, false> {
145 static const Scalar eps () {
153 template<
class MatrixType>
158 Importer_ = Teuchos::null;
159 Diagonal_ = Teuchos::null;
160 isInitialized_ =
false;
162 diagOffsets_ = Kokkos::View<size_t*, typename node_type::device_type> ();
163 savedDiagOffsets_ =
false;
164 hasBlockCrsMatrix_ =
false;
166 IsParallel_ = (A->getRowMap ()->getComm ()->getSize () > 1);
173 template<
class MatrixType>
180 DampingFactor_ (
STS::one ()),
181 IsParallel_ ((A.
is_null () || A->getRowMap ().
is_null () || A->getRowMap ()->getComm ().
is_null ()) ?
183 A->getRowMap ()->getComm ()->getSize () > 1),
184 ZeroStartingSolution_ (true),
185 DoBackwardGS_ (false),
188 MinDiagonalValue_ (
STS::
zero ()),
189 fixTinyDiagEntries_ (false),
190 checkDiagEntries_ (false),
191 isInitialized_ (false),
196 InitializeTime_ (0.0),
201 globalMinMagDiagEntryMag_ (
STM::
zero ()),
202 globalMaxMagDiagEntryMag_ (
STM::
zero ()),
203 globalNumSmallDiagEntries_ (0),
204 globalNumZeroDiagEntries_ (0),
205 globalNumNegDiagEntries_ (0),
207 savedDiagOffsets_ (false),
208 hasBlockCrsMatrix_ (false)
210 this->setObjectLabel (
"Ifpack2::Relaxation");
214 template<
class MatrixType>
218 template<
class MatrixType>
224 using Teuchos::parameterList;
227 using Teuchos::rcp_const_cast;
228 using Teuchos::rcp_implicit_cast;
229 using Teuchos::setStringToIntegralParameter;
232 if (validParams_.is_null ()) {
233 RCP<ParameterList> pl = parameterList (
"Ifpack2::Relaxation");
238 precTypes[0] =
"Jacobi";
239 precTypes[1] =
"Gauss-Seidel";
240 precTypes[2] =
"Symmetric Gauss-Seidel";
241 precTypes[3] =
"MT Gauss-Seidel";
242 precTypes[4] =
"MT Symmetric Gauss-Seidel";
243 Array<Details::RelaxationType> precTypeEnums (5);
244 precTypeEnums[0] = Details::JACOBI;
245 precTypeEnums[1] = Details::GS;
246 precTypeEnums[2] = Details::SGS;
247 precTypeEnums[3] = Details::MTGS;
248 precTypeEnums[4] = Details::MTSGS;
249 const std::string defaultPrecType (
"Jacobi");
250 setStringToIntegralParameter<Details::RelaxationType> (
"relaxation: type",
251 defaultPrecType,
"Relaxation method", precTypes (), precTypeEnums (),
254 const int numSweeps = 1;
255 RCP<PEV> numSweepsValidator =
256 rcp_implicit_cast<PEV> (
rcp (
new NonnegativeIntValidator));
257 pl->set (
"relaxation: sweeps", numSweeps,
"Number of relaxation sweeps",
258 rcp_const_cast<const PEV> (numSweepsValidator));
261 pl->set (
"relaxation: damping factor", dampingFactor);
263 const bool zeroStartingSolution =
true;
264 pl->set (
"relaxation: zero starting solution", zeroStartingSolution);
266 const bool doBackwardGS =
false;
267 pl->set (
"relaxation: backward mode", doBackwardGS);
269 const bool doL1Method =
false;
270 pl->set (
"relaxation: use l1", doL1Method);
272 const magnitude_type l1eta = (STM::one() + STM::one() + STM::one()) /
273 (STM::one() + STM::one());
274 pl->set (
"relaxation: l1 eta", l1eta);
277 pl->set (
"relaxation: min diagonal value", minDiagonalValue);
279 const bool fixTinyDiagEntries =
false;
280 pl->set (
"relaxation: fix tiny diagonal entries", fixTinyDiagEntries);
282 const bool checkDiagEntries =
false;
283 pl->set (
"relaxation: check diagonal entries", checkDiagEntries);
286 pl->set(
"relaxation: local smoothing indices", localSmoothingIndices);
288 validParams_ = rcp_const_cast<
const ParameterList> (pl);
294 template<
class MatrixType>
297 using Teuchos::getIntegralValue;
300 typedef scalar_type ST;
304 const Details::RelaxationType precType =
305 getIntegralValue<Details::RelaxationType> (pl,
"relaxation: type");
306 const int numSweeps = pl.get<
int> (
"relaxation: sweeps");
307 const ST dampingFactor = pl.get<ST> (
"relaxation: damping factor");
308 const bool zeroStartSol = pl.get<
bool> (
"relaxation: zero starting solution");
309 const bool doBackwardGS = pl.get<
bool> (
"relaxation: backward mode");
310 const bool doL1Method = pl.get<
bool> (
"relaxation: use l1");
311 const magnitude_type l1Eta = pl.get<magnitude_type> (
"relaxation: l1 eta");
312 const ST minDiagonalValue = pl.get<ST> (
"relaxation: min diagonal value");
313 const bool fixTinyDiagEntries = pl.get<
bool> (
"relaxation: fix tiny diagonal entries");
314 const bool checkDiagEntries = pl.get<
bool> (
"relaxation: check diagonal entries");
319 PrecType_ = precType;
320 NumSweeps_ = numSweeps;
321 DampingFactor_ = dampingFactor;
322 ZeroStartingSolution_ = zeroStartSol;
323 DoBackwardGS_ = doBackwardGS;
324 DoL1Method_ = doL1Method;
326 MinDiagonalValue_ = minDiagonalValue;
327 fixTinyDiagEntries_ = fixTinyDiagEntries;
328 checkDiagEntries_ = checkDiagEntries;
329 localSmoothingIndices_ = localSmoothingIndices;
333 template<
class MatrixType>
338 this->setParametersImpl(const_cast<Teuchos::ParameterList&>(pl));
342 template<
class MatrixType>
346 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::getComm: " 347 "The input matrix A is null. Please call setMatrix() with a nonnull " 348 "input matrix before calling this method.");
349 return A_->getRowMap ()->getComm ();
353 template<
class MatrixType>
360 template<
class MatrixType>
361 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
362 typename MatrixType::global_ordinal_type,
363 typename MatrixType::node_type> >
366 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::getDomainMap: " 367 "The input matrix A is null. Please call setMatrix() with a nonnull " 368 "input matrix before calling this method.");
369 return A_->getDomainMap ();
373 template<
class MatrixType>
374 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
375 typename MatrixType::global_ordinal_type,
376 typename MatrixType::node_type> >
379 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::getRangeMap: " 380 "The input matrix A is null. Please call setMatrix() with a nonnull " 381 "input matrix before calling this method.");
382 return A_->getRangeMap ();
386 template<
class MatrixType>
392 template<
class MatrixType>
394 return(NumInitialize_);
398 template<
class MatrixType>
404 template<
class MatrixType>
410 template<
class MatrixType>
412 return(InitializeTime_);
416 template<
class MatrixType>
418 return(ComputeTime_);
422 template<
class MatrixType>
428 template<
class MatrixType>
430 return(ComputeFlops_);
434 template<
class MatrixType>
440 template<
class MatrixType>
443 apply (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
444 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
452 using Teuchos::rcpFromRef;
456 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::apply: " 457 "The input matrix A is null. Please call setMatrix() with a nonnull " 458 "input matrix, then call compute(), before calling this method.");
462 "Ifpack2::Relaxation::apply: You must call compute() on this Ifpack2 " 463 "preconditioner instance before you may call apply(). You may call " 464 "isComputed() to find out if compute() has been called already.");
466 X.getNumVectors() != Y.getNumVectors(),
468 "Ifpack2::Relaxation::apply: X and Y have different numbers of columns. " 469 "X has " << X.getNumVectors() <<
" columns, but Y has " 470 << Y.getNumVectors() <<
" columns.");
472 beta != STS::zero (), std::logic_error,
473 "Ifpack2::Relaxation::apply: beta = " << beta <<
" != 0 case not " 481 if (alpha == STS::zero ()) {
483 Y.putScalar (STS::zero ());
491 auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
492 auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
493 if (X_lcl_host.ptr_on_device () == Y_lcl_host.ptr_on_device ()) {
494 Xcopy =
rcp (
new MV (X, Teuchos::Copy));
496 Xcopy = rcpFromRef (X);
504 case Ifpack2::Details::JACOBI:
505 ApplyInverseJacobi(*Xcopy,Y);
507 case Ifpack2::Details::GS:
508 ApplyInverseGS(*Xcopy,Y);
510 case Ifpack2::Details::SGS:
511 ApplyInverseSGS(*Xcopy,Y);
513 case Ifpack2::Details::MTSGS:
514 ApplyInverseMTSGS_CrsMatrix(*Xcopy,Y);
516 case Ifpack2::Details::MTGS:
517 ApplyInverseMTGS_CrsMatrix(*Xcopy,Y);
522 "Ifpack2::Relaxation::apply: Invalid preconditioner type enum value " 523 << PrecType_ <<
". Please report this bug to the Ifpack2 developers.");
525 if (alpha != STS::one ()) {
527 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
528 const double numVectors = as<double> (Y.getNumVectors ());
529 ApplyFlops_ += numGlobalRows * numVectors;
533 ApplyTime_ += Time_->totalElapsedTime ();
538 template<
class MatrixType>
541 applyMat (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
542 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
546 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::applyMat: " 547 "The input matrix A is null. Please call setMatrix() with a nonnull " 548 "input matrix, then call compute(), before calling this method.");
550 ! isComputed (), std::runtime_error,
"Ifpack2::Relaxation::applyMat: " 551 "isComputed() must be true before you may call applyMat(). " 552 "Please call compute() before calling this method.");
554 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
555 "Ifpack2::Relaxation::applyMat: X.getNumVectors() = " << X.getNumVectors ()
556 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
557 A_->apply (X, Y, mode);
561 template<
class MatrixType>
566 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::initialize: " 567 "The input matrix A is null. Please call setMatrix() with a nonnull " 568 "input matrix before calling this method.");
573 hasBlockCrsMatrix_ =
false;
577 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
579 hasBlockCrsMatrix_ =
false;
582 hasBlockCrsMatrix_ =
true;
587 #ifdef HAVE_IFPACK2_EXPERIMENTAL_KOKKOSKERNELS_FEATURES 589 if (PrecType_ == Ifpack2::Details::MTGS || PrecType_ == Ifpack2::Details::MTSGS) {
590 const crs_matrix_type* crsMat =
dynamic_cast<const crs_matrix_type*
> (&(*A_));
592 crsMat == NULL, std::runtime_error,
"Ifpack2::Relaxation::compute: " 593 "MT methods works for CRSMatrix Only.");
595 #ifdef HAVE_IFPACK2_DUMP_MTX_MATRIX 596 Tpetra::MatrixMarket::Writer<crs_matrix_type> crs_writer;
597 std::string file_name =
"Ifpack2_MT_GS.mtx";
599 crs_writer.writeSparseFile(file_name, rcp_crs_mat);
602 this->kh = Teuchos::rcp(
new KernelHandle());
603 if (kh->get_gs_handle() == NULL){
604 kh->create_gs_handle();
606 kokkos_csr_matrix kcsr = crsMat->getLocalMatrix ();
608 bool is_symmetric =
false;
609 if (PrecType_ == Ifpack2::Details::MTSGS){
612 KokkosKernels::Experimental::Graph::gauss_seidel_symbolic
613 <KernelHandle, lno_row_view_t, lno_nonzero_view_t>
614 (kh.getRawPtr(), A_->getNodeNumRows(),
615 A_->getNodeNumCols(),
624 InitializeTime_ += Time_->totalElapsedTime ();
626 isInitialized_ =
true;
631 template <
typename BlockDiagView>
632 struct InvertDiagBlocks {
633 typedef int value_type;
634 typedef typename BlockDiagView::size_type Size;
637 typedef Kokkos::MemoryTraits<Kokkos::Unmanaged> Unmanaged;
638 template <
typename View>
639 using UnmanagedView = Kokkos::View<
typename View::data_type,
typename View::array_layout,
640 typename View::device_type, Unmanaged>;
642 typedef typename BlockDiagView::non_const_value_type Scalar;
643 typedef typename BlockDiagView::device_type Device;
644 typedef Kokkos::View<Scalar**, Kokkos::LayoutRight, Device> RWrk;
645 typedef Kokkos::View<int**, Kokkos::LayoutRight, Device> IWrk;
647 UnmanagedView<BlockDiagView> block_diag_;
651 UnmanagedView<RWrk> rwrk_;
653 UnmanagedView<IWrk> iwrk_;
656 InvertDiagBlocks (BlockDiagView& block_diag)
657 : block_diag_(block_diag)
659 const auto blksz = block_diag.dimension_1();
660 Kokkos::resize(rwrk_buf_, block_diag_.dimension_0(), blksz);
662 Kokkos::resize(iwrk_buf_, block_diag_.dimension_0(), blksz);
666 KOKKOS_INLINE_FUNCTION
667 void operator() (
const Size i,
int& jinfo)
const {
668 auto D_cur = Kokkos::subview(block_diag_, i, Kokkos::ALL(), Kokkos::ALL());
669 auto ipiv = Kokkos::subview(iwrk_, i, Kokkos::ALL());
670 auto work = Kokkos::subview(rwrk_, i, Kokkos::ALL());
672 Tpetra::Experimental::GETF2(D_cur, ipiv, info);
677 Tpetra::Experimental::GETRI(D_cur, ipiv, work, info);
682 KOKKOS_INLINE_FUNCTION
683 void join (
volatile value_type& dst,
volatile value_type
const& src)
const { dst += src; }
687 template<
class MatrixType>
688 void Relaxation<MatrixType>::computeBlockCrs ()
698 using Teuchos::REDUCE_MAX;
699 using Teuchos::REDUCE_MIN;
700 using Teuchos::REDUCE_SUM;
701 using Teuchos::rcp_dynamic_cast;
703 typedef local_ordinal_type LO;
704 typedef typename node_type::device_type device_type;
712 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::" 713 "computeBlockCrs: The input matrix A is null. Please call setMatrix() " 714 "with a nonnull input matrix, then call initialize(), before calling " 716 const block_crs_matrix_type* blockCrsA =
717 dynamic_cast<const block_crs_matrix_type*
> (A_.getRawPtr ());
719 blockCrsA == NULL, std::logic_error,
"Ifpack2::Relaxation::" 720 "computeBlockCrs: A_ is not a BlockCrsMatrix, but it should be if we " 721 "got this far. Please report this bug to the Ifpack2 developers.");
723 const scalar_type one = STS::one ();
728 const LO lclNumMeshRows =
729 blockCrsA->getCrsGraph ().getNodeNumRows ();
730 const LO blockSize = blockCrsA->getBlockSize ();
732 if (! savedDiagOffsets_) {
733 blockDiag_ = block_diag_type ();
734 blockDiag_ = block_diag_type (
"Ifpack2::Relaxation::blockDiag_",
735 lclNumMeshRows, blockSize, blockSize);
736 if (Teuchos::as<LO>(diagOffsets_.dimension_0 () ) < lclNumMeshRows) {
739 diagOffsets_ = Kokkos::View<size_t*, device_type> ();
740 diagOffsets_ = Kokkos::View<size_t*, device_type> (
"offsets", lclNumMeshRows);
742 blockCrsA->getCrsGraph ().getLocalDiagOffsets (diagOffsets_);
744 (static_cast<size_t> (diagOffsets_.dimension_0 ()) !=
745 static_cast<size_t> (blockDiag_.dimension_0 ()),
746 std::logic_error,
"diagOffsets_.dimension_0() = " <<
747 diagOffsets_.dimension_0 () <<
" != blockDiag_.dimension_0() = " 748 << blockDiag_.dimension_0 () <<
749 ". Please report this bug to the Ifpack2 developers.");
750 savedDiagOffsets_ =
true;
752 blockCrsA->getLocalDiagCopy (blockDiag_, diagOffsets_);
759 unmanaged_block_diag_type blockDiag = blockDiag_;
761 if (DoL1Method_ && IsParallel_) {
762 const scalar_type two = one + one;
763 const size_t maxLength = A_->getNodeMaxNumRowEntries ();
764 Array<LO> indices (maxLength);
765 Array<scalar_type> values (maxLength * blockSize * blockSize);
766 size_t numEntries = 0;
768 for (LO i = 0; i < lclNumMeshRows; ++i) {
770 blockCrsA->getLocalRowCopy (i, indices (), values (), numEntries);
772 auto diagBlock = Kokkos::subview (blockDiag, i, ALL (), ALL ());
773 for (LO subRow = 0; subRow < blockSize; ++subRow) {
774 magnitude_type diagonal_boost = STM::zero ();
775 for (
size_t k = 0 ; k < numEntries ; ++k) {
776 if (indices[k] > lclNumMeshRows) {
777 const size_t offset = blockSize*blockSize*k + subRow*blockSize;
778 for (LO subCol = 0; subCol < blockSize; ++subCol) {
779 diagonal_boost += STS::magnitude (values[offset+subCol] / two);
783 if (STS::magnitude (diagBlock(subRow, subRow)) < L1Eta_ * diagonal_boost) {
784 diagBlock(subRow, subRow) += diagonal_boost;
792 Impl::InvertDiagBlocks<unmanaged_block_diag_type> idb(blockDiag);
793 Kokkos::parallel_reduce(lclNumMeshRows, idb, info);
795 (info > 0, std::runtime_error,
"GETF2 or GETRI failed on = " << info
796 <<
" diagonal blocks.");
801 #ifdef HAVE_IFPACK2_DEBUG 802 const int numResults = 2;
804 int lclResults[2], gblResults[2];
805 lclResults[0] = info;
806 lclResults[1] = -info;
809 reduceAll<int, int> (* (A_->getGraph ()->getComm ()), REDUCE_MIN,
810 numResults, lclResults, gblResults);
812 (gblResults[0] != 0 || gblResults[1] != 0, std::runtime_error,
813 "Ifpack2::Relaxation::compute: When processing the input " 814 "Tpetra::BlockCrsMatrix, one or more diagonal block LU factorizations " 815 "failed on one or more (MPI) processes.");
816 #endif // HAVE_IFPACK2_DEBUG 818 Importer_ = A_->getGraph ()->getImporter ();
821 ComputeTime_ += Time_->totalElapsedTime ();
826 template<
class MatrixType>
836 using Teuchos::REDUCE_MAX;
837 using Teuchos::REDUCE_MIN;
838 using Teuchos::REDUCE_SUM;
839 using Teuchos::rcp_dynamic_cast;
840 using Teuchos::reduceAll;
843 typedef typename vector_type::device_type device_type;
848 if (! isInitialized ()) {
852 if (hasBlockCrsMatrix_) {
865 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::compute: " 866 "The input matrix A is null. Please call setMatrix() with a nonnull " 867 "input matrix, then call initialize(), before calling this method.");
873 NumSweeps_ < 0, std::logic_error,
874 "Ifpack2::Relaxation::compute: NumSweeps_ = " << NumSweeps_ <<
" < 0. " 875 "Please report this bug to the Ifpack2 developers.");
877 Diagonal_ =
rcp (
new vector_type (A_->getRowMap ()));
895 const crs_matrix_type* crsMat =
896 dynamic_cast<const crs_matrix_type*
> (A_.getRawPtr ());
897 if (crsMat == NULL || ! crsMat->isStaticGraph ()) {
898 A_->getLocalDiagCopy (*Diagonal_);
900 if (! savedDiagOffsets_) {
901 const size_t lclNumRows = A_->getRowMap ()->getNodeNumElements ();
902 if (diagOffsets_.dimension_0 () < lclNumRows) {
903 typedef typename node_type::device_type DT;
904 diagOffsets_ = Kokkos::View<size_t*, DT> ();
905 diagOffsets_ = Kokkos::View<size_t*, DT> (
"offsets", lclNumRows);
907 crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
908 savedDiagOffsets_ =
true;
910 crsMat->getLocalDiagCopy (*Diagonal_, diagOffsets_);
911 #ifdef HAVE_IFPACK2_DEBUG 913 vector_type D_copy (A_->getRowMap ());
914 A_->getLocalDiagCopy (D_copy);
915 D_copy.update (STS::one (), *Diagonal_, -STS::one ());
920 err != STM::zero(), std::logic_error,
"Ifpack2::Relaxation::compute: " 921 "\"fast-path\" diagonal computation failed. \\|D1 - D2\\|_inf = " 923 #endif // HAVE_IFPACK2_DEBUG 929 RCP<vector_type> origDiag;
930 if (checkDiagEntries_) {
931 origDiag =
rcp (
new vector_type (A_->getRowMap ()));
932 Tpetra::deep_copy (*origDiag, *Diagonal_);
935 const size_t numMyRows = A_->getNodeNumRows ();
938 Diagonal_->template sync<Kokkos::HostSpace> ();
939 Diagonal_->template modify<Kokkos::HostSpace> ();
940 auto diag_2d = Diagonal_->template getLocalView<Kokkos::HostSpace> ();
941 auto diag_1d = Kokkos::subview (diag_2d, Kokkos::ALL (), 0);
952 if (DoL1Method_ && IsParallel_) {
954 const size_t maxLength = A_->getNodeMaxNumRowEntries ();
955 Array<local_ordinal_type> indices (maxLength);
956 Array<scalar_type> values (maxLength);
959 for (
size_t i = 0; i < numMyRows; ++i) {
960 A_->getLocalRowCopy (i, indices (), values (), numEntries);
962 for (
size_t k = 0 ; k < numEntries ; ++k) {
963 if (static_cast<size_t> (indices[k]) > numMyRows) {
964 diagonal_boost += STS::magnitude (values[k] / two);
967 if (STS::magnitude (diag[i]) < L1Eta_ * diagonal_boost) {
968 diag[i] += diagonal_boost;
985 one / SmallTraits<scalar_type>::eps () :
986 one / MinDiagonalValue_;
988 const magnitude_type minDiagValMag = STS::magnitude (MinDiagonalValue_);
990 if (checkDiagEntries_) {
996 size_t numSmallDiagEntries = 0;
997 size_t numZeroDiagEntries = 0;
998 size_t numNegDiagEntries = 0;
1009 if (numMyRows > 0) {
1014 minMagDiagEntryMag = d_0_mag;
1015 maxMagDiagEntryMag = d_0_mag;
1023 for (
size_t i = 0 ; i < numMyRows; ++i) {
1030 if (d_i_real < STM::zero ()) {
1031 ++numNegDiagEntries;
1033 if (d_i_mag < minMagDiagEntryMag) {
1035 minMagDiagEntryMag = d_i_mag;
1037 if (d_i_mag > maxMagDiagEntryMag) {
1039 maxMagDiagEntryMag = d_i_mag;
1042 if (fixTinyDiagEntries_) {
1044 if (d_i_mag <= minDiagValMag) {
1045 ++numSmallDiagEntries;
1046 if (d_i_mag == STM::zero ()) {
1047 ++numZeroDiagEntries;
1049 diag[i] = oneOverMinDiagVal;
1051 diag[i] = one / d_i;
1056 if (d_i_mag <= minDiagValMag) {
1057 ++numSmallDiagEntries;
1058 if (d_i_mag == STM::zero ()) {
1059 ++numZeroDiagEntries;
1062 diag[i] = one / d_i;
1069 if (STS::isComplex) {
1070 ComputeFlops_ += 4.0 * numMyRows;
1072 ComputeFlops_ += numMyRows;
1089 const Comm<int>& comm = * (A_->getRowMap ()->getComm ());
1092 Array<magnitude_type> localVals (2);
1093 localVals[0] = minMagDiagEntryMag;
1095 localVals[1] = -maxMagDiagEntryMag;
1096 Array<magnitude_type> globalVals (2);
1097 reduceAll<int, magnitude_type> (comm, REDUCE_MIN, 2,
1098 localVals.getRawPtr (),
1099 globalVals.getRawPtr ());
1100 globalMinMagDiagEntryMag_ = globalVals[0];
1101 globalMaxMagDiagEntryMag_ = -globalVals[1];
1103 Array<size_t> localCounts (3);
1104 localCounts[0] = numSmallDiagEntries;
1105 localCounts[1] = numZeroDiagEntries;
1106 localCounts[2] = numNegDiagEntries;
1107 Array<size_t> globalCounts (3);
1108 reduceAll<int, size_t> (comm, REDUCE_SUM, 3,
1109 localCounts.getRawPtr (),
1110 globalCounts.getRawPtr ());
1111 globalNumSmallDiagEntries_ = globalCounts[0];
1112 globalNumZeroDiagEntries_ = globalCounts[1];
1113 globalNumNegDiagEntries_ = globalCounts[2];
1125 vector_type diff (A_->getRowMap ());
1126 diff.reciprocal (*origDiag);
1127 Diagonal_->template sync<device_type> ();
1128 diff.update (-one, *Diagonal_, one);
1129 globalDiagNormDiff_ = diff.norm2 ();
1132 if (fixTinyDiagEntries_) {
1136 for (
size_t i = 0 ; i < numMyRows; ++i) {
1141 if (d_i_mag <= minDiagValMag) {
1142 diag[i] = oneOverMinDiagVal;
1144 diag[i] = one / d_i;
1149 for (
size_t i = 0 ; i < numMyRows; ++i) {
1150 diag[i] = one / diag[i];
1153 if (STS::isComplex) {
1154 ComputeFlops_ += 4.0 * numMyRows;
1156 ComputeFlops_ += numMyRows;
1160 if (IsParallel_ && (PrecType_ == Ifpack2::Details::GS ||
1161 PrecType_ == Ifpack2::Details::SGS)) {
1165 Importer_ = A_->getGraph ()->getImporter ();
1166 Diagonal_->template sync<device_type> ();
1168 #ifdef HAVE_IFPACK2_EXPERIMENTAL_KOKKOSKERNELS_FEATURES 1170 if (PrecType_ == Ifpack2::Details::MTGS || PrecType_ == Ifpack2::Details::MTSGS) {
1171 const crs_matrix_type* crsMat =
dynamic_cast<const crs_matrix_type*
> (&(*A_));
1173 crsMat == NULL, std::runtime_error,
"Ifpack2::Relaxation::compute: " 1174 "MT methods works for CRSMatrix Only.");
1176 kokkos_csr_matrix kcsr = crsMat->getLocalMatrix ();
1178 bool is_symmetric =
false;
1179 if (PrecType_ == Ifpack2::Details::MTSGS){
1180 is_symmetric =
true;
1182 KokkosKernels::Experimental::Graph::gauss_seidel_numeric
1183 <KernelHandle, lno_row_view_t, lno_nonzero_view_t, scalar_nonzero_view_t>
1184 (kh.getRawPtr(), A_->getNodeNumRows(), A_->getNodeNumCols(), kcsr.graph.row_map,
1185 kcsr.graph.entries, kcsr.values, is_symmetric);
1190 ComputeTime_ += Time_->totalElapsedTime ();
1196 template<
class MatrixType>
1199 ApplyInverseJacobi (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1200 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const 1203 typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
1204 global_ordinal_type, node_type> MV;
1206 if (hasBlockCrsMatrix_) {
1207 ApplyInverseJacobi_BlockCrsMatrix (X, Y);
1211 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1212 const double numVectors = as<double> (X.getNumVectors ());
1213 if (ZeroStartingSolution_) {
1218 Y.elementWiseMultiply (DampingFactor_, *Diagonal_, X, STS::zero ());
1224 double flopUpdate = 0.0;
1225 if (DampingFactor_ == STS::one ()) {
1227 flopUpdate = numGlobalRows * numVectors;
1231 flopUpdate = 2.0 * numGlobalRows * numVectors;
1233 ApplyFlops_ += flopUpdate;
1234 if (NumSweeps_ == 1) {
1240 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1243 MV A_times_Y (Y.getMap (), as<size_t>(numVectors),
false);
1244 for (
int j = startSweep; j < NumSweeps_; ++j) {
1246 applyMat (Y, A_times_Y);
1247 A_times_Y.update (STS::one (), X, -STS::one ());
1248 Y.elementWiseMultiply (DampingFactor_, *Diagonal_, A_times_Y, STS::one ());
1262 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1263 const double dampingFlops = (DampingFactor_ == STS::one ()) ? 0.0 : 1.0;
1264 ApplyFlops_ += as<double> (NumSweeps_ - startSweep) * numVectors *
1265 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1268 template<
class MatrixType>
1270 Relaxation<MatrixType>::
1271 ApplyInverseJacobi_BlockCrsMatrix (
const Tpetra::MultiVector<scalar_type,
1273 global_ordinal_type,
1275 Tpetra::MultiVector<scalar_type,
1277 global_ordinal_type,
1278 node_type>& Y)
const 1280 typedef Tpetra::Experimental::BlockMultiVector<scalar_type,
1281 local_ordinal_type, global_ordinal_type, node_type> BMV;
1283 const block_crs_matrix_type* blockMatConst =
1284 dynamic_cast<const block_crs_matrix_type*
> (A_.getRawPtr ());
1286 (blockMatConst == NULL, std::logic_error,
"This method should never be " 1287 "called if the matrix A_ is not a BlockCrsMatrix. Please report this " 1288 "bug to the Ifpack2 developers.");
1293 block_crs_matrix_type* blockMat =
1294 const_cast<block_crs_matrix_type*
> (blockMatConst);
1296 auto meshRowMap = blockMat->getRowMap ();
1297 auto meshColMap = blockMat->getColMap ();
1298 const local_ordinal_type blockSize = blockMat->getBlockSize ();
1300 BMV X_blk (X, *meshColMap, blockSize);
1301 BMV Y_blk (Y, *meshRowMap, blockSize);
1303 if (ZeroStartingSolution_) {
1307 Y_blk.blockWiseMultiply (DampingFactor_, blockDiag_, X_blk);
1308 if (NumSweeps_ == 1) {
1313 auto pointRowMap = Y.getMap ();
1314 const size_t numVecs = X.getNumVectors ();
1318 BMV A_times_Y (*meshRowMap, *pointRowMap, blockSize, numVecs);
1322 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1324 for (
int j = startSweep; j < NumSweeps_; ++j) {
1325 blockMat->applyBlock (Y_blk, A_times_Y);
1327 Y_blk.blockJacobiUpdate (DampingFactor_, blockDiag_,
1328 X_blk, A_times_Y, STS::one ());
1332 template<
class MatrixType>
1334 Relaxation<MatrixType>::
1335 ApplyInverseGS (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1336 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const 1338 typedef Relaxation<MatrixType> this_type;
1348 const block_crs_matrix_type* blockCrsMat =
1349 dynamic_cast<const block_crs_matrix_type*
> (A_.getRawPtr ());
1350 const crs_matrix_type* crsMat =
1351 dynamic_cast<const crs_matrix_type*
> (A_.getRawPtr ());
1352 if (blockCrsMat != NULL) {
1353 const_cast<this_type*
> (
this)->ApplyInverseGS_BlockCrsMatrix (*blockCrsMat, X, Y);
1354 }
else if (crsMat != NULL) {
1355 ApplyInverseGS_CrsMatrix (*crsMat, X, Y);
1357 ApplyInverseGS_RowMatrix (X, Y);
1362 template<
class MatrixType>
1364 Relaxation<MatrixType>::
1365 ApplyInverseGS_RowMatrix (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1366 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const 1374 using Teuchos::rcpFromRef;
1375 typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
1380 if (ZeroStartingSolution_) {
1381 Y.putScalar (STS::zero ());
1384 const size_t NumVectors = X.getNumVectors ();
1385 const size_t maxLength = A_->getNodeMaxNumRowEntries ();
1386 Array<local_ordinal_type> Indices (maxLength);
1387 Array<scalar_type> Values (maxLength);
1390 const size_t numMyRows = A_->getNodeNumRows();
1391 const local_ordinal_type* rowInd = 0;
1392 size_t numActive = numMyRows;
1393 bool do_local = ! localSmoothingIndices_.is_null ();
1395 rowInd = localSmoothingIndices_.getRawPtr ();
1396 numActive = localSmoothingIndices_.size ();
1401 if (Importer_.is_null ()) {
1403 Y2 =
rcp (
new MV (Y.getMap (), NumVectors,
false));
1408 Y2 =
rcp (
new MV (Importer_->getTargetMap (), NumVectors));
1412 Y2 = rcpFromRef (Y);
1416 ArrayRCP<const scalar_type> d_rcp = Diagonal_->get1dView ();
1417 ArrayView<const scalar_type> d_ptr = d_rcp();
1420 bool constant_stride = X.isConstantStride() && Y2->isConstantStride();
1422 if (constant_stride) {
1424 size_t x_stride = X.getStride();
1425 size_t y2_stride = Y2->getStride();
1426 ArrayRCP<scalar_type> y2_rcp = Y2->get1dViewNonConst();
1427 ArrayRCP<const scalar_type> x_rcp = X.get1dView();
1428 ArrayView<scalar_type> y2_ptr = y2_rcp();
1429 ArrayView<const scalar_type> x_ptr = x_rcp();
1430 Array<scalar_type> dtemp(NumVectors,STS::zero());
1432 for (
int j = 0; j < NumSweeps_; ++j) {
1435 if (Importer_.is_null ()) {
1438 Y2->doImport (Y, *Importer_, Tpetra::INSERT);
1442 if (! DoBackwardGS_) {
1443 for (
size_t ii = 0; ii < numActive; ++ii) {
1444 local_ordinal_type i = as<local_ordinal_type> (do_local ? rowInd[ii] : ii);
1446 A_->getLocalRowCopy (i, Indices (), Values (), NumEntries);
1447 dtemp.assign(NumVectors,STS::zero());
1449 for (
size_t k = 0; k < NumEntries; ++k) {
1450 const local_ordinal_type col = Indices[k];
1451 for (
size_t m = 0; m < NumVectors; ++m) {
1452 dtemp[m] += Values[k] * y2_ptr[col + y2_stride*m];
1456 for (
size_t m = 0; m < NumVectors; ++m) {
1457 y2_ptr[i + y2_stride*m] += DampingFactor_ * d_ptr[i] * (x_ptr[i + x_stride*m] - dtemp[m]);
1464 for (ptrdiff_t ii = as<ptrdiff_t> (numActive) - 1; ii >= 0; --ii) {
1465 local_ordinal_type i = as<local_ordinal_type> (do_local ? rowInd[ii] : ii);
1467 A_->getLocalRowCopy (i, Indices (), Values (), NumEntries);
1468 dtemp.assign (NumVectors, STS::zero ());
1470 for (
size_t k = 0; k < NumEntries; ++k) {
1471 const local_ordinal_type col = Indices[k];
1472 for (
size_t m = 0; m < NumVectors; ++m) {
1473 dtemp[m] += Values[k] * y2_ptr[col + y2_stride*m];
1477 for (
size_t m = 0; m < NumVectors; ++m) {
1478 y2_ptr[i + y2_stride*m] += DampingFactor_ * d_ptr[i] * (x_ptr[i + x_stride*m] - dtemp[m]);
1484 Tpetra::deep_copy (Y, *Y2);
1490 ArrayRCP<ArrayRCP<scalar_type> > y2_ptr = Y2->get2dViewNonConst ();
1491 ArrayRCP<ArrayRCP<const scalar_type> > x_ptr = X.get2dView ();
1493 for (
int j = 0; j < NumSweeps_; ++j) {
1496 if (Importer_.is_null ()) {
1499 Y2->doImport (Y, *Importer_, Tpetra::INSERT);
1503 if (! DoBackwardGS_) {
1504 for (
size_t ii = 0; ii < numActive; ++ii) {
1505 local_ordinal_type i = as<local_ordinal_type> (do_local ? rowInd[ii] : ii);
1507 A_->getLocalRowCopy (i, Indices (), Values (), NumEntries);
1509 for (
size_t m = 0; m < NumVectors; ++m) {
1510 scalar_type dtemp = STS::zero ();
1511 ArrayView<const scalar_type> x_local = (x_ptr())[m]();
1512 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
1514 for (
size_t k = 0; k < NumEntries; ++k) {
1515 const local_ordinal_type col = Indices[k];
1516 dtemp += Values[k] * y2_local[col];
1518 y2_local[i] += DampingFactor_ * d_ptr[i] * (x_local[i] - dtemp);
1525 for (ptrdiff_t ii = as<ptrdiff_t> (numActive) - 1; ii >= 0; --ii) {
1526 local_ordinal_type i = as<local_ordinal_type> (do_local ? rowInd[ii] : ii);
1529 A_->getLocalRowCopy (i, Indices (), Values (), NumEntries);
1531 for (
size_t m = 0; m < NumVectors; ++m) {
1532 scalar_type dtemp = STS::zero ();
1533 ArrayView<const scalar_type> x_local = (x_ptr())[m]();
1534 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
1536 for (
size_t k = 0; k < NumEntries; ++k) {
1537 const local_ordinal_type col = Indices[k];
1538 dtemp += Values[k] * y2_local[col];
1540 y2_local[i] += DampingFactor_ * d_ptr[i] * (x_local[i] - dtemp);
1547 Tpetra::deep_copy (Y, *Y2);
1554 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1555 const double numVectors = as<double> (X.getNumVectors ());
1556 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1557 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1558 ApplyFlops_ += 2.0 * NumSweeps_ * numVectors *
1559 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1563 template<
class MatrixType>
1565 Relaxation<MatrixType>::
1566 ApplyInverseGS_CrsMatrix (
const crs_matrix_type& A,
1567 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1568 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const 1571 const Tpetra::ESweepDirection direction =
1572 DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward;
1573 if (localSmoothingIndices_.is_null ()) {
1574 A.gaussSeidelCopy (Y, X, *Diagonal_, DampingFactor_, direction,
1575 NumSweeps_, ZeroStartingSolution_);
1578 A.reorderedGaussSeidelCopy (Y, X, *Diagonal_, localSmoothingIndices_ (),
1579 DampingFactor_, direction,
1580 NumSweeps_, ZeroStartingSolution_);
1594 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1595 const double numVectors = as<double> (X.getNumVectors ());
1596 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1597 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1598 ApplyFlops_ += NumSweeps_ * numVectors *
1599 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1602 template<
class MatrixType>
1604 Relaxation<MatrixType>::
1605 ApplyInverseGS_BlockCrsMatrix (
const block_crs_matrix_type& A,
1606 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1607 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
1611 using Teuchos::rcpFromRef;
1612 typedef Tpetra::Experimental::BlockMultiVector<scalar_type,
1613 local_ordinal_type, global_ordinal_type, node_type> BMV;
1614 typedef Tpetra::MultiVector<scalar_type,
1615 local_ordinal_type, global_ordinal_type, node_type> MV;
1624 BMV yBlock (Y, * (A.getGraph ()->getDomainMap ()), A.getBlockSize ());
1625 const BMV xBlock (X, * (A.getColMap ()), A.getBlockSize ());
1627 bool performImport =
false;
1629 if (Importer_.is_null ()) {
1630 yBlockCol = rcpFromRef (yBlock);
1633 if (yBlockColumnPointMap_.is_null () ||
1634 yBlockColumnPointMap_->getNumVectors () != yBlock.getNumVectors () ||
1635 yBlockColumnPointMap_->getBlockSize () != yBlock.getBlockSize ()) {
1636 yBlockColumnPointMap_ =
1637 rcp (
new BMV (* (A.getColMap ()), A.getBlockSize (),
1638 static_cast<local_ordinal_type
> (yBlock.getNumVectors ())));
1640 yBlockCol = yBlockColumnPointMap_;
1641 performImport =
true;
1644 if (ZeroStartingSolution_) {
1645 yBlockCol->putScalar (STS::zero ());
1647 else if (performImport) {
1648 yBlockCol->doImport (yBlock, *Importer_, Tpetra::INSERT);
1651 const Tpetra::ESweepDirection direction =
1652 DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward;
1654 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
1655 if (performImport && sweep > 0) {
1656 yBlockCol->doImport(yBlock, *Importer_, Tpetra::INSERT);
1658 A.localGaussSeidel (xBlock, *yBlockCol, blockDiag_,
1659 DampingFactor_, direction);
1660 if (performImport) {
1661 RCP<const MV> yBlockColPointDomain =
1662 yBlockCol->getMultiVectorView ().offsetView (A.getDomainMap (), 0);
1663 Tpetra::deep_copy (Y, *yBlockColPointDomain);
1671 template<
class MatrixType>
1672 void Relaxation<MatrixType>::MTGaussSeidel (
1673 const crs_matrix_type* crsMat,
1674 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1675 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& B,
1676 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& D,
1677 const scalar_type& dampingFactor,
1678 const Tpetra::ESweepDirection direction,
1679 const int numSweeps,
1680 const bool zeroInitialGuess)
const 1682 #ifdef HAVE_IFPACK2_EXPERIMENTAL_KOKKOSKERNELS_FEATURES 1683 using Teuchos::null;
1686 using Teuchos::rcpFromRef;
1687 using Teuchos::rcp_const_cast;
1689 typedef scalar_type Scalar;
1690 typedef local_ordinal_type LocalOrdinal;
1691 typedef global_ordinal_type GlobalOrdinal;
1692 typedef node_type Node;
1695 const char prefix[] =
"Ifpack2::Relaxation::(reordered)MTGaussSeidel: ";
1699 ! crsMat->isFillComplete (), std::runtime_error,
1700 prefix <<
"The matrix is not fill complete.");
1702 numSweeps < 0, std::invalid_argument,
1703 prefix <<
"The number of sweeps must be nonnegative, " 1704 "but you provided numSweeps = " << numSweeps <<
" < 0.");
1708 if (numSweeps == 0) {
1713 typedef typename Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> MV;
1714 typedef typename crs_matrix_type::import_type import_type;
1715 typedef typename crs_matrix_type::export_type export_type;
1716 typedef typename crs_matrix_type::map_type map_type;
1718 RCP<const import_type> importer = crsMat->getGraph ()->getImporter ();
1719 RCP<const export_type> exporter = crsMat->getGraph ()->getExporter ();
1721 ! exporter.is_null (), std::runtime_error,
1722 "This method's implementation currently requires that the matrix's row, " 1723 "domain, and range Maps be the same. This cannot be the case, because " 1724 "the matrix has a nontrivial Export object.");
1726 RCP<const map_type> domainMap = crsMat->getDomainMap ();
1727 RCP<const map_type> rangeMap = crsMat->getRangeMap ();
1728 RCP<const map_type> rowMap = crsMat->getGraph ()->getRowMap ();
1729 RCP<const map_type> colMap = crsMat->getGraph ()->getColMap ();
1732 #ifdef HAVE_IFPACK2_DEBUG 1738 ! X.getMap ()->isSameAs (*domainMap), std::runtime_error,
1739 "Ifpack2::Relaxation::MTGaussSeidel requires that the input " 1740 "multivector X be in the domain Map of the matrix.");
1742 ! B.getMap ()->isSameAs (*rangeMap), std::runtime_error,
1743 "Ifpack2::Relaxation::MTGaussSeidel requires that the input " 1744 "B be in the range Map of the matrix.");
1746 ! D.getMap ()->isSameAs (*rowMap), std::runtime_error,
1747 "Ifpack2::Relaxation::MTGaussSeidel requires that the input " 1748 "D be in the row Map of the matrix.");
1750 ! rowMap->isSameAs (*rangeMap), std::runtime_error,
1751 "Ifpack2::Relaxation::MTGaussSeidel requires that the row Map and the " 1752 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
1754 ! domainMap->isSameAs (*rangeMap), std::runtime_error,
1755 "Ifpack2::Relaxation::MTGaussSeidel requires that the domain Map and " 1756 "the range Map of the matrix be the same.");
1762 #endif // HAVE_IFPACK2_DEBUG 1773 RCP<MV> X_domainMap;
1774 bool copyBackOutput =
false;
1775 if (importer.is_null ()) {
1776 if (X.isConstantStride ()) {
1777 X_colMap = rcpFromRef (X);
1778 X_domainMap = rcpFromRef (X);
1785 if (zeroInitialGuess) {
1786 X_colMap->putScalar (ZERO);
1796 X_colMap =
rcp (
new MV (colMap, X.getNumVectors ()));
1800 X_domainMap = X_colMap;
1801 if (! zeroInitialGuess) {
1803 deep_copy (*X_domainMap , X);
1804 }
catch (std::exception& e) {
1805 std::ostringstream os;
1806 os <<
"Ifpack2::Relaxation::MTGaussSeidel: " 1807 "deep_copy(*X_domainMap, X) threw an exception: " 1808 << e.what () <<
".";
1812 copyBackOutput =
true;
1827 X_colMap =
rcp (
new MV (colMap, X.getNumVectors ()));
1829 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
1831 #ifdef HAVE_IFPACK2_DEBUG 1832 auto X_colMap_host_view = X_colMap->template getLocalView<Kokkos::HostSpace> ();
1833 auto X_domainMap_host_view = X_domainMap->template getLocalView<Kokkos::HostSpace> ();
1835 if (X_colMap->getLocalLength () != 0 && X_domainMap->getLocalLength ()) {
1837 X_colMap_host_view.ptr_on_device () != X_domainMap_host_view.ptr_on_device (),
1838 std::logic_error,
"Ifpack2::Relaxation::MTGaussSeidel: " 1839 "Pointer to start of column Map view of X is not equal to pointer to " 1840 "start of (domain Map view of) X. This may mean that " 1841 "Tpetra::MultiVector::offsetViewNonConst is broken. " 1842 "Please report this bug to the Tpetra developers.");
1846 X_colMap_host_view.dimension_0 () < X_domainMap_host_view.dimension_0 () ||
1847 X_colMap->getLocalLength () < X_domainMap->getLocalLength (),
1848 std::logic_error,
"Ifpack2::Relaxation::MTGaussSeidel: " 1849 "X_colMap has fewer local rows than X_domainMap. " 1850 "X_colMap_host_view.dimension_0() = " << X_colMap_host_view.dimension_0 ()
1851 <<
", X_domainMap_host_view.dimension_0() = " 1852 << X_domainMap_host_view.dimension_0 ()
1853 <<
", X_colMap->getLocalLength() = " << X_colMap->getLocalLength ()
1854 <<
", and X_domainMap->getLocalLength() = " 1855 << X_domainMap->getLocalLength ()
1856 <<
". This means that Tpetra::MultiVector::offsetViewNonConst " 1857 "is broken. Please report this bug to the Tpetra developers.");
1860 X_colMap->getNumVectors () != X_domainMap->getNumVectors (),
1861 std::logic_error,
"Ifpack2::Relaxation::MTGaussSeidel: " 1862 "X_colMap has a different number of columns than X_domainMap. " 1863 "X_colMap->getNumVectors() = " << X_colMap->getNumVectors ()
1864 <<
" != X_domainMap->getNumVectors() = " 1865 << X_domainMap->getNumVectors ()
1866 <<
". This means that Tpetra::MultiVector::offsetViewNonConst " 1867 "is broken. Please report this bug to the Tpetra developers.");
1868 #endif // HAVE_IFPACK2_DEBUG 1870 if (zeroInitialGuess) {
1872 X_colMap->putScalar (ZERO);
1882 X_colMap->doImport (X, *importer, Tpetra::CombineMode::INSERT);
1884 copyBackOutput =
true;
1891 if (B.isConstantStride ()) {
1892 B_in = rcpFromRef (B);
1899 RCP<MV> B_in_nonconst =
rcp (
new MV (rowMap, B.getNumVectors()));
1901 deep_copy (*B_in_nonconst, B);
1902 }
catch (std::exception& e) {
1903 std::ostringstream os;
1904 os <<
"Ifpack2::Relaxation::MTGaussSeidel: " 1905 "deep_copy(*B_in_nonconst, B) threw an exception: " 1906 << e.what () <<
".";
1909 B_in = rcp_const_cast<
const MV> (B_in_nonconst);
1923 kokkos_csr_matrix kcsr = crsMat->getLocalMatrix ();
1924 const size_t NumVectors = X.getNumVectors ();
1926 bool update_y_vector =
true;
1928 bool zero_x_vector =
false;
1929 for (
int sweep = 0; sweep < numSweeps; ++sweep) {
1930 if (! importer.is_null () && sweep > 0) {
1933 X_colMap->doImport (*X_domainMap, *importer, Tpetra::CombineMode::INSERT);
1987 for (
size_t indVec = 0; indVec < NumVectors; ++indVec){
1988 if (direction == Tpetra::Symmetric) {
1989 KokkosKernels::Experimental::Graph::symmetric_gauss_seidel_apply
1990 (kh.getRawPtr(), A_->getNodeNumRows(), A_->getNodeNumCols(),
1991 kcsr.graph.row_map, kcsr.graph.entries, kcsr.values,
1992 Kokkos::subview(X_colMap->template getLocalView<MyExecSpace> (), Kokkos::ALL (), indVec),
1993 Kokkos::subview(B_in->template getLocalView<MyExecSpace> (), Kokkos::ALL (), indVec),
1994 zero_x_vector, update_y_vector);
1996 else if (direction == Tpetra::Forward) {
1997 KokkosKernels::Experimental::Graph::forward_sweep_gauss_seidel_apply
1998 (kh.getRawPtr(), A_->getNodeNumRows(), A_->getNodeNumCols(),
1999 kcsr.graph.row_map,kcsr.graph.entries, kcsr.values,
2000 Kokkos::subview(X_colMap->template getLocalView<MyExecSpace> (), Kokkos::ALL (), indVec ),
2001 Kokkos::subview(B_in->template getLocalView<MyExecSpace> (), Kokkos::ALL (), indVec),
2002 zero_x_vector, update_y_vector);
2004 else if (direction == Tpetra::Backward) {
2005 KokkosKernels::Experimental::Graph::backward_sweep_gauss_seidel_apply
2006 (kh.getRawPtr(), A_->getNodeNumRows(), A_->getNodeNumCols(),
2007 kcsr.graph.row_map,kcsr.graph.entries, kcsr.values,
2008 Kokkos::subview(X_colMap->template getLocalView<MyExecSpace> (), Kokkos::ALL (), indVec ),
2009 Kokkos::subview(B_in->template getLocalView<MyExecSpace> (), Kokkos::ALL (), indVec),
2010 zero_x_vector, update_y_vector);
2014 true, std::invalid_argument,
2015 prefix <<
"The 'direction' enum does not have any of its valid " 2016 "values: Forward, Backward, or Symmetric.");
2020 if (NumVectors > 1){
2021 update_y_vector =
true;
2024 update_y_vector =
false;
2051 if (copyBackOutput) {
2053 deep_copy (X , *X_domainMap);
2054 }
catch (std::exception& e) {
2056 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) " 2057 "threw an exception: " << e.what ());
2063 template<
class MatrixType>
2064 void Relaxation<MatrixType>::ApplyInverseMTSGS_CrsMatrix (
2065 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& B,
2066 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X)
const {
2068 const crs_matrix_type* crsMat =
dynamic_cast<const crs_matrix_type*
> (&(*A_));
2070 crsMat == NULL, std::runtime_error,
"Ifpack2::Relaxation::compute: " 2071 "MT methods works for CRSMatrix Only.");
2074 const Tpetra::ESweepDirection direction = Tpetra::Symmetric;
2077 if (!localSmoothingIndices_.is_null ()) {
2078 std::cerr <<
"MT GaussSeidel ignores the given order" << std::endl;
2080 this->MTGaussSeidel (
2085 direction, NumSweeps_,
2086 ZeroStartingSolution_);
2088 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
2089 const double numVectors = as<double> (X.getNumVectors ());
2090 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
2091 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
2092 ApplyFlops_ += 2.0 * NumSweeps_ * numVectors *
2093 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
2097 template<
class MatrixType>
2098 void Relaxation<MatrixType>::ApplyInverseMTGS_CrsMatrix (
2099 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& B,
2100 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X)
const {
2102 const crs_matrix_type* crsMat =
dynamic_cast<const crs_matrix_type*
> (&(*A_));
2104 crsMat == NULL, std::runtime_error,
"Ifpack2::Relaxation::compute: " 2105 "MT methods works for CRSMatrix Only.");
2108 const Tpetra::ESweepDirection direction =
2109 DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward;
2112 if (!localSmoothingIndices_.is_null ()) {
2113 std::cerr <<
"MT GaussSeidel ignores the given order" << std::endl;
2115 this->MTGaussSeidel (
2120 direction, NumSweeps_,
2121 ZeroStartingSolution_);
2123 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
2124 const double numVectors = as<double> (X.getNumVectors ());
2125 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
2126 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
2127 ApplyFlops_ += NumSweeps_ * numVectors *
2128 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
2131 template<
class MatrixType>
2133 Relaxation<MatrixType>::
2134 ApplyInverseSGS (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
2135 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const 2137 typedef Relaxation<MatrixType> this_type;
2147 const block_crs_matrix_type* blockCrsMat =
dynamic_cast<const block_crs_matrix_type*
> (A_.getRawPtr());
2148 const crs_matrix_type* crsMat =
dynamic_cast<const crs_matrix_type*
> (&(*A_));
2149 if (blockCrsMat != NULL) {
2150 const_cast<this_type*
> (
this)->ApplyInverseSGS_BlockCrsMatrix(*blockCrsMat, X, Y);
2152 else if (crsMat != NULL) {
2153 ApplyInverseSGS_CrsMatrix (*crsMat, X, Y);
2155 ApplyInverseSGS_RowMatrix (X, Y);
2160 template<
class MatrixType>
2162 Relaxation<MatrixType>::
2163 ApplyInverseSGS_RowMatrix (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
2164 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const 2172 using Teuchos::rcpFromRef;
2173 typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
2178 if (ZeroStartingSolution_) {
2179 Y.putScalar (STS::zero ());
2182 const size_t NumVectors = X.getNumVectors ();
2183 const size_t maxLength = A_->getNodeMaxNumRowEntries ();
2184 Array<local_ordinal_type> Indices (maxLength);
2185 Array<scalar_type> Values (maxLength);
2188 const size_t numMyRows = A_->getNodeNumRows();
2189 const local_ordinal_type * rowInd = 0;
2190 size_t numActive = numMyRows;
2191 bool do_local = localSmoothingIndices_.is_null();
2193 rowInd = localSmoothingIndices_.getRawPtr();
2194 numActive = localSmoothingIndices_.size();
2200 if (Importer_.is_null ()) {
2202 Y2 =
rcp (
new MV (Y.getMap (), NumVectors,
false));
2207 Y2 =
rcp (
new MV (Importer_->getTargetMap (), NumVectors));
2211 Y2 = rcpFromRef (Y);
2215 ArrayRCP<const scalar_type> d_rcp = Diagonal_->get1dView();
2216 ArrayView<const scalar_type> d_ptr = d_rcp();
2219 bool constant_stride = X.isConstantStride() && Y2->isConstantStride();
2221 if(constant_stride) {
2223 size_t x_stride = X.getStride();
2224 size_t y2_stride = Y2->getStride();
2225 ArrayRCP<scalar_type> y2_rcp = Y2->get1dViewNonConst();
2226 ArrayRCP<const scalar_type> x_rcp = X.get1dView();
2227 ArrayView<scalar_type> y2_ptr = y2_rcp();
2228 ArrayView<const scalar_type> x_ptr = x_rcp();
2229 Array<scalar_type> dtemp(NumVectors,STS::zero());
2230 for (
int iter = 0; iter < NumSweeps_; ++iter) {
2233 if (Importer_.is_null ()) {
2235 Tpetra::deep_copy (*Y2, Y);
2237 Y2->doImport (Y, *Importer_, Tpetra::INSERT);
2241 for (
int j = 0; j < NumSweeps_; j++) {
2244 if (Importer_.is_null ()) {
2246 Tpetra::deep_copy (*Y2, Y);
2248 Y2->doImport (Y, *Importer_, Tpetra::INSERT);
2252 for (
size_t ii = 0; ii < numActive; ++ii) {
2253 local_ordinal_type i = as<local_ordinal_type>(do_local ? rowInd[ii] : ii);
2255 A_->getLocalRowCopy (i, Indices (), Values (), NumEntries);
2256 dtemp.assign(NumVectors,STS::zero());
2258 for (
size_t k = 0; k < NumEntries; ++k) {
2259 const local_ordinal_type col = Indices[k];
2260 for (
size_t m = 0; m < NumVectors; ++m) {
2261 dtemp[m] += Values[k] * y2_ptr[col + y2_stride*m];
2265 for (
size_t m = 0; m < NumVectors; ++m) {
2266 y2_ptr[i + y2_stride*m] += DampingFactor_ * d_ptr[i] * (x_ptr[i + x_stride*m] - dtemp[m]);
2272 for (ptrdiff_t ii = as<ptrdiff_t> (numActive) - 1; ii >= 0; --ii) {
2273 local_ordinal_type i = as<local_ordinal_type>(do_local ? rowInd[ii] : ii);
2275 A_->getLocalRowCopy (i, Indices (), Values (), NumEntries);
2276 dtemp.assign(NumVectors,STS::zero());
2278 for (
size_t k = 0; k < NumEntries; ++k) {
2279 const local_ordinal_type col = Indices[k];
2280 for (
size_t m = 0; m < NumVectors; ++m) {
2281 dtemp[m] += Values[k] * y2_ptr[col + y2_stride*m];
2285 for (
size_t m = 0; m < NumVectors; ++m) {
2286 y2_ptr[i + y2_stride*m] += DampingFactor_ * d_ptr[i] * (x_ptr[i + x_stride*m] - dtemp[m]);
2293 Tpetra::deep_copy (Y, *Y2);
2299 ArrayRCP<ArrayRCP<scalar_type> > y2_ptr = Y2->get2dViewNonConst ();
2300 ArrayRCP<ArrayRCP<const scalar_type> > x_ptr = X.get2dView ();
2302 for (
int iter = 0; iter < NumSweeps_; ++iter) {
2305 if (Importer_.is_null ()) {
2307 Tpetra::deep_copy (*Y2, Y);
2309 Y2->doImport (Y, *Importer_, Tpetra::INSERT);
2313 for (
size_t ii = 0; ii < numActive; ++ii) {
2314 local_ordinal_type i = as<local_ordinal_type>(do_local ? rowInd[ii] : ii);
2315 const scalar_type diag = d_ptr[i];
2317 A_->getLocalRowCopy (as<local_ordinal_type> (i), Indices (), Values (), NumEntries);
2319 for (
size_t m = 0; m < NumVectors; ++m) {
2320 scalar_type dtemp = STS::zero ();
2321 ArrayView<const scalar_type> x_local = (x_ptr())[m]();
2322 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
2324 for (
size_t k = 0; k < NumEntries; ++k) {
2325 const local_ordinal_type col = Indices[k];
2326 dtemp += Values[k] * y2_local[col];
2328 y2_local[i] += DampingFactor_ * (x_local[i] - dtemp) * diag;
2334 for (ptrdiff_t ii = as<ptrdiff_t> (numActive) - 1; ii >= 0; --ii) {
2335 local_ordinal_type i = as<local_ordinal_type>(do_local ? rowInd[ii] : ii);
2336 const scalar_type diag = d_ptr[i];
2338 A_->getLocalRowCopy (as<local_ordinal_type> (i), Indices (), Values (), NumEntries);
2340 for (
size_t m = 0; m < NumVectors; ++m) {
2341 scalar_type dtemp = STS::zero ();
2342 ArrayView<const scalar_type> x_local = (x_ptr())[m]();
2343 ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
2345 for (
size_t k = 0; k < NumEntries; ++k) {
2346 const local_ordinal_type col = Indices[k];
2347 dtemp += Values[k] * y2_local[col];
2349 y2_local[i] += DampingFactor_ * (x_local[i] - dtemp) * diag;
2355 Tpetra::deep_copy (Y, *Y2);
2361 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
2362 const double numVectors = as<double> (X.getNumVectors ());
2363 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
2364 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
2365 ApplyFlops_ += 2.0 * NumSweeps_ * numVectors *
2366 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
2370 template<
class MatrixType>
2372 Relaxation<MatrixType>::
2373 ApplyInverseSGS_CrsMatrix (
const crs_matrix_type& A,
2374 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
2375 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const 2378 const Tpetra::ESweepDirection direction = Tpetra::Symmetric;
2379 if (localSmoothingIndices_.is_null ()) {
2380 A.gaussSeidelCopy (Y, X, *Diagonal_, DampingFactor_, direction,
2381 NumSweeps_, ZeroStartingSolution_);
2384 A.reorderedGaussSeidelCopy (Y, X, *Diagonal_, localSmoothingIndices_ (),
2385 DampingFactor_, direction,
2386 NumSweeps_, ZeroStartingSolution_);
2403 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
2404 const double numVectors = as<double> (X.getNumVectors ());
2405 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
2406 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
2407 ApplyFlops_ += 2.0 * NumSweeps_ * numVectors *
2408 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
2411 template<
class MatrixType>
2413 Relaxation<MatrixType>::
2414 ApplyInverseSGS_BlockCrsMatrix (
const block_crs_matrix_type& A,
2415 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
2416 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
2420 using Teuchos::rcpFromRef;
2421 typedef Tpetra::Experimental::BlockMultiVector<scalar_type,
2422 local_ordinal_type, global_ordinal_type, node_type> BMV;
2423 typedef Tpetra::MultiVector<scalar_type,
2424 local_ordinal_type, global_ordinal_type, node_type> MV;
2433 BMV yBlock (Y, * (A.getGraph ()->getDomainMap ()), A.getBlockSize ());
2434 const BMV xBlock (X, * (A.getColMap ()), A.getBlockSize ());
2436 bool performImport =
false;
2438 if (Importer_.is_null ()) {
2439 yBlockCol = Teuchos::rcpFromRef (yBlock);
2442 if (yBlockColumnPointMap_.is_null () ||
2443 yBlockColumnPointMap_->getNumVectors () != yBlock.getNumVectors () ||
2444 yBlockColumnPointMap_->getBlockSize () != yBlock.getBlockSize ()) {
2445 yBlockColumnPointMap_ =
2446 rcp (
new BMV (* (A.getColMap ()), A.getBlockSize (),
2447 static_cast<local_ordinal_type
> (yBlock.getNumVectors ())));
2449 yBlockCol = yBlockColumnPointMap_;
2450 performImport =
true;
2453 if (ZeroStartingSolution_) {
2454 yBlockCol->putScalar (STS::zero ());
2456 else if (performImport) {
2457 yBlockCol->doImport (yBlock, *Importer_, Tpetra::INSERT);
2461 const Tpetra::ESweepDirection direction = Tpetra::Symmetric;
2463 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
2464 if (performImport && sweep > 0) {
2465 yBlockCol->doImport (yBlock, *Importer_, Tpetra::INSERT);
2467 A.localGaussSeidel (xBlock, *yBlockCol, blockDiag_,
2468 DampingFactor_, direction);
2469 if (performImport) {
2470 RCP<const MV> yBlockColPointDomain =
2471 yBlockCol->getMultiVectorView ().offsetView (A.getDomainMap (), 0);
2472 MV yBlockView = yBlock.getMultiVectorView ();
2473 Tpetra::deep_copy (yBlockView, *yBlockColPointDomain);
2479 template<
class MatrixType>
2482 std::ostringstream os;
2487 os <<
"\"Ifpack2::Relaxation\": {";
2489 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", " 2490 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
2496 if (PrecType_ == Ifpack2::Details::JACOBI) {
2498 }
else if (PrecType_ == Ifpack2::Details::GS) {
2499 os <<
"Gauss-Seidel";
2500 }
else if (PrecType_ == Ifpack2::Details::SGS) {
2501 os <<
"Symmetric Gauss-Seidel";
2502 }
else if (PrecType_ == Ifpack2::Details::MTGS) {
2503 os <<
"MT Gauss-Seidel";
2504 }
else if (PrecType_ == Ifpack2::Details::MTSGS) {
2505 os <<
"MT Symmetric Gauss-Seidel";
2511 os <<
", " <<
"sweeps: " << NumSweeps_ <<
", " 2512 <<
"damping factor: " << DampingFactor_ <<
", ";
2514 os <<
"use l1: " << DoL1Method_ <<
", " 2515 <<
"l1 eta: " << L1Eta_ <<
", ";
2518 if (A_.is_null ()) {
2519 os <<
"Matrix: null";
2522 os <<
"Global matrix dimensions: [" 2523 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]" 2524 <<
", Global nnz: " << A_->getGlobalNumEntries();
2532 template<
class MatrixType>
2551 const int myRank = this->getComm ()->getRank ();
2563 out <<
"\"Ifpack2::Relaxation\":" << endl;
2565 out <<
"MatrixType: \"" << TypeNameTraits<MatrixType>::name () <<
"\"" 2567 if (this->getObjectLabel () !=
"") {
2568 out <<
"Label: " << this->getObjectLabel () << endl;
2570 out <<
"Initialized: " << (isInitialized () ?
"true" :
"false") << endl
2571 <<
"Computed: " << (isComputed () ?
"true" :
"false") << endl
2572 <<
"Parameters: " << endl;
2575 out <<
"\"relaxation: type\": ";
2576 if (PrecType_ == Ifpack2::Details::JACOBI) {
2578 }
else if (PrecType_ == Ifpack2::Details::GS) {
2579 out <<
"Gauss-Seidel";
2580 }
else if (PrecType_ == Ifpack2::Details::SGS) {
2581 out <<
"Symmetric Gauss-Seidel";
2582 }
else if (PrecType_ == Ifpack2::Details::MTGS) {
2583 out <<
"MT Gauss-Seidel";
2584 }
else if (PrecType_ == Ifpack2::Details::MTSGS) {
2585 out <<
"MT Symmetric Gauss-Seidel";
2592 <<
"\"relaxation: sweeps\": " << NumSweeps_ << endl
2593 <<
"\"relaxation: damping factor\": " << DampingFactor_ << endl
2594 <<
"\"relaxation: min diagonal value\": " << MinDiagonalValue_ << endl
2595 <<
"\"relaxation: zero starting solution\": " << ZeroStartingSolution_ << endl
2596 <<
"\"relaxation: backward mode\": " << DoBackwardGS_ << endl
2597 <<
"\"relaxation: use l1\": " << DoL1Method_ << endl
2598 <<
"\"relaxation: l1 eta\": " << L1Eta_ << endl;
2600 out <<
"Computed quantities:" << endl;
2603 out <<
"Global number of rows: " << A_->getGlobalNumRows () << endl
2604 <<
"Global number of columns: " << A_->getGlobalNumCols () << endl;
2606 if (checkDiagEntries_ && isComputed ()) {
2607 out <<
"Properties of input diagonal entries:" << endl;
2610 out <<
"Magnitude of minimum-magnitude diagonal entry: " 2611 << globalMinMagDiagEntryMag_ << endl
2612 <<
"Magnitude of maximum-magnitude diagonal entry: " 2613 << globalMaxMagDiagEntryMag_ << endl
2614 <<
"Number of diagonal entries with small magnitude: " 2615 << globalNumSmallDiagEntries_ << endl
2616 <<
"Number of zero diagonal entries: " 2617 << globalNumZeroDiagEntries_ << endl
2618 <<
"Number of diagonal entries with negative real part: " 2619 << globalNumNegDiagEntries_ << endl
2620 <<
"Abs 2-norm diff between computed and actual inverse " 2621 <<
"diagonal: " << globalDiagNormDiff_ << endl;
2624 if (isComputed ()) {
2625 out <<
"Saved diagonal offsets: " 2626 << (savedDiagOffsets_ ?
"true" :
"false") << endl;
2628 out <<
"Call counts and total times (in seconds): " << endl;
2631 out <<
"initialize: " << endl;
2634 out <<
"Call count: " << NumInitialize_ << endl;
2635 out <<
"Total time: " << InitializeTime_ << endl;
2637 out <<
"compute: " << endl;
2640 out <<
"Call count: " << NumCompute_ << endl;
2641 out <<
"Total time: " << ComputeTime_ << endl;
2643 out <<
"apply: " << endl;
2646 out <<
"Call count: " << NumApply_ << endl;
2647 out <<
"Total time: " << ApplyTime_ << endl;
2655 #define IFPACK2_RELAXATION_INSTANT(S,LO,GO,N) \ 2656 template class Ifpack2::Relaxation< Tpetra::RowMatrix<S, LO, GO, N> >; 2658 #endif // IFPACK2_RELAXATION_DEF_HPP
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_Relaxation_def.hpp:2480
bool is_null(const boost::shared_ptr< T > &p)
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_Relaxation_def.hpp:377
int getNumInitialize() const
Total number of calls to initialize().
Definition: Ifpack2_Relaxation_def.hpp:393
basic_OSTab< char > OSTab
static magnitudeType eps()
double getApplyTime() const
Total time in seconds spent in all calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:423
Relaxation(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_Relaxation_def.hpp:175
void compute()
Compute the preconditioner.
Definition: Ifpack2_Relaxation_def.hpp:827
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
Apply the preconditioner to X, returning the result in Y.
Definition: Ifpack2_Relaxation_def.hpp:443
virtual void printDoc(std::string const &docString, std::ostream &out) const=0
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool hasTransposeApply() const
Whether apply() and applyMat() let you apply the transpose or conjugate transpose.
Definition: Ifpack2_Relaxation_def.hpp:387
virtual const std::string getXMLTypeName() const=0
int getNumApply() const
Total number of calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:405
static std::ostream & printLines(std::ostream &os, const std::string &linePrefix, const std::string &lines)
std::string * getRawPtr()
double getInitializeTime() const
Total time in seconds spent in all calls to initialize().
Definition: Ifpack2_Relaxation_def.hpp:411
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Relaxation_def.hpp:155
Ifpack2 implementation details.
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_Relaxation_def.hpp:364
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix to be preconditioned.
Definition: Ifpack2_Relaxation_def.hpp:355
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:250
void setParameters(const Teuchos::ParameterList ¶ms)
Set the relaxation / preconditioner parameters.
Definition: Ifpack2_Relaxation_def.hpp:334
LinearOp zero(const VectorSpace &vs)
File for utility functions.
double getApplyFlops() const
Total number of floating-point operations over all calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:435
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a list of all the parameters that this class accepts.
Definition: Ifpack2_Relaxation_def.hpp:220
virtual ValidStringsList validStringValues() const=0
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:247
const std::type_info & type() const
any & getAny(bool activeQry=true)
double getComputeFlops() const
Total number of floating-point operations over all calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:429
double getComputeTime() const
Total time in seconds spent in all calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:417
int getNumCompute() const
Total number of calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:399
std::string typeName() const
TypeTo as(const TypeFrom &t)
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:244
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix and vectors are distributed.
Definition: Ifpack2_Relaxation_def.hpp:344
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:241
virtual ~Relaxation()
Destructor.
Definition: Ifpack2_Relaxation_def.hpp:215
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
Apply the preconditioner to X, returning the result in Y.
Definition: Ifpack2_Relaxation_def.hpp:541
Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices.
Definition: Ifpack2_Relaxation_decl.hpp:226
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object's attributes to the given output stream.
Definition: Ifpack2_Relaxation_def.hpp:2535
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
void getValidParameters(Teuchos::ParameterList ¶ms)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:50
void initialize()
Initialize the preconditioner.
Definition: Ifpack2_Relaxation_def.hpp:562
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_Relaxation_decl.hpp:253
virtual void validate(ParameterEntry const &entry, std::string const ¶mName, std::string const &sublistName) const=0