41 #ifndef IFPACK2_RELAXATION_DEF_HPP
42 #define IFPACK2_RELAXATION_DEF_HPP
44 #include "Teuchos_StandardParameterEntryValidators.hpp"
46 #include "Tpetra_CrsMatrix.hpp"
47 #include "Tpetra_BlockCrsMatrix.hpp"
48 #include "Tpetra_BlockView.hpp"
50 #include "Ifpack2_Details_getCrsMatrix.hpp"
51 #include "MatrixMarket_Tpetra.hpp"
52 #include "Tpetra_Details_residual.hpp"
56 #include "KokkosSparse_gauss_seidel.hpp"
63 NonnegativeIntValidator () {}
73 const std::string& paramName,
74 const std::string& sublistName)
const
81 anyVal.
type () !=
typeid (int),
83 "Parameter \"" << paramName <<
"\" in sublist \"" << sublistName
84 <<
"\" has the wrong type." << endl <<
"Parameter: " << paramName
85 << endl <<
"Type specified: " << entryName << endl
86 <<
"Type required: int" << endl);
88 const int val = Teuchos::any_cast<
int> (anyVal);
91 "Parameter \"" << paramName <<
"\" in sublist \"" << sublistName
92 <<
"\" is negative." << endl <<
"Parameter: " << paramName
93 << endl <<
"Value specified: " << val << endl
94 <<
"Required range: [0, INT_MAX]" << endl);
99 return "NonnegativeIntValidator";
104 printDoc (
const std::string& docString,
105 std::ostream &out)
const
108 out <<
"#\tValidator Used: " << std::endl;
109 out <<
"#\t\tNonnegativeIntValidator" << std::endl;
116 template<class Scalar, const bool isOrdinal=Teuchos::ScalarTraits<Scalar>::isOrdinal>
120 static const Scalar eps ();
124 template<
class Scalar>
125 class SmallTraits<Scalar, true> {
127 static const Scalar eps () {
133 template<
class Scalar>
134 class SmallTraits<Scalar, false> {
136 static const Scalar eps () {
142 template<
class Scalar,
144 struct RealTraits {};
146 template<
class Scalar>
147 struct RealTraits<Scalar, false> {
148 using val_type = Scalar;
149 using mag_type = Scalar;
150 static KOKKOS_INLINE_FUNCTION mag_type real (
const val_type& z) {
155 template<
class Scalar>
156 struct RealTraits<Scalar, true> {
157 using val_type = Scalar;
159 static KOKKOS_INLINE_FUNCTION mag_type real (
const val_type& z) {
160 return Kokkos::ArithTraits<val_type>::real (z);
164 template<
class Scalar>
165 KOKKOS_INLINE_FUNCTION
typename RealTraits<Scalar>::mag_type
166 getRealValue (
const Scalar& z) {
167 return RealTraits<Scalar>::real (z);
174 template<
class MatrixType>
176 Relaxation<MatrixType>::
177 updateCachedMultiVector (
const Teuchos::RCP<
const Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>>& map,
178 size_t numVecs)
const
182 if (cachedMV_.is_null () ||
183 map.get () != cachedMV_->getMap ().get () ||
184 cachedMV_->getNumVectors () != numVecs) {
185 using MV = Tpetra::MultiVector<scalar_type, local_ordinal_type,
186 global_ordinal_type, node_type>;
187 cachedMV_ =
Teuchos::rcp (
new MV (map, numVecs,
false));
191 template<
class MatrixType>
196 Importer_ = Teuchos::null;
197 pointImporter_ = Teuchos::null;
198 Diagonal_ = Teuchos::null;
199 isInitialized_ =
false;
201 diagOffsets_ = Kokkos::View<size_t*, device_type>();
202 savedDiagOffsets_ =
false;
203 hasBlockCrsMatrix_ =
false;
204 serialGaussSeidel_ = Teuchos::null;
206 IsParallel_ = (A->getRowMap ()->getComm ()->getSize () > 1);
212 template<
class MatrixType>
216 IsParallel_ ((A.
is_null () || A->getRowMap ().
is_null () || A->getRowMap ()->getComm ().
is_null ()) ?
218 A->getRowMap ()->getComm ()->getSize () > 1)
220 this->setObjectLabel (
"Ifpack2::Relaxation");
224 template<
class MatrixType>
230 using Teuchos::parameterList;
233 using Teuchos::rcp_const_cast;
234 using Teuchos::rcp_implicit_cast;
235 using Teuchos::setStringToIntegralParameter;
238 if (validParams_.is_null ()) {
239 RCP<ParameterList> pl = parameterList (
"Ifpack2::Relaxation");
243 Array<std::string> precTypes (8);
244 precTypes[0] =
"Jacobi";
245 precTypes[1] =
"Gauss-Seidel";
246 precTypes[2] =
"Symmetric Gauss-Seidel";
247 precTypes[3] =
"MT Gauss-Seidel";
248 precTypes[4] =
"MT Symmetric Gauss-Seidel";
249 precTypes[5] =
"Richardson";
250 precTypes[6] =
"Two-stage Gauss-Seidel";
251 precTypes[7] =
"Two-stage Symmetric Gauss-Seidel";
252 Array<Details::RelaxationType> precTypeEnums (8);
253 precTypeEnums[0] = Details::JACOBI;
254 precTypeEnums[1] = Details::GS;
255 precTypeEnums[2] = Details::SGS;
256 precTypeEnums[3] = Details::MTGS;
257 precTypeEnums[4] = Details::MTSGS;
258 precTypeEnums[5] = Details::RICHARDSON;
259 precTypeEnums[6] = Details::GS2;
260 precTypeEnums[7] = Details::SGS2;
261 const std::string defaultPrecType (
"Jacobi");
262 setStringToIntegralParameter<Details::RelaxationType> (
"relaxation: type",
263 defaultPrecType,
"Relaxation method", precTypes (), precTypeEnums (),
266 const int numSweeps = 1;
267 RCP<PEV> numSweepsValidator =
268 rcp_implicit_cast<PEV> (
rcp (
new NonnegativeIntValidator));
269 pl->set (
"relaxation: sweeps", numSweeps,
"Number of relaxation sweeps",
270 rcp_const_cast<const PEV> (numSweepsValidator));
273 const int numOuterSweeps = 1;
274 RCP<PEV> numOuterSweepsValidator =
275 rcp_implicit_cast<PEV> (
rcp (
new NonnegativeIntValidator));
276 pl->set (
"relaxation: outer sweeps", numOuterSweeps,
"Number of outer local relaxation sweeps for two-stage GS",
277 rcp_const_cast<const PEV> (numOuterSweepsValidator));
279 const int numInnerSweeps = 1;
280 RCP<PEV> numInnerSweepsValidator =
281 rcp_implicit_cast<PEV> (
rcp (
new NonnegativeIntValidator));
282 pl->set (
"relaxation: inner sweeps", numInnerSweeps,
"Number of inner local relaxation sweeps for two-stage GS",
283 rcp_const_cast<const PEV> (numInnerSweepsValidator));
285 const scalar_type innerDampingFactor = STS::one ();
286 pl->set (
"relaxation: inner damping factor", innerDampingFactor,
"Damping factor for the inner sweep of two-stage GS");
288 const bool innerSpTrsv =
false;
289 pl->set (
"relaxation: inner sparse-triangular solve", innerSpTrsv,
"Specify whether to use sptrsv instead of JR iterations for two-stage GS");
291 const bool compactForm =
false;
292 pl->set (
"relaxation: compact form", compactForm,
"Specify whether to use compact form of recurrence for two-stage GS");
295 pl->set (
"relaxation: damping factor", dampingFactor);
297 const bool zeroStartingSolution =
true;
298 pl->set (
"relaxation: zero starting solution", zeroStartingSolution);
300 const bool doBackwardGS =
false;
301 pl->set (
"relaxation: backward mode", doBackwardGS);
303 const bool doL1Method =
false;
304 pl->set (
"relaxation: use l1", doL1Method);
306 const magnitude_type l1eta = (STM::one() + STM::one() + STM::one()) /
307 (STM::one() + STM::one());
308 pl->set (
"relaxation: l1 eta", l1eta);
311 pl->set (
"relaxation: min diagonal value", minDiagonalValue);
313 const bool fixTinyDiagEntries =
false;
314 pl->set (
"relaxation: fix tiny diagonal entries", fixTinyDiagEntries);
316 const bool checkDiagEntries =
false;
317 pl->set (
"relaxation: check diagonal entries", checkDiagEntries);
320 pl->set(
"relaxation: local smoothing indices", localSmoothingIndices);
322 const bool is_matrix_structurally_symmetric =
false;
323 pl->set(
"relaxation: symmetric matrix structure", is_matrix_structurally_symmetric);
325 const bool ifpack2_dump_matrix =
false;
326 pl->set(
"relaxation: ifpack2 dump matrix", ifpack2_dump_matrix);
328 const int cluster_size = 1;
329 pl->set(
"relaxation: mtgs cluster size", cluster_size);
331 pl->set(
"relaxation: mtgs coloring algorithm",
"Default");
333 const int long_row_threshold = 0;
334 pl->set(
"relaxation: long row threshold", long_row_threshold);
336 const bool timer_for_apply =
true;
337 pl->set(
"timer for apply", timer_for_apply);
339 validParams_ = rcp_const_cast<
const ParameterList> (pl);
345 template<
class MatrixType>
348 using Teuchos::getIntegralValue;
349 using Teuchos::getStringValue;
352 typedef scalar_type ST;
354 if (pl.
isType<
double>(
"relaxation: damping factor")) {
357 ST df = pl.
get<
double>(
"relaxation: damping factor");
358 pl.
remove(
"relaxation: damping factor");
359 pl.
set(
"relaxation: damping factor",df);
364 const Details::RelaxationType precType =
365 getIntegralValue<Details::RelaxationType> (pl,
"relaxation: type");
366 const std::string precTypeStr = getStringValue<Details::RelaxationType>(pl,
"relaxation: type");
368 pl.set<std::string>(
"relaxation: type", precTypeStr);
369 pl.get<std::string>(
"relaxation: type");
370 const int numSweeps = pl.get<
int> (
"relaxation: sweeps");
371 const ST dampingFactor = pl.get<ST> (
"relaxation: damping factor");
372 const bool zeroStartSol = pl.get<
bool> (
"relaxation: zero starting solution");
373 const bool doBackwardGS = pl.get<
bool> (
"relaxation: backward mode");
374 const bool doL1Method = pl.get<
bool> (
"relaxation: use l1");
375 const magnitude_type l1Eta = pl.get<magnitude_type> (
"relaxation: l1 eta");
376 const ST minDiagonalValue = pl.get<ST> (
"relaxation: min diagonal value");
377 const bool fixTinyDiagEntries = pl.get<
bool> (
"relaxation: fix tiny diagonal entries");
378 const bool checkDiagEntries = pl.get<
bool> (
"relaxation: check diagonal entries");
379 const bool is_matrix_structurally_symmetric = pl.get<
bool> (
"relaxation: symmetric matrix structure");
380 const bool ifpack2_dump_matrix = pl.get<
bool> (
"relaxation: ifpack2 dump matrix");
381 const bool timer_for_apply = pl.get<
bool> (
"timer for apply");
382 int cluster_size = 1;
383 if(pl.isParameter (
"relaxation: mtgs cluster size"))
384 cluster_size = pl.get<
int> (
"relaxation: mtgs cluster size");
385 int long_row_threshold = 0;
386 if(pl.isParameter (
"relaxation: long row threshold"))
387 long_row_threshold = pl.get<
int> (
"relaxation: long row threshold");
388 std::string color_algo_name = pl.get<std::string>(
"relaxation: mtgs coloring algorithm");
390 for(
char& c : color_algo_name)
392 if(color_algo_name ==
"default")
393 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_DEFAULT;
394 else if(color_algo_name ==
"serial")
395 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_SERIAL;
396 else if(color_algo_name ==
"vb")
397 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VB;
398 else if(color_algo_name ==
"vbbit")
399 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBBIT;
400 else if(color_algo_name ==
"vbcs")
401 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBCS;
402 else if(color_algo_name ==
"vbd")
403 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBD;
404 else if(color_algo_name ==
"vbdbit")
405 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBDBIT;
406 else if(color_algo_name ==
"eb")
407 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_EB;
410 std::ostringstream msg;
411 msg <<
"Ifpack2::Relaxation: 'relaxation: mtgs coloring algorithm' = '" << color_algo_name <<
"' is not valid.\n";
412 msg <<
"Choices (not case sensitive) are: Default, Serial, VB, VBBIT, VBCS, VBD, VBDBIT, EB.";
414 true, std::invalid_argument, msg.str());
420 if (!std::is_same<double, ST>::value && pl.isType<
double>(
"relaxation: inner damping factor")) {
423 ST df = pl.
get<
double>(
"relaxation: inner damping factor");
424 pl.remove(
"relaxation: inner damping factor");
425 pl.set(
"relaxation: inner damping factor",df);
428 if (long_row_threshold > 0) {
430 cluster_size != 1, std::invalid_argument,
"Ifpack2::Relaxation: "
431 "Requested long row MTGS/MTSGS algorithm and cluster GS/SGS, but those are not compatible.");
433 precType != Details::RelaxationType::MTGS && precType != Details::RelaxationType::MTSGS,
434 std::invalid_argument,
"Ifpack2::Relaxation: "
435 "Requested long row MTGS/MTSGS algorithm, but this is only compatible with preconditioner types "
436 "'MT Gauss-Seidel' and 'MT Symmetric Gauss-Seidel'.");
439 const ST innerDampingFactor = pl.get<ST> (
"relaxation: inner damping factor");
440 const int numInnerSweeps = pl.get<
int> (
"relaxation: inner sweeps");
441 const int numOuterSweeps = pl.get<
int> (
"relaxation: outer sweeps");
442 const bool innerSpTrsv = pl.get<
bool> (
"relaxation: inner sparse-triangular solve");
443 const bool compactForm = pl.get<
bool> (
"relaxation: compact form");
446 PrecType_ = precType;
447 NumSweeps_ = numSweeps;
448 DampingFactor_ = dampingFactor;
449 ZeroStartingSolution_ = zeroStartSol;
450 DoBackwardGS_ = doBackwardGS;
451 DoL1Method_ = doL1Method;
453 MinDiagonalValue_ = minDiagonalValue;
454 fixTinyDiagEntries_ = fixTinyDiagEntries;
455 checkDiagEntries_ = checkDiagEntries;
456 clusterSize_ = cluster_size;
457 longRowThreshold_ = long_row_threshold;
458 is_matrix_structurally_symmetric_ = is_matrix_structurally_symmetric;
459 ifpack2_dump_matrix_ = ifpack2_dump_matrix;
460 localSmoothingIndices_ = localSmoothingIndices;
462 NumInnerSweeps_ = numInnerSweeps;
463 NumOuterSweeps_ = numOuterSweeps;
464 InnerSpTrsv_ = innerSpTrsv;
465 InnerDampingFactor_ = innerDampingFactor;
466 CompactForm_ = compactForm;
467 TimerForApply_ = timer_for_apply;
471 template<
class MatrixType>
476 this->setParametersImpl(const_cast<Teuchos::ParameterList&>(pl));
480 template<
class MatrixType>
484 A_.
is_null (), std::runtime_error,
"Ifpack2::Relaxation::getComm: "
485 "The input matrix A is null. Please call setMatrix() with a nonnull "
486 "input matrix before calling this method.");
487 return A_->getRowMap ()->getComm ();
491 template<
class MatrixType>
498 template<
class MatrixType>
499 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
500 typename MatrixType::global_ordinal_type,
501 typename MatrixType::node_type> >
504 A_.
is_null (), std::runtime_error,
"Ifpack2::Relaxation::getDomainMap: "
505 "The input matrix A is null. Please call setMatrix() with a nonnull "
506 "input matrix before calling this method.");
507 return A_->getDomainMap ();
511 template<
class MatrixType>
512 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
513 typename MatrixType::global_ordinal_type,
514 typename MatrixType::node_type> >
517 A_.
is_null (), std::runtime_error,
"Ifpack2::Relaxation::getRangeMap: "
518 "The input matrix A is null. Please call setMatrix() with a nonnull "
519 "input matrix before calling this method.");
520 return A_->getRangeMap ();
524 template<
class MatrixType>
530 template<
class MatrixType>
532 return(NumInitialize_);
536 template<
class MatrixType>
542 template<
class MatrixType>
548 template<
class MatrixType>
550 return(InitializeTime_);
554 template<
class MatrixType>
556 return(ComputeTime_);
560 template<
class MatrixType>
566 template<
class MatrixType>
568 return(ComputeFlops_);
572 template<
class MatrixType>
579 template<
class MatrixType>
582 A_.
is_null (), std::runtime_error,
"Ifpack2::Relaxation::getNodeSmootherComplexity: "
583 "The input matrix A is null. Please call setMatrix() with a nonnull "
584 "input matrix, then call compute(), before calling this method.");
586 return A_->getLocalNumRows() + A_->getLocalNumEntries();
590 template<
class MatrixType>
593 apply (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
594 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
602 using Teuchos::rcpFromRef;
606 A_.
is_null (), std::runtime_error,
"Ifpack2::Relaxation::apply: "
607 "The input matrix A is null. Please call setMatrix() with a nonnull "
608 "input matrix, then call compute(), before calling this method.");
612 "Ifpack2::Relaxation::apply: You must call compute() on this Ifpack2 "
613 "preconditioner instance before you may call apply(). You may call "
614 "isComputed() to find out if compute() has been called already.");
615 TEUCHOS_TEST_FOR_EXCEPTION(
616 X.getNumVectors() != Y.getNumVectors(),
618 "Ifpack2::Relaxation::apply: X and Y have different numbers of columns. "
619 "X has " << X.getNumVectors() <<
" columns, but Y has "
620 << Y.getNumVectors() <<
" columns.");
621 TEUCHOS_TEST_FOR_EXCEPTION(
622 beta != STS::zero (), std::logic_error,
623 "Ifpack2::Relaxation::apply: beta = " << beta <<
" != 0 case not "
627 const std::string timerName (
"Ifpack2::Relaxation::apply");
628 if (TimerForApply_) {
644 if (alpha == STS::zero ()) {
646 Y.putScalar (STS::zero ());
656 Xcopy = rcpFromRef (X);
664 case Ifpack2::Details::JACOBI:
665 ApplyInverseJacobi(*Xcopy,Y);
667 case Ifpack2::Details::GS:
668 ApplyInverseSerialGS(*Xcopy, Y, DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward);
670 case Ifpack2::Details::SGS:
671 ApplyInverseSerialGS(*Xcopy, Y, Tpetra::Symmetric);
673 case Ifpack2::Details::MTGS:
674 case Ifpack2::Details::GS2:
675 ApplyInverseMTGS_CrsMatrix(*Xcopy, Y, DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward);
677 case Ifpack2::Details::MTSGS:
678 case Ifpack2::Details::SGS2:
679 ApplyInverseMTGS_CrsMatrix(*Xcopy, Y, Tpetra::Symmetric);
681 case Ifpack2::Details::RICHARDSON:
682 ApplyInverseRichardson(*Xcopy,Y);
686 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
687 "Ifpack2::Relaxation::apply: Invalid preconditioner type enum value "
688 << PrecType_ <<
". Please report this bug to the Ifpack2 developers.");
690 if (alpha != STS::one ()) {
692 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
693 const double numVectors = as<double> (Y.getNumVectors ());
694 ApplyFlops_ += numGlobalRows * numVectors;
698 ApplyTime_ += (time.
wallTime() - startTime);
703 template<
class MatrixType>
706 applyMat (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
707 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
711 A_.
is_null (), std::runtime_error,
"Ifpack2::Relaxation::applyMat: "
712 "The input matrix A is null. Please call setMatrix() with a nonnull "
713 "input matrix, then call compute(), before calling this method.");
715 ! isComputed (), std::runtime_error,
"Ifpack2::Relaxation::applyMat: "
716 "isComputed() must be true before you may call applyMat(). "
717 "Please call compute() before calling this method.");
718 TEUCHOS_TEST_FOR_EXCEPTION(
719 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
720 "Ifpack2::Relaxation::applyMat: X.getNumVectors() = " << X.getNumVectors ()
721 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
722 A_->apply (X, Y, mode);
726 template<
class MatrixType>
729 const char methodName[] =
"Ifpack2::Relaxation::initialize";
732 (A_.
is_null (), std::runtime_error, methodName <<
": The "
733 "input matrix A is null. Please call setMatrix() with "
734 "a nonnull input matrix before calling this method.");
739 double startTime = timer->wallTime();
743 isInitialized_ =
false;
746 auto rowMap = A_->getRowMap ();
747 auto comm = rowMap.
is_null () ? Teuchos::null : rowMap->getComm ();
748 IsParallel_ = ! comm.is_null () && comm->getSize () > 1;
758 Importer_ = IsParallel_ ? A_->getGraph ()->getImporter () :
763 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
764 hasBlockCrsMatrix_ = ! A_bcrs.
is_null ();
767 serialGaussSeidel_ = Teuchos::null;
769 if (PrecType_ == Details::MTGS || PrecType_ == Details::MTSGS ||
770 PrecType_ == Details::GS2 || PrecType_ == Details::SGS2) {
771 auto crsMat = Details::getCrsMatrix(A_);
773 (crsMat.is_null(), std::logic_error, methodName <<
": "
774 "Multithreaded Gauss-Seidel methods currently only work "
775 "when the input matrix is a Tpetra::CrsMatrix.");
780 if (ifpack2_dump_matrix_) {
781 static int sequence_number = 0;
782 const std::string file_name =
"Ifpack2_MT_GS_" +
783 std::to_string (sequence_number++) +
".mtx";
784 using writer_type = Tpetra::MatrixMarket::Writer<crs_matrix_type>;
785 writer_type::writeSparseFile (file_name, crsMat);
788 this->mtKernelHandle_ =
Teuchos::rcp (
new mt_kernel_handle_type ());
789 if (mtKernelHandle_->get_gs_handle () ==
nullptr) {
790 if (PrecType_ == Details::GS2 || PrecType_ == Details::SGS2)
791 mtKernelHandle_->create_gs_handle (KokkosSparse::GS_TWOSTAGE);
792 else if(this->clusterSize_ == 1) {
793 mtKernelHandle_->create_gs_handle (KokkosSparse::GS_DEFAULT, this->mtColoringAlgorithm_);
794 mtKernelHandle_->get_point_gs_handle()->set_long_row_threshold(longRowThreshold_);
797 mtKernelHandle_->create_gs_handle (KokkosSparse::CLUSTER_DEFAULT, this->clusterSize_, this->mtColoringAlgorithm_);
799 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice ();
800 if (PrecType_ == Details::GS2 || PrecType_ == Details::SGS2) {
802 mtKernelHandle_->set_gs_set_num_inner_sweeps (NumInnerSweeps_);
803 mtKernelHandle_->set_gs_set_num_outer_sweeps (NumOuterSweeps_);
804 mtKernelHandle_->set_gs_set_inner_damp_factor (InnerDampingFactor_);
805 mtKernelHandle_->set_gs_twostage (!InnerSpTrsv_, A_->getLocalNumRows ());
806 mtKernelHandle_->set_gs_twostage_compact_form (CompactForm_);
809 KokkosSparse::Experimental::gauss_seidel_symbolic(
810 mtKernelHandle_.getRawPtr (),
811 A_->getLocalNumRows (),
812 A_->getLocalNumCols (),
815 is_matrix_structurally_symmetric_);
819 InitializeTime_ += (timer->wallTime() - startTime);
821 isInitialized_ =
true;
825 template <
typename BlockDiagView>
826 struct InvertDiagBlocks {
827 typedef typename BlockDiagView::size_type Size;
830 typedef Kokkos::MemoryTraits<Kokkos::Unmanaged> Unmanaged;
831 template <
typename View>
832 using UnmanagedView = Kokkos::View<
typename View::data_type,
typename View::array_layout,
833 typename View::device_type, Unmanaged>;
835 typedef typename BlockDiagView::non_const_value_type Scalar;
836 typedef typename BlockDiagView::device_type Device;
837 typedef Kokkos::View<Scalar**, Kokkos::LayoutRight, Device> RWrk;
838 typedef Kokkos::View<int**, Kokkos::LayoutRight, Device> IWrk;
840 UnmanagedView<BlockDiagView> block_diag_;
844 UnmanagedView<RWrk> rwrk_;
846 UnmanagedView<IWrk> iwrk_;
849 InvertDiagBlocks (BlockDiagView& block_diag)
850 : block_diag_(block_diag)
852 const auto blksz = block_diag.extent(1);
853 Kokkos::resize(rwrk_buf_, block_diag_.extent(0), blksz);
855 Kokkos::resize(iwrk_buf_, block_diag_.extent(0), blksz);
859 KOKKOS_INLINE_FUNCTION
860 void operator() (
const Size i,
int& jinfo)
const {
861 auto D_cur = Kokkos::subview(block_diag_, i, Kokkos::ALL(), Kokkos::ALL());
862 auto ipiv = Kokkos::subview(iwrk_, i, Kokkos::ALL());
863 auto work = Kokkos::subview(rwrk_, i, Kokkos::ALL());
865 Tpetra::GETF2(D_cur, ipiv, info);
870 Tpetra::GETRI(D_cur, ipiv, work, info);
876 template<
class MatrixType>
877 void Relaxation<MatrixType>::computeBlockCrs ()
890 using Teuchos::rcp_dynamic_cast;
891 using Teuchos::reduceAll;
892 typedef local_ordinal_type LO;
894 const std::string timerName (
"Ifpack2::Relaxation::computeBlockCrs");
899 double startTime = timer->
wallTime();
904 A_.
is_null (), std::runtime_error,
"Ifpack2::Relaxation::"
905 "computeBlockCrs: The input matrix A is null. Please call setMatrix() "
906 "with a nonnull input matrix, then call initialize(), before calling "
908 auto blockCrsA = rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
910 blockCrsA.is_null(), std::logic_error,
"Ifpack2::Relaxation::"
911 "computeBlockCrs: A_ is not a BlockCrsMatrix, but it should be if we "
912 "got this far. Please report this bug to the Ifpack2 developers.");
914 const scalar_type one = STS::one ();
919 const LO lclNumMeshRows =
920 blockCrsA->getCrsGraph ().getLocalNumRows ();
921 const LO blockSize = blockCrsA->getBlockSize ();
923 if (! savedDiagOffsets_) {
924 blockDiag_ = block_diag_type ();
925 blockDiag_ = block_diag_type (
"Ifpack2::Relaxation::blockDiag_",
926 lclNumMeshRows, blockSize, blockSize);
927 if (Teuchos::as<LO>(diagOffsets_.extent (0) ) < lclNumMeshRows) {
930 diagOffsets_ = Kokkos::View<size_t*, device_type> ();
931 diagOffsets_ = Kokkos::View<size_t*, device_type> (
"offsets", lclNumMeshRows);
933 blockCrsA->getCrsGraph ().getLocalDiagOffsets (diagOffsets_);
935 (static_cast<size_t> (diagOffsets_.extent (0)) !=
936 static_cast<size_t> (blockDiag_.extent (0)),
937 std::logic_error,
"diagOffsets_.extent(0) = " <<
938 diagOffsets_.extent (0) <<
" != blockDiag_.extent(0) = "
939 << blockDiag_.extent (0) <<
940 ". Please report this bug to the Ifpack2 developers.");
941 savedDiagOffsets_ =
true;
943 blockCrsA->getLocalDiagCopy (blockDiag_, diagOffsets_);
950 unmanaged_block_diag_type blockDiag = blockDiag_;
952 if (DoL1Method_ && IsParallel_) {
953 const scalar_type two = one + one;
954 const size_t maxLength = A_->getLocalMaxNumRowEntries ();
955 nonconst_local_inds_host_view_type indices (
"indices",maxLength);
956 nonconst_values_host_view_type values_ (
"values",maxLength * blockSize * blockSize);
957 size_t numEntries = 0;
959 for (LO i = 0; i < lclNumMeshRows; ++i) {
961 blockCrsA->getLocalRowCopy (i, indices, values_, numEntries);
962 scalar_type * values =
reinterpret_cast<scalar_type*
>(values_.data());
964 auto diagBlock = Kokkos::subview (blockDiag, i, ALL (), ALL ());
965 for (LO subRow = 0; subRow < blockSize; ++subRow) {
966 magnitude_type diagonal_boost = STM::zero ();
967 for (
size_t k = 0 ; k < numEntries ; ++k) {
968 if (indices[k] >= lclNumMeshRows) {
969 const size_t offset = blockSize*blockSize*k + subRow*blockSize;
970 for (LO subCol = 0; subCol < blockSize; ++subCol) {
971 diagonal_boost += STS::magnitude (values[offset+subCol] / two);
975 if (STS::magnitude (diagBlock(subRow, subRow)) < L1Eta_ * diagonal_boost) {
976 diagBlock(subRow, subRow) += diagonal_boost;
984 Impl::InvertDiagBlocks<unmanaged_block_diag_type> idb(blockDiag);
985 typedef typename unmanaged_block_diag_type::execution_space exec_space;
986 typedef Kokkos::RangePolicy<exec_space, LO> range_type;
988 Kokkos::parallel_reduce (range_type (0, lclNumMeshRows), idb, info);
994 (info > 0, std::runtime_error,
"GETF2 or GETRI failed on = " << info
995 <<
" diagonal blocks.");
1000 #ifdef HAVE_IFPACK2_DEBUG
1001 const int numResults = 2;
1003 int lclResults[2], gblResults[2];
1004 lclResults[0] = info;
1005 lclResults[1] = -info;
1008 reduceAll<int, int> (* (A_->getGraph ()->getComm ()), REDUCE_MIN,
1009 numResults, lclResults, gblResults);
1011 (gblResults[0] != 0 || gblResults[1] != 0, std::runtime_error,
1012 "Ifpack2::Relaxation::compute: When processing the input "
1013 "Tpetra::BlockCrsMatrix, one or more diagonal block LU factorizations "
1014 "failed on one or more (MPI) processes.");
1015 #endif // HAVE_IFPACK2_DEBUG
1016 serialGaussSeidel_ =
rcp(
new SerialGaussSeidel(blockCrsA, blockDiag_, localSmoothingIndices_, DampingFactor_));
1019 ComputeTime_ += (timer->
wallTime() - startTime);
1024 template<
class MatrixType>
1037 using Teuchos::reduceAll;
1041 using IST =
typename vector_type::impl_scalar_type;
1042 using KAT = Kokkos::ArithTraits<IST>;
1044 const char methodName[] =
"Ifpack2::Relaxation::compute";
1045 const scalar_type zero = STS::zero ();
1046 const scalar_type one = STS::one ();
1049 (A_.
is_null (), std::runtime_error, methodName <<
": "
1050 "The input matrix A is null. Please call setMatrix() with a nonnull "
1051 "input matrix, then call initialize(), before calling this method.");
1054 if (! isInitialized ()) {
1058 if (hasBlockCrsMatrix_) {
1065 double startTime = timer->
wallTime();
1077 IST oneOverMinDiagVal = KAT::one () /
static_cast<IST
> (SmallTraits<scalar_type>::eps ());
1078 if ( MinDiagonalValue_ != zero)
1079 oneOverMinDiagVal = KAT::one () /
static_cast<IST
> (MinDiagonalValue_);
1082 const magnitude_type minDiagValMag = STS::magnitude (MinDiagonalValue_);
1084 const LO numMyRows =
static_cast<LO
> (A_->getLocalNumRows ());
1087 (NumSweeps_ < 0, std::logic_error, methodName
1088 <<
": NumSweeps_ = " << NumSweeps_ <<
" < 0. "
1089 "Please report this bug to the Ifpack2 developers.");
1090 IsComputed_ =
false;
1092 if (Diagonal_.is_null()) {
1095 Diagonal_ =
rcp (
new vector_type (A_->getRowMap (),
false));
1098 if (checkDiagEntries_) {
1104 size_t numSmallDiagEntries = 0;
1105 size_t numZeroDiagEntries = 0;
1106 size_t numNegDiagEntries = 0;
1113 A_->getLocalDiagCopy (*Diagonal_);
1114 std::unique_ptr<vector_type> origDiag;
1115 origDiag = std::unique_ptr<vector_type> (
new vector_type (*Diagonal_,
Teuchos::Copy));
1116 auto diag2d = Diagonal_->getLocalViewHost(Tpetra::Access::ReadWrite);
1117 auto diag = Kokkos::subview(diag2d, Kokkos::ALL(), 0);
1124 if (numMyRows != 0) {
1126 minMagDiagEntryMag = d_0_mag;
1127 maxMagDiagEntryMag = d_0_mag;
1136 for (LO i = 0; i < numMyRows; ++i) {
1137 const IST d_i = diag(i);
1141 const auto d_i_real = getRealValue (d_i);
1145 if (d_i_real < STM::zero ()) {
1146 ++numNegDiagEntries;
1148 if (d_i_mag < minMagDiagEntryMag) {
1149 minMagDiagEntryMag = d_i_mag;
1151 if (d_i_mag > maxMagDiagEntryMag) {
1152 maxMagDiagEntryMag = d_i_mag;
1155 if (fixTinyDiagEntries_) {
1157 if (d_i_mag <= minDiagValMag) {
1158 ++numSmallDiagEntries;
1159 if (d_i_mag == STM::zero ()) {
1160 ++numZeroDiagEntries;
1162 diag(i) = oneOverMinDiagVal;
1165 diag(i) = KAT::one () / d_i;
1170 if (d_i_mag <= minDiagValMag) {
1171 ++numSmallDiagEntries;
1172 if (d_i_mag == STM::zero ()) {
1173 ++numZeroDiagEntries;
1176 diag(i) = KAT::one () / d_i;
1194 const Comm<int>& comm = * (A_->getRowMap ()->getComm ());
1197 Array<magnitude_type> localVals (2);
1198 localVals[0] = minMagDiagEntryMag;
1200 localVals[1] = -maxMagDiagEntryMag;
1201 Array<magnitude_type> globalVals (2);
1202 reduceAll<int, magnitude_type> (comm, REDUCE_MIN, 2,
1203 localVals.getRawPtr (),
1204 globalVals.getRawPtr ());
1205 globalMinMagDiagEntryMag_ = globalVals[0];
1206 globalMaxMagDiagEntryMag_ = -globalVals[1];
1208 Array<size_t> localCounts (3);
1209 localCounts[0] = numSmallDiagEntries;
1210 localCounts[1] = numZeroDiagEntries;
1211 localCounts[2] = numNegDiagEntries;
1212 Array<size_t> globalCounts (3);
1213 reduceAll<int, size_t> (comm, REDUCE_SUM, 3,
1214 localCounts.getRawPtr (),
1215 globalCounts.getRawPtr ());
1216 globalNumSmallDiagEntries_ = globalCounts[0];
1217 globalNumZeroDiagEntries_ = globalCounts[1];
1218 globalNumNegDiagEntries_ = globalCounts[2];
1222 vector_type diff (A_->getRowMap ());
1223 diff.reciprocal (*origDiag);
1224 diff.update (-one, *Diagonal_, one);
1225 globalDiagNormDiff_ = diff.norm2 ();
1232 bool debugAgainstSlowPath =
false;
1234 auto crsMat = Details::getCrsMatrix(A_);
1236 if (crsMat.get() && crsMat->isFillComplete ()) {
1241 if (invDiagKernel_.is_null())
1244 invDiagKernel_->setMatrix(crsMat);
1245 invDiagKernel_->compute(*Diagonal_,
1246 DoL1Method_ && IsParallel_, L1Eta_,
1247 fixTinyDiagEntries_, minDiagValMag);
1248 savedDiagOffsets_ =
true;
1252 #ifdef HAVE_IFPACK2_DEBUG
1253 debugAgainstSlowPath =
true;
1257 if (crsMat.is_null() || ! crsMat->isFillComplete () || debugAgainstSlowPath) {
1264 RCP<vector_type> Diagonal;
1265 if (debugAgainstSlowPath)
1266 Diagonal =
rcp(
new vector_type(A_->getRowMap ()));
1268 Diagonal = Diagonal_;
1270 A_->getLocalDiagCopy (*Diagonal);
1280 if (DoL1Method_ && IsParallel_) {
1282 auto diag = Diagonal->getLocalViewHost(Tpetra::Access::ReadWrite);
1284 const size_t maxLength = A_row.getLocalMaxNumRowEntries ();
1285 nonconst_local_inds_host_view_type indices(
"indices",maxLength);
1286 nonconst_values_host_view_type values(
"values",maxLength);
1289 for (LO i = 0; i < numMyRows; ++i) {
1290 A_row.getLocalRowCopy (i, indices, values, numEntries);
1292 for (
size_t k = 0 ; k < numEntries; ++k) {
1293 if (indices[k] >= numMyRows) {
1294 diagonal_boost += STS::magnitude (values[k] / two);
1297 if (KAT::magnitude (diag(i, 0)) < L1Eta_ * diagonal_boost) {
1298 diag(i, 0) += diagonal_boost;
1306 if (fixTinyDiagEntries_) {
1310 auto localDiag = Diagonal->getLocalViewDevice(Tpetra::Access::ReadWrite);
1311 Kokkos::parallel_for(Kokkos::RangePolicy<MyExecSpace>(0, localDiag.extent(0)),
1312 KOKKOS_LAMBDA (local_ordinal_type i) {
1313 auto d_i = localDiag(i, 0);
1316 if (d_i_mag <= minDiagValMag) {
1317 d_i = oneOverMinDiagVal;
1322 d_i = IST (KAT::one () / d_i);
1324 localDiag(i, 0) = d_i;
1328 Diagonal->reciprocal (*Diagonal);
1331 if (debugAgainstSlowPath) {
1333 Diagonal->update (STS::one (), *Diagonal_, -STS::one ());
1338 (err > 100*STM::eps(), std::logic_error, methodName <<
": "
1339 <<
"\"fast-path\" diagonal computation failed. "
1340 "\\|D1 - D2\\|_inf = " << err <<
".");
1347 if (STS::isComplex) {
1348 ComputeFlops_ += 4.0 * numMyRows;
1351 ComputeFlops_ += numMyRows;
1354 if (PrecType_ == Ifpack2::Details::MTGS ||
1355 PrecType_ == Ifpack2::Details::MTSGS ||
1356 PrecType_ == Ifpack2::Details::GS2 ||
1357 PrecType_ == Ifpack2::Details::SGS2) {
1360 using scalar_view_t =
typename local_matrix_device_type::values_type;
1363 (crsMat.is_null(), std::logic_error, methodName <<
": "
1364 "Multithreaded Gauss-Seidel methods currently only work "
1365 "when the input matrix is a Tpetra::CrsMatrix.");
1366 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice ();
1370 auto diagView_2d = Diagonal_->getLocalViewDevice (Tpetra::Access::ReadWrite);
1371 scalar_view_t diagView_1d = Kokkos::subview (diagView_2d, Kokkos::ALL (), 0);
1372 KokkosSparse::Experimental::gauss_seidel_numeric(
1373 mtKernelHandle_.getRawPtr (),
1374 A_->getLocalNumRows (),
1375 A_->getLocalNumCols (),
1380 is_matrix_structurally_symmetric_);
1382 else if(PrecType_ == Ifpack2::Details::GS || PrecType_ == Ifpack2::Details::SGS) {
1384 serialGaussSeidel_ =
rcp(
new SerialGaussSeidel(crsMat, Diagonal_, localSmoothingIndices_, DampingFactor_));
1386 serialGaussSeidel_ =
rcp(
new SerialGaussSeidel(A_, Diagonal_, localSmoothingIndices_, DampingFactor_));
1390 ComputeTime_ += (timer->
wallTime() - startTime);
1396 template<
class MatrixType>
1399 ApplyInverseRichardson (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1400 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const
1403 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1404 const double numVectors = as<double> (X.getNumVectors ());
1405 if (ZeroStartingSolution_) {
1409 Y.scale(DampingFactor_,X);
1415 double flopUpdate = 0.0;
1416 if (DampingFactor_ == STS::one ()) {
1418 flopUpdate = numGlobalRows * numVectors;
1422 flopUpdate = numGlobalRows * numVectors;
1424 ApplyFlops_ += flopUpdate;
1425 if (NumSweeps_ == 1) {
1431 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1434 updateCachedMultiVector(Y.getMap(),as<size_t>(numVectors));
1436 for (
int j = startSweep; j < NumSweeps_; ++j) {
1438 Tpetra::Details::residual(*A_,Y,X,*cachedMV_);
1439 Y.update(DampingFactor_,*cachedMV_,STS::one());
1453 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1454 const double dampingFlops = (DampingFactor_ == STS::one ()) ? 0.0 : 1.0;
1455 ApplyFlops_ += as<double> (NumSweeps_ - startSweep) * numVectors *
1456 (2.0 * numGlobalNonzeros + dampingFlops);
1460 template<
class MatrixType>
1462 Relaxation<MatrixType>::
1463 ApplyInverseJacobi (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1464 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const
1467 if (hasBlockCrsMatrix_) {
1468 ApplyInverseJacobi_BlockCrsMatrix (X, Y);
1472 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1473 const double numVectors = as<double> (X.getNumVectors ());
1474 if (ZeroStartingSolution_) {
1479 Y.elementWiseMultiply (DampingFactor_, *Diagonal_, X, STS::zero ());
1485 double flopUpdate = 0.0;
1486 if (DampingFactor_ == STS::one ()) {
1488 flopUpdate = numGlobalRows * numVectors;
1492 flopUpdate = 2.0 * numGlobalRows * numVectors;
1494 ApplyFlops_ += flopUpdate;
1495 if (NumSweeps_ == 1) {
1501 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1504 updateCachedMultiVector(Y.getMap(),as<size_t>(numVectors));
1506 for (
int j = startSweep; j < NumSweeps_; ++j) {
1508 Tpetra::Details::residual(*A_,Y,X,*cachedMV_);
1509 Y.elementWiseMultiply (DampingFactor_, *Diagonal_, *cachedMV_, STS::one ());
1523 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1524 const double dampingFlops = (DampingFactor_ == STS::one ()) ? 0.0 : 1.0;
1525 ApplyFlops_ += as<double> (NumSweeps_ - startSweep) * numVectors *
1526 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1529 template<
class MatrixType>
1531 Relaxation<MatrixType>::
1532 ApplyInverseJacobi_BlockCrsMatrix (
const Tpetra::MultiVector<scalar_type,
1534 global_ordinal_type,
1536 Tpetra::MultiVector<scalar_type,
1538 global_ordinal_type,
1539 node_type>& Y)
const
1541 using Tpetra::BlockMultiVector;
1542 using BMV = BlockMultiVector<scalar_type, local_ordinal_type,
1543 global_ordinal_type, node_type>;
1545 const block_crs_matrix_type* blockMatConst =
1546 dynamic_cast<const block_crs_matrix_type*
> (A_.
getRawPtr ());
1548 (blockMatConst ==
nullptr, std::logic_error,
"This method should "
1549 "never be called if the matrix A_ is not a BlockCrsMatrix. "
1550 "Please report this bug to the Ifpack2 developers.");
1555 block_crs_matrix_type* blockMat =
1556 const_cast<block_crs_matrix_type*
> (blockMatConst);
1558 auto meshRowMap = blockMat->getRowMap ();
1559 auto meshColMap = blockMat->getColMap ();
1560 const local_ordinal_type blockSize = blockMat->getBlockSize ();
1562 BMV X_blk (X, *meshColMap, blockSize);
1563 BMV Y_blk (Y, *meshRowMap, blockSize);
1565 if (ZeroStartingSolution_) {
1569 Y_blk.blockWiseMultiply (DampingFactor_, blockDiag_, X_blk);
1570 if (NumSweeps_ == 1) {
1575 auto pointRowMap = Y.getMap ();
1576 const size_t numVecs = X.getNumVectors ();
1580 BMV A_times_Y (*meshRowMap, *pointRowMap, blockSize, numVecs);
1584 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1586 for (
int j = startSweep; j < NumSweeps_; ++j) {
1587 blockMat->applyBlock (Y_blk, A_times_Y);
1589 Y_blk.blockJacobiUpdate (DampingFactor_, blockDiag_,
1590 X_blk, A_times_Y, STS::one ());
1594 template<
class MatrixType>
1596 Relaxation<MatrixType>::
1597 ApplyInverseSerialGS (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1598 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1599 Tpetra::ESweepDirection direction)
const
1601 using this_type = Relaxation<MatrixType>;
1611 auto blockCrsMat = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (A_);
1612 auto crsMat = Details::getCrsMatrix(A_);
1613 if (blockCrsMat.get()) {
1614 const_cast<this_type&
> (*this).ApplyInverseSerialGS_BlockCrsMatrix (*blockCrsMat, X, Y, direction);
1616 else if (crsMat.get()) {
1617 ApplyInverseSerialGS_CrsMatrix (*crsMat, X, Y, direction);
1620 ApplyInverseSerialGS_RowMatrix (X, Y, direction);
1625 template<
class MatrixType>
1627 Relaxation<MatrixType>::
1628 ApplyInverseSerialGS_RowMatrix (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1629 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1630 Tpetra::ESweepDirection direction)
const {
1637 using Teuchos::rcpFromRef;
1638 typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
1643 if (ZeroStartingSolution_) {
1644 Y.putScalar (STS::zero ());
1647 size_t NumVectors = X.getNumVectors();
1651 if (Importer_.is_null ()) {
1652 updateCachedMultiVector (Y.getMap (), NumVectors);
1655 updateCachedMultiVector (Importer_->getTargetMap (), NumVectors);
1660 Y2 = rcpFromRef (Y);
1663 for (
int j = 0; j < NumSweeps_; ++j) {
1666 if (Importer_.is_null ()) {
1672 Y2->doImport (Y, *Importer_, Tpetra::INSERT);
1675 serialGaussSeidel_->apply(*Y2, X, direction);
1679 Tpetra::deep_copy (Y, *Y2);
1684 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1685 const double numVectors = as<double> (X.getNumVectors ());
1686 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1687 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1688 ApplyFlops_ += 2.0 * NumSweeps_ * numVectors *
1689 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1693 template<
class MatrixType>
1695 Relaxation<MatrixType>::
1696 ApplyInverseSerialGS_CrsMatrix(
const crs_matrix_type& A,
1697 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& B,
1698 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1699 Tpetra::ESweepDirection direction)
const
1701 using Teuchos::null;
1704 using Teuchos::rcpFromRef;
1705 using Teuchos::rcp_const_cast;
1706 typedef scalar_type Scalar;
1707 const char prefix[] =
"Ifpack2::Relaxation::SerialGS: ";
1711 ! A.isFillComplete (), std::runtime_error,
1712 prefix <<
"The matrix is not fill complete.");
1714 NumSweeps_ < 0, std::invalid_argument,
1715 prefix <<
"The number of sweeps must be nonnegative, "
1716 "but you provided numSweeps = " << NumSweeps_ <<
" < 0.");
1718 if (NumSweeps_ == 0) {
1722 RCP<const import_type> importer = A.getGraph ()->getImporter ();
1724 RCP<const map_type> domainMap = A.getDomainMap ();
1725 RCP<const map_type> rangeMap = A.getRangeMap ();
1726 RCP<const map_type> rowMap = A.getGraph ()->getRowMap ();
1727 RCP<const map_type> colMap = A.getGraph ()->getColMap ();
1729 #ifdef HAVE_IFPACK2_DEBUG
1735 ! X.getMap ()->isSameAs (*domainMap), std::runtime_error,
1736 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1737 "multivector X be in the domain Map of the matrix.");
1739 ! B.getMap ()->isSameAs (*rangeMap), std::runtime_error,
1740 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1741 "B be in the range Map of the matrix.");
1743 ! Diagonal_->getMap ()->isSameAs (*rowMap), std::runtime_error,
1744 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1745 "D be in the row Map of the matrix.");
1747 ! rowMap->isSameAs (*rangeMap), std::runtime_error,
1748 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the row Map and the "
1749 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
1751 ! domainMap->isSameAs (*rangeMap), std::runtime_error,
1752 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the domain Map and "
1753 "the range Map of the matrix be the same.");
1765 RCP<multivector_type> X_colMap;
1766 RCP<multivector_type> X_domainMap;
1767 bool copyBackOutput =
false;
1768 if (importer.is_null ()) {
1769 X_colMap = Teuchos::rcpFromRef (X);
1770 X_domainMap = Teuchos::rcpFromRef (X);
1776 if (ZeroStartingSolution_) {
1777 X_colMap->putScalar (ZERO);
1782 updateCachedMultiVector(colMap, X.getNumVectors());
1783 X_colMap = cachedMV_;
1784 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
1786 if (ZeroStartingSolution_) {
1788 X_colMap->putScalar (ZERO);
1798 X_colMap->doImport (X, *importer, Tpetra::INSERT);
1800 copyBackOutput =
true;
1803 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
1804 if (! importer.is_null () && sweep > 0) {
1807 X_colMap->doImport (*X_domainMap, *importer, Tpetra::INSERT);
1810 serialGaussSeidel_->apply(*X_colMap, B, direction);
1813 if (copyBackOutput) {
1815 deep_copy (X , *X_domainMap);
1816 }
catch (std::exception& e) {
1818 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) "
1819 "threw an exception: " << e.what ());
1834 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1835 const double numVectors = X.getNumVectors ();
1836 const double numGlobalRows = A_->getGlobalNumRows ();
1837 const double numGlobalNonzeros = A_->getGlobalNumEntries ();
1838 ApplyFlops_ += NumSweeps_ * numVectors *
1839 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1842 template<
class MatrixType>
1844 Relaxation<MatrixType>::
1845 ApplyInverseSerialGS_BlockCrsMatrix (
const block_crs_matrix_type& A,
1846 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1847 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1848 Tpetra::ESweepDirection direction)
1850 using Tpetra::INSERT;
1853 using Teuchos::rcpFromRef;
1862 block_multivector_type yBlock(Y, *(A.getGraph ()->getDomainMap()), A.getBlockSize());
1863 const block_multivector_type xBlock(X, *(A.getColMap ()), A.getBlockSize());
1865 bool performImport =
false;
1866 RCP<block_multivector_type> yBlockCol;
1867 if (Importer_.is_null()) {
1868 yBlockCol = rcpFromRef(yBlock);
1871 if (yBlockColumnPointMap_.is_null() ||
1872 yBlockColumnPointMap_->getNumVectors() != yBlock.getNumVectors() ||
1873 yBlockColumnPointMap_->getBlockSize() != yBlock.getBlockSize()) {
1874 yBlockColumnPointMap_ =
1875 rcp(
new block_multivector_type(*(A.getColMap()), A.getBlockSize(),
1876 static_cast<local_ordinal_type
>(yBlock.getNumVectors())));
1878 yBlockCol = yBlockColumnPointMap_;
1879 if (pointImporter_.is_null()) {
1880 auto srcMap =
rcp(
new map_type(yBlock.getPointMap()));
1881 auto tgtMap =
rcp(
new map_type(yBlockCol->getPointMap()));
1882 pointImporter_ =
rcp(
new import_type(srcMap, tgtMap));
1884 performImport =
true;
1887 multivector_type yBlock_mv;
1888 multivector_type yBlockCol_mv;
1889 RCP<const multivector_type> yBlockColPointDomain;
1890 if (performImport) {
1891 yBlock_mv = yBlock.getMultiVectorView();
1892 yBlockCol_mv = yBlockCol->getMultiVectorView();
1893 yBlockColPointDomain =
1894 yBlockCol_mv.offsetView(A.getDomainMap(), 0);
1897 if (ZeroStartingSolution_) {
1898 yBlockCol->putScalar(STS::zero ());
1900 else if (performImport) {
1901 yBlockCol_mv.doImport(yBlock_mv, *pointImporter_, Tpetra::INSERT);
1904 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
1905 if (performImport && sweep > 0) {
1906 yBlockCol_mv.doImport(yBlock_mv, *pointImporter_, Tpetra::INSERT);
1908 serialGaussSeidel_->applyBlock(*yBlockCol, xBlock, direction);
1909 if (performImport) {
1910 Tpetra::deep_copy(Y, *yBlockColPointDomain);
1915 template<
class MatrixType>
1917 Relaxation<MatrixType>::
1918 ApplyInverseMTGS_CrsMatrix(
1919 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& B,
1920 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1921 Tpetra::ESweepDirection direction)
const
1923 using Teuchos::null;
1926 using Teuchos::rcpFromRef;
1927 using Teuchos::rcp_const_cast;
1930 typedef scalar_type Scalar;
1932 const char prefix[] =
"Ifpack2::Relaxation::(reordered)MTGaussSeidel: ";
1935 auto crsMat = Details::getCrsMatrix(A_);
1937 (crsMat.is_null(), std::logic_error,
"Ifpack2::Relaxation::apply: "
1938 "Multithreaded Gauss-Seidel methods currently only work when the "
1939 "input matrix is a Tpetra::CrsMatrix.");
1943 (! localSmoothingIndices_.is_null (), std::logic_error,
1944 "Ifpack2's implementation of Multithreaded Gauss-Seidel does not "
1945 "implement the case where the user supplies an iteration order. "
1946 "This error used to appear as \"MT GaussSeidel ignores the given "
1948 "I tried to add more explanation, but I didn't implement \"MT "
1949 "GaussSeidel\" [sic]. "
1950 "You'll have to ask the person who did.");
1953 (! crsMat->isFillComplete (), std::runtime_error, prefix <<
"The "
1954 "input CrsMatrix is not fill complete. Please call fillComplete "
1955 "on the matrix, then do setup again, before calling apply(). "
1956 "\"Do setup\" means that if only the matrix's values have changed "
1957 "since last setup, you need only call compute(). If the matrix's "
1958 "structure may also have changed, you must first call initialize(), "
1959 "then call compute(). If you have not set up this preconditioner "
1960 "for this matrix before, you must first call initialize(), then "
1963 (NumSweeps_ < 0, std::logic_error, prefix <<
": NumSweeps_ = "
1964 << NumSweeps_ <<
" < 0. Please report this bug to the Ifpack2 "
1966 if (NumSweeps_ == 0) {
1970 RCP<const import_type> importer = crsMat->getGraph ()->getImporter ();
1972 ! crsMat->getGraph ()->getExporter ().is_null (), std::runtime_error,
1973 "This method's implementation currently requires that the matrix's row, "
1974 "domain, and range Maps be the same. This cannot be the case, because "
1975 "the matrix has a nontrivial Export object.");
1977 RCP<const map_type> domainMap = crsMat->getDomainMap ();
1978 RCP<const map_type> rangeMap = crsMat->getRangeMap ();
1979 RCP<const map_type> rowMap = crsMat->getGraph ()->getRowMap ();
1980 RCP<const map_type> colMap = crsMat->getGraph ()->getColMap ();
1982 #ifdef HAVE_IFPACK2_DEBUG
1988 ! X.getMap ()->isSameAs (*domainMap), std::runtime_error,
1989 "Ifpack2::Relaxation::MTGaussSeidel requires that the input "
1990 "multivector X be in the domain Map of the matrix.");
1992 ! B.getMap ()->isSameAs (*rangeMap), std::runtime_error,
1993 "Ifpack2::Relaxation::MTGaussSeidel requires that the input "
1994 "B be in the range Map of the matrix.");
1996 ! rowMap->isSameAs (*rangeMap), std::runtime_error,
1997 "Ifpack2::Relaxation::MTGaussSeidel requires that the row Map and the "
1998 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
2000 ! domainMap->isSameAs (*rangeMap), std::runtime_error,
2001 "Ifpack2::Relaxation::MTGaussSeidel requires that the domain Map and "
2002 "the range Map of the matrix be the same.");
2008 #endif // HAVE_IFPACK2_DEBUG
2018 RCP<multivector_type> X_colMap;
2019 RCP<multivector_type> X_domainMap;
2020 bool copyBackOutput =
false;
2021 if (importer.is_null ()) {
2022 if (X.isConstantStride ()) {
2023 X_colMap = rcpFromRef (X);
2024 X_domainMap = rcpFromRef (X);
2031 if (ZeroStartingSolution_) {
2032 X_colMap->putScalar (ZERO);
2041 updateCachedMultiVector(colMap,as<size_t>(X.getNumVectors()));
2042 X_colMap = cachedMV_;
2046 X_domainMap = X_colMap;
2047 if (! ZeroStartingSolution_) {
2049 deep_copy (*X_domainMap , X);
2050 }
catch (std::exception& e) {
2051 std::ostringstream os;
2052 os <<
"Ifpack2::Relaxation::MTGaussSeidel: "
2053 "deep_copy(*X_domainMap, X) threw an exception: "
2054 << e.what () <<
".";
2058 copyBackOutput =
true;
2072 updateCachedMultiVector(colMap,as<size_t>(X.getNumVectors()));
2073 X_colMap = cachedMV_;
2075 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
2077 if (ZeroStartingSolution_) {
2079 X_colMap->putScalar (ZERO);
2089 X_colMap->doImport (X, *importer, Tpetra::CombineMode::INSERT);
2091 copyBackOutput =
true;
2097 RCP<const multivector_type> B_in;
2098 if (B.isConstantStride ()) {
2099 B_in = rcpFromRef (B);
2105 RCP<multivector_type> B_in_nonconst =
rcp (
new multivector_type(rowMap, B.getNumVectors()));
2107 deep_copy (*B_in_nonconst, B);
2108 }
catch (std::exception& e) {
2109 std::ostringstream os;
2110 os <<
"Ifpack2::Relaxation::MTGaussSeidel: "
2111 "deep_copy(*B_in_nonconst, B) threw an exception: "
2112 << e.what () <<
".";
2115 B_in = rcp_const_cast<
const multivector_type> (B_in_nonconst);
2129 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice ();
2131 bool update_y_vector =
true;
2133 bool zero_x_vector =
false;
2135 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
2136 if (! importer.is_null () && sweep > 0) {
2139 X_colMap->doImport (*X_domainMap, *importer, Tpetra::CombineMode::INSERT);
2142 if (direction == Tpetra::Symmetric) {
2143 KokkosSparse::Experimental::symmetric_gauss_seidel_apply
2144 (mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2145 kcsr.graph.row_map, kcsr.graph.entries, kcsr.values,
2146 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2147 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2148 zero_x_vector, update_y_vector, DampingFactor_, 1);
2150 else if (direction == Tpetra::Forward) {
2151 KokkosSparse::Experimental::forward_sweep_gauss_seidel_apply
2152 (mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2153 kcsr.graph.row_map,kcsr.graph.entries, kcsr.values,
2154 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2155 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2156 zero_x_vector, update_y_vector, DampingFactor_, 1);
2158 else if (direction == Tpetra::Backward) {
2159 KokkosSparse::Experimental::backward_sweep_gauss_seidel_apply
2160 (mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2161 kcsr.graph.row_map,kcsr.graph.entries, kcsr.values,
2162 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2163 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2164 zero_x_vector, update_y_vector, DampingFactor_, 1);
2168 true, std::invalid_argument,
2169 prefix <<
"The 'direction' enum does not have any of its valid "
2170 "values: Forward, Backward, or Symmetric.");
2172 update_y_vector =
false;
2175 if (copyBackOutput) {
2177 deep_copy (X , *X_domainMap);
2178 }
catch (std::exception& e) {
2180 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) "
2181 "threw an exception: " << e.what ());
2185 const double dampingFlops = (DampingFactor_ == STS::one ()) ? 0.0 : 1.0;
2186 const double numVectors = as<double> (X.getNumVectors ());
2187 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
2188 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
2189 double ApplyFlops = NumSweeps_ * numVectors *
2190 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
2191 if (direction == Tpetra::Symmetric)
2193 ApplyFlops_ += ApplyFlops;
2196 template<
class MatrixType>
2199 std::ostringstream os;
2204 os <<
"\"Ifpack2::Relaxation\": {";
2206 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", "
2207 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
2213 if (PrecType_ == Ifpack2::Details::JACOBI) {
2215 }
else if (PrecType_ == Ifpack2::Details::GS) {
2216 os <<
"Gauss-Seidel";
2217 }
else if (PrecType_ == Ifpack2::Details::SGS) {
2218 os <<
"Symmetric Gauss-Seidel";
2219 }
else if (PrecType_ == Ifpack2::Details::MTGS) {
2220 os <<
"MT Gauss-Seidel";
2221 }
else if (PrecType_ == Ifpack2::Details::MTSGS) {
2222 os <<
"MT Symmetric Gauss-Seidel";
2223 }
else if (PrecType_ == Ifpack2::Details::GS2) {
2224 os <<
"Two-stage Gauss-Seidel";
2225 }
else if (PrecType_ == Ifpack2::Details::SGS2) {
2226 os <<
"Two-stage Symmetric Gauss-Seidel";
2231 if(hasBlockCrsMatrix_)
2234 os <<
", " <<
"sweeps: " << NumSweeps_ <<
", "
2235 <<
"damping factor: " << DampingFactor_ <<
", ";
2237 if (PrecType_ == Ifpack2::Details::MTGS || PrecType_ == Ifpack2::Details::MTSGS) {
2238 os <<
"\"relaxation: mtgs cluster size\": " << clusterSize_ <<
", ";
2239 os <<
"\"relaxation: long row threshold\": " << longRowThreshold_ <<
", ";
2240 os <<
"\"relaxation: symmetric matrix structure\": " << (is_matrix_structurally_symmetric_ ?
"true" :
"false") <<
", ";
2241 os <<
"\"relaxation: relaxation: mtgs coloring algorithm\": ";
2242 switch(mtColoringAlgorithm_)
2244 case KokkosGraph::COLORING_DEFAULT:
2245 os <<
"DEFAULT";
break;
2246 case KokkosGraph::COLORING_SERIAL:
2247 os <<
"SERIAL";
break;
2248 case KokkosGraph::COLORING_VB:
2250 case KokkosGraph::COLORING_VBBIT:
2251 os <<
"VBBIT";
break;
2252 case KokkosGraph::COLORING_VBCS:
2253 os <<
"VBCS";
break;
2254 case KokkosGraph::COLORING_VBD:
2256 case KokkosGraph::COLORING_VBDBIT:
2257 os <<
"VBDBIT";
break;
2258 case KokkosGraph::COLORING_EB:
2266 if (PrecType_ == Ifpack2::Details::GS2 ||
2267 PrecType_ == Ifpack2::Details::SGS2) {
2268 os <<
"outer sweeps: " << NumOuterSweeps_ <<
", "
2269 <<
"inner sweeps: " << NumInnerSweeps_ <<
", "
2270 <<
"inner damping factor: " << InnerDampingFactor_ <<
", ";
2274 os <<
"use l1: " << DoL1Method_ <<
", "
2275 <<
"l1 eta: " << L1Eta_ <<
", ";
2279 os <<
"Matrix: null";
2282 os <<
"Global matrix dimensions: ["
2283 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]"
2284 <<
", Global nnz: " << A_->getGlobalNumEntries();
2292 template<
class MatrixType>
2311 const int myRank = this->getComm ()->getRank ();
2323 out <<
"\"Ifpack2::Relaxation\":" << endl;
2325 out <<
"MatrixType: \"" << TypeNameTraits<MatrixType>::name () <<
"\""
2327 if (this->getObjectLabel () !=
"") {
2328 out <<
"Label: " << this->getObjectLabel () << endl;
2330 out <<
"Initialized: " << (isInitialized () ?
"true" :
"false") << endl
2331 <<
"Computed: " << (isComputed () ?
"true" :
"false") << endl
2332 <<
"Parameters: " << endl;
2335 out <<
"\"relaxation: type\": ";
2336 if (PrecType_ == Ifpack2::Details::JACOBI) {
2338 }
else if (PrecType_ == Ifpack2::Details::GS) {
2339 out <<
"Gauss-Seidel";
2340 }
else if (PrecType_ == Ifpack2::Details::SGS) {
2341 out <<
"Symmetric Gauss-Seidel";
2342 }
else if (PrecType_ == Ifpack2::Details::MTGS) {
2343 out <<
"MT Gauss-Seidel";
2344 }
else if (PrecType_ == Ifpack2::Details::MTSGS) {
2345 out <<
"MT Symmetric Gauss-Seidel";
2346 }
else if (PrecType_ == Ifpack2::Details::GS2) {
2347 out <<
"Two-stage Gauss-Seidel";
2348 }
else if (PrecType_ == Ifpack2::Details::SGS2) {
2349 out <<
"Two-stage Symmetric Gauss-Seidel";
2356 <<
"\"relaxation: sweeps\": " << NumSweeps_ << endl
2357 <<
"\"relaxation: damping factor\": " << DampingFactor_ << endl
2358 <<
"\"relaxation: min diagonal value\": " << MinDiagonalValue_ << endl
2359 <<
"\"relaxation: zero starting solution\": " << ZeroStartingSolution_ << endl
2360 <<
"\"relaxation: backward mode\": " << DoBackwardGS_ << endl
2361 <<
"\"relaxation: use l1\": " << DoL1Method_ << endl
2362 <<
"\"relaxation: l1 eta\": " << L1Eta_ << endl;
2363 if (PrecType_ == Ifpack2::Details::MTGS || PrecType_ == Ifpack2::Details::MTSGS) {
2364 out <<
"\"relaxation: mtgs cluster size\": " << clusterSize_ << endl;
2365 out <<
"\"relaxation: long row threshold\": " << longRowThreshold_ << endl;
2366 out <<
"\"relaxation: symmetric matrix structure\": " << (is_matrix_structurally_symmetric_ ?
"true" :
"false") << endl;
2367 out <<
"\"relaxation: relaxation: mtgs coloring algorithm\": ";
2368 switch(mtColoringAlgorithm_)
2370 case KokkosGraph::COLORING_DEFAULT:
2371 out <<
"DEFAULT";
break;
2372 case KokkosGraph::COLORING_SERIAL:
2373 out <<
"SERIAL";
break;
2374 case KokkosGraph::COLORING_VB:
2376 case KokkosGraph::COLORING_VBBIT:
2377 out <<
"VBBIT";
break;
2378 case KokkosGraph::COLORING_VBCS:
2379 out <<
"VBCS";
break;
2380 case KokkosGraph::COLORING_VBD:
2381 out <<
"VBD";
break;
2382 case KokkosGraph::COLORING_VBDBIT:
2383 out <<
"VBDBIT";
break;
2384 case KokkosGraph::COLORING_EB:
2391 if (PrecType_ == Ifpack2::Details::GS2 || PrecType_ == Ifpack2::Details::SGS2) {
2392 out <<
"\"relaxation: inner damping factor\": " << InnerDampingFactor_ << endl;
2393 out <<
"\"relaxation: outer sweeps\" : " << NumOuterSweeps_ << endl;
2394 out <<
"\"relaxation: inner sweeps\" : " << NumInnerSweeps_ << endl;
2397 out <<
"Computed quantities:" << endl;
2400 out <<
"Global number of rows: " << A_->getGlobalNumRows () << endl
2401 <<
"Global number of columns: " << A_->getGlobalNumCols () << endl;
2403 if (checkDiagEntries_ && isComputed ()) {
2404 out <<
"Properties of input diagonal entries:" << endl;
2407 out <<
"Magnitude of minimum-magnitude diagonal entry: "
2408 << globalMinMagDiagEntryMag_ << endl
2409 <<
"Magnitude of maximum-magnitude diagonal entry: "
2410 << globalMaxMagDiagEntryMag_ << endl
2411 <<
"Number of diagonal entries with small magnitude: "
2412 << globalNumSmallDiagEntries_ << endl
2413 <<
"Number of zero diagonal entries: "
2414 << globalNumZeroDiagEntries_ << endl
2415 <<
"Number of diagonal entries with negative real part: "
2416 << globalNumNegDiagEntries_ << endl
2417 <<
"Abs 2-norm diff between computed and actual inverse "
2418 <<
"diagonal: " << globalDiagNormDiff_ << endl;
2421 if (isComputed ()) {
2422 out <<
"Saved diagonal offsets: "
2423 << (savedDiagOffsets_ ?
"true" :
"false") << endl;
2425 out <<
"Call counts and total times (in seconds): " << endl;
2428 out <<
"initialize: " << endl;
2431 out <<
"Call count: " << NumInitialize_ << endl;
2432 out <<
"Total time: " << InitializeTime_ << endl;
2434 out <<
"compute: " << endl;
2437 out <<
"Call count: " << NumCompute_ << endl;
2438 out <<
"Total time: " << ComputeTime_ << endl;
2440 out <<
"apply: " << endl;
2443 out <<
"Call count: " << NumApply_ << endl;
2444 out <<
"Total time: " << ApplyTime_ << endl;
2453 #define IFPACK2_RELAXATION_INSTANT(S,LO,GO,N) \
2454 template class Ifpack2::Relaxation< Tpetra::RowMatrix<S, LO, GO, N> >;
2456 #endif // IFPACK2_RELAXATION_DEF_HPP
bool hasTransposeApply() const
Whether apply() and applyMat() let you apply the transpose or conjugate transpose.
Definition: Ifpack2_Relaxation_def.hpp:525
double getComputeFlops() const
Total number of floating-point operations over all calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:567
double getApplyFlops() const
Total number of floating-point operations over all calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:573
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:2295
bool is_null(const boost::shared_ptr< T > &p)
basic_OSTab< char > OSTab
static magnitudeType eps()
Relaxation(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_Relaxation_def.hpp:214
T & get(const std::string &name, T def_value)
void compute()
Compute the preconditioner ("numeric setup");.
Definition: Ifpack2_Relaxation_def.hpp:1025
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual void printDoc(std::string const &docString, std::ostream &out) const =0
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_Relaxation_def.hpp:580
virtual const std::string getXMLTypeName() const =0
static std::ostream & printLines(std::ostream &os, const std::string &linePrefix, const std::string &lines)
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix and vectors are distributed.
Definition: Ifpack2_Relaxation_def.hpp:482
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Relaxation_def.hpp:193
int getNumApply() const
Total number of calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:543
bool remove(std::string const &name, bool throwIfNotExists=true)
double getComputeTime() const
Total time in seconds spent in all calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:555
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:263
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:502
void setParameters(const Teuchos::ParameterList ¶ms)
Set the relaxation / preconditioner parameters.
Definition: Ifpack2_Relaxation_def.hpp:472
File for utility functions.
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_Relaxation_def.hpp:2197
int getNumInitialize() const
Total number of calls to initialize().
Definition: Ifpack2_Relaxation_def.hpp:531
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
virtual ValidStringsList validStringValues() const =0
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:515
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:260
double getInitializeTime() const
Total time in seconds spent in all calls to initialize().
Definition: Ifpack2_Relaxation_def.hpp:549
any & getAny(bool activeQry=true)
Compute scaled damped residual for Chebyshev.
Definition: Ifpack2_Details_InverseDiagonalKernel_decl.hpp:77
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 input matrix to X, returning the result in Y.
Definition: Ifpack2_Relaxation_def.hpp:706
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix to be preconditioned.
Definition: Ifpack2_Relaxation_def.hpp:493
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:257
bool isType(const std::string &name) const
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:254
std::string typeName() const
const std::type_info & type() const
int getNumCompute() const
Total number of calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:537
Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices.
Definition: Ifpack2_Relaxation_decl.hpp:237
double getApplyTime() const
Total time in seconds spent in all calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:561
void getValidParameters(Teuchos::ParameterList ¶ms)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:51
void initialize()
Initialize the preconditioner ("symbolic setup").
Definition: Ifpack2_Relaxation_def.hpp:727
Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > row_matrix_type
Tpetra::RowMatrix specialization used by this class.
Definition: Ifpack2_Relaxation_decl.hpp:273
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:593
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a list of all the parameters that this class accepts.
Definition: Ifpack2_Relaxation_def.hpp:226
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_Relaxation_decl.hpp:269
virtual void validate(ParameterEntry const &entry, std::string const ¶mName, std::string const &sublistName) const =0