10 #ifndef IFPACK2_RELAXATION_DEF_HPP
11 #define IFPACK2_RELAXATION_DEF_HPP
13 #include "Teuchos_StandardParameterEntryValidators.hpp"
15 #include "Tpetra_CrsMatrix.hpp"
16 #include "Tpetra_BlockCrsMatrix.hpp"
17 #include "Tpetra_BlockView.hpp"
19 #include "Ifpack2_Details_getCrsMatrix.hpp"
20 #include "MatrixMarket_Tpetra.hpp"
21 #include "Tpetra_Details_residual.hpp"
25 #include "KokkosSparse_gauss_seidel.hpp"
32 NonnegativeIntValidator() {}
42 const std::string& paramName,
43 const std::string& sublistName)
const {
49 anyVal.
type() !=
typeid(int),
51 "Parameter \"" << paramName <<
"\" in sublist \"" << sublistName
52 <<
"\" has the wrong type." << endl
53 <<
"Parameter: " << paramName
55 <<
"Type specified: " << entryName << endl
56 <<
"Type required: int" << endl);
58 const int val = Teuchos::any_cast<
int>(anyVal);
61 "Parameter \"" << paramName <<
"\" in sublist \"" << sublistName
62 <<
"\" is negative." << endl
63 <<
"Parameter: " << paramName
65 <<
"Value specified: " << val << endl
66 <<
"Required range: [0, INT_MAX]" << endl);
71 return "NonnegativeIntValidator";
76 printDoc(
const std::string& docString,
77 std::ostream& out)
const {
79 out <<
"#\tValidator Used: " << std::endl;
80 out <<
"#\t\tNonnegativeIntValidator" << std::endl;
87 template <class Scalar, const bool isOrdinal = Teuchos::ScalarTraits<Scalar>::isOrdinal>
91 static const Scalar eps();
95 template <
class Scalar>
96 class SmallTraits<Scalar, true> {
98 static const Scalar eps() {
104 template <
class Scalar>
105 class SmallTraits<Scalar, false> {
107 static const Scalar eps() {
113 template <
class Scalar,
115 struct RealTraits {};
117 template <
class Scalar>
118 struct RealTraits<Scalar, false> {
119 using val_type = Scalar;
120 using mag_type = Scalar;
121 static KOKKOS_INLINE_FUNCTION mag_type real(
const val_type& z) {
126 template <
class Scalar>
127 struct RealTraits<Scalar, true> {
128 using val_type = Scalar;
130 static KOKKOS_INLINE_FUNCTION mag_type real(
const val_type& z) {
131 #if KOKKOS_VERSION >= 40799
132 return KokkosKernels::ArithTraits<val_type>::real(z);
134 return Kokkos::ArithTraits<val_type>::real(z);
139 template <
class Scalar>
140 KOKKOS_INLINE_FUNCTION
typename RealTraits<Scalar>::mag_type
141 getRealValue(
const Scalar& z) {
142 return RealTraits<Scalar>::real(z);
149 template <
class MatrixType>
150 void Relaxation<MatrixType>::
151 updateCachedMultiVector(
const Teuchos::RCP<
const Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>>& map,
152 size_t numVecs)
const {
155 if (cachedMV_.is_null() ||
156 map.get() != cachedMV_->getMap().get() ||
157 cachedMV_->getNumVectors() != numVecs) {
158 using MV = Tpetra::MultiVector<scalar_type, local_ordinal_type,
159 global_ordinal_type, node_type>;
164 template <
class MatrixType>
168 Importer_ = Teuchos::null;
169 pointImporter_ = Teuchos::null;
170 Diagonal_ = Teuchos::null;
171 isInitialized_ =
false;
173 diagOffsets_ = Kokkos::View<size_t*, device_type>();
174 savedDiagOffsets_ =
false;
175 hasBlockCrsMatrix_ =
false;
176 serialGaussSeidel_ = Teuchos::null;
178 IsParallel_ = (A->getRowMap()->getComm()->getSize() > 1);
184 template <
class MatrixType>
188 , IsParallel_((A.
is_null() || A->getRowMap().
is_null() || A->getRowMap()->getComm().
is_null()) ? false :
189 A->getRowMap()->getComm()->getSize() > 1) {
190 this->setObjectLabel(
"Ifpack2::Relaxation");
193 template <
class MatrixType>
198 using Teuchos::parameterList;
201 using Teuchos::rcp_const_cast;
202 using Teuchos::rcp_implicit_cast;
203 using Teuchos::setStringToIntegralParameter;
206 if (validParams_.is_null()) {
207 RCP<ParameterList> pl = parameterList(
"Ifpack2::Relaxation");
211 Array<std::string> precTypes(8);
212 precTypes[0] =
"Jacobi";
213 precTypes[1] =
"Gauss-Seidel";
214 precTypes[2] =
"Symmetric Gauss-Seidel";
215 precTypes[3] =
"MT Gauss-Seidel";
216 precTypes[4] =
"MT Symmetric Gauss-Seidel";
217 precTypes[5] =
"Richardson";
218 precTypes[6] =
"Two-stage Gauss-Seidel";
219 precTypes[7] =
"Two-stage Symmetric Gauss-Seidel";
220 Array<Details::RelaxationType> precTypeEnums(8);
221 precTypeEnums[0] = Details::JACOBI;
222 precTypeEnums[1] = Details::GS;
223 precTypeEnums[2] = Details::SGS;
224 precTypeEnums[3] = Details::MTGS;
225 precTypeEnums[4] = Details::MTSGS;
226 precTypeEnums[5] = Details::RICHARDSON;
227 precTypeEnums[6] = Details::GS2;
228 precTypeEnums[7] = Details::SGS2;
229 const std::string defaultPrecType(
"Jacobi");
230 setStringToIntegralParameter<Details::RelaxationType>(
"relaxation: type",
231 defaultPrecType,
"Relaxation method", precTypes(), precTypeEnums(),
234 const int numSweeps = 1;
235 RCP<PEV> numSweepsValidator =
236 rcp_implicit_cast<PEV>(
rcp(
new NonnegativeIntValidator));
237 pl->set(
"relaxation: sweeps", numSweeps,
"Number of relaxation sweeps",
238 rcp_const_cast<const PEV>(numSweepsValidator));
241 const int numOuterSweeps = 1;
242 RCP<PEV> numOuterSweepsValidator =
243 rcp_implicit_cast<PEV>(
rcp(
new NonnegativeIntValidator));
244 pl->set(
"relaxation: outer sweeps", numOuterSweeps,
"Number of outer local relaxation sweeps for two-stage GS",
245 rcp_const_cast<const PEV>(numOuterSweepsValidator));
247 const int numInnerSweeps = 1;
248 RCP<PEV> numInnerSweepsValidator =
249 rcp_implicit_cast<PEV>(
rcp(
new NonnegativeIntValidator));
250 pl->set(
"relaxation: inner sweeps", numInnerSweeps,
"Number of inner local relaxation sweeps for two-stage GS",
251 rcp_const_cast<const PEV>(numInnerSweepsValidator));
254 pl->set(
"relaxation: inner damping factor", innerDampingFactor,
"Damping factor for the inner sweep of two-stage GS");
256 const bool innerSpTrsv =
false;
257 pl->set(
"relaxation: inner sparse-triangular solve", innerSpTrsv,
"Specify whether to use sptrsv instead of JR iterations for two-stage GS");
259 const bool compactForm =
false;
260 pl->set(
"relaxation: compact form", compactForm,
"Specify whether to use compact form of recurrence for two-stage GS");
263 pl->set(
"relaxation: damping factor", dampingFactor);
265 const bool zeroStartingSolution =
true;
266 pl->set(
"relaxation: zero starting solution", zeroStartingSolution);
268 const bool doBackwardGS =
false;
269 pl->set(
"relaxation: backward mode", doBackwardGS);
271 const bool doL1Method =
false;
272 pl->set(
"relaxation: use l1", doL1Method);
274 const magnitude_type l1eta = (STM::one() + STM::one() + STM::one()) /
275 (STM::one() + STM::one());
276 pl->set(
"relaxation: l1 eta", l1eta);
279 pl->set(
"relaxation: min diagonal value", minDiagonalValue);
281 const bool fixTinyDiagEntries =
false;
282 pl->set(
"relaxation: fix tiny diagonal entries", fixTinyDiagEntries);
284 const bool checkDiagEntries =
false;
285 pl->set(
"relaxation: check diagonal entries", checkDiagEntries);
288 pl->set(
"relaxation: local smoothing indices", localSmoothingIndices);
290 const bool is_matrix_structurally_symmetric =
false;
291 pl->set(
"relaxation: symmetric matrix structure", is_matrix_structurally_symmetric);
293 const bool ifpack2_dump_matrix =
false;
294 pl->set(
"relaxation: ifpack2 dump matrix", ifpack2_dump_matrix);
296 const int cluster_size = 1;
297 pl->set(
"relaxation: mtgs cluster size", cluster_size);
299 pl->set(
"relaxation: mtgs coloring algorithm",
"Default");
301 const int long_row_threshold = 0;
302 pl->set(
"relaxation: long row threshold", long_row_threshold);
304 const bool timer_for_apply =
true;
305 pl->set(
"timer for apply", timer_for_apply);
307 validParams_ = rcp_const_cast<
const ParameterList>(pl);
312 template <
class MatrixType>
314 using Teuchos::getIntegralValue;
315 using Teuchos::getStringValue;
318 typedef scalar_type ST;
320 if (pl.
isType<
double>(
"relaxation: damping factor")) {
323 ST df = pl.
get<
double>(
"relaxation: damping factor");
324 pl.
remove(
"relaxation: damping factor");
325 pl.
set(
"relaxation: damping factor", df);
330 const Details::RelaxationType precType =
331 getIntegralValue<Details::RelaxationType>(pl,
"relaxation: type");
332 const std::string precTypeStr = getStringValue<Details::RelaxationType>(pl,
"relaxation: type");
334 pl.set<std::string>(
"relaxation: type", precTypeStr);
335 pl.get<std::string>(
"relaxation: type");
336 const int numSweeps = pl.get<
int>(
"relaxation: sweeps");
337 const ST dampingFactor = pl.get<ST>(
"relaxation: damping factor");
338 const bool zeroStartSol = pl.get<
bool>(
"relaxation: zero starting solution");
339 const bool doBackwardGS = pl.get<
bool>(
"relaxation: backward mode");
340 const bool doL1Method = pl.get<
bool>(
"relaxation: use l1");
341 const magnitude_type l1Eta = pl.get<magnitude_type>(
"relaxation: l1 eta");
342 const ST minDiagonalValue = pl.get<ST>(
"relaxation: min diagonal value");
343 const bool fixTinyDiagEntries = pl.get<
bool>(
"relaxation: fix tiny diagonal entries");
344 const bool checkDiagEntries = pl.get<
bool>(
"relaxation: check diagonal entries");
345 const bool is_matrix_structurally_symmetric = pl.get<
bool>(
"relaxation: symmetric matrix structure");
346 const bool ifpack2_dump_matrix = pl.get<
bool>(
"relaxation: ifpack2 dump matrix");
347 const bool timer_for_apply = pl.get<
bool>(
"timer for apply");
348 int cluster_size = 1;
349 if (pl.isParameter(
"relaxation: mtgs cluster size"))
350 cluster_size = pl.get<
int>(
"relaxation: mtgs cluster size");
351 int long_row_threshold = 0;
352 if (pl.isParameter(
"relaxation: long row threshold"))
353 long_row_threshold = pl.get<
int>(
"relaxation: long row threshold");
354 std::string color_algo_name = pl.get<std::string>(
"relaxation: mtgs coloring algorithm");
356 for (
char& c : color_algo_name)
358 if (color_algo_name ==
"default")
359 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_DEFAULT;
360 else if (color_algo_name ==
"serial")
361 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_SERIAL;
362 else if (color_algo_name ==
"vb")
363 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VB;
364 else if (color_algo_name ==
"vbbit")
365 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBBIT;
366 else if (color_algo_name ==
"vbcs")
367 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBCS;
368 else if (color_algo_name ==
"vbd")
369 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBD;
370 else if (color_algo_name ==
"vbdbit")
371 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBDBIT;
372 else if (color_algo_name ==
"eb")
373 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_EB;
375 std::ostringstream msg;
376 msg <<
"Ifpack2::Relaxation: 'relaxation: mtgs coloring algorithm' = '" << color_algo_name <<
"' is not valid.\n";
377 msg <<
"Choices (not case sensitive) are: Default, Serial, VB, VBBIT, VBCS, VBD, VBDBIT, EB.";
379 true, std::invalid_argument, msg.str());
385 if (!std::is_same<double, ST>::value && pl.isType<
double>(
"relaxation: inner damping factor")) {
388 ST df = pl.
get<
double>(
"relaxation: inner damping factor");
389 pl.remove(
"relaxation: inner damping factor");
390 pl.set(
"relaxation: inner damping factor", df);
393 if (long_row_threshold > 0) {
395 cluster_size != 1, std::invalid_argument,
396 "Ifpack2::Relaxation: "
397 "Requested long row MTGS/MTSGS algorithm and cluster GS/SGS, but those are not compatible.");
399 precType != Details::RelaxationType::MTGS && precType != Details::RelaxationType::MTSGS,
400 std::invalid_argument,
401 "Ifpack2::Relaxation: "
402 "Requested long row MTGS/MTSGS algorithm, but this is only compatible with preconditioner types "
403 "'MT Gauss-Seidel' and 'MT Symmetric Gauss-Seidel'.");
406 const ST innerDampingFactor = pl.get<ST>(
"relaxation: inner damping factor");
407 const int numInnerSweeps = pl.get<
int>(
"relaxation: inner sweeps");
408 const int numOuterSweeps = pl.get<
int>(
"relaxation: outer sweeps");
409 const bool innerSpTrsv = pl.get<
bool>(
"relaxation: inner sparse-triangular solve");
410 const bool compactForm = pl.get<
bool>(
"relaxation: compact form");
413 PrecType_ = precType;
414 NumSweeps_ = numSweeps;
415 DampingFactor_ = dampingFactor;
416 ZeroStartingSolution_ = zeroStartSol;
417 DoBackwardGS_ = doBackwardGS;
418 DoL1Method_ = doL1Method;
420 MinDiagonalValue_ = minDiagonalValue;
421 fixTinyDiagEntries_ = fixTinyDiagEntries;
422 checkDiagEntries_ = checkDiagEntries;
423 clusterSize_ = cluster_size;
424 longRowThreshold_ = long_row_threshold;
425 is_matrix_structurally_symmetric_ = is_matrix_structurally_symmetric;
426 ifpack2_dump_matrix_ = ifpack2_dump_matrix;
427 localSmoothingIndices_ = localSmoothingIndices;
429 NumInnerSweeps_ = numInnerSweeps;
430 NumOuterSweeps_ = numOuterSweeps;
431 InnerSpTrsv_ = innerSpTrsv;
432 InnerDampingFactor_ = innerDampingFactor;
433 CompactForm_ = compactForm;
434 TimerForApply_ = timer_for_apply;
437 template <
class MatrixType>
441 this->setParametersImpl(const_cast<Teuchos::ParameterList&>(pl));
444 template <
class MatrixType>
448 A_.
is_null(), std::runtime_error,
449 "Ifpack2::Relaxation::getComm: "
450 "The input matrix A is null. Please call setMatrix() with a nonnull "
451 "input matrix before calling this method.");
452 return A_->getRowMap()->getComm();
455 template <
class MatrixType>
461 template <
class MatrixType>
462 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
463 typename MatrixType::global_ordinal_type,
464 typename MatrixType::node_type>>
467 A_.
is_null(), std::runtime_error,
468 "Ifpack2::Relaxation::getDomainMap: "
469 "The input matrix A is null. Please call setMatrix() with a nonnull "
470 "input matrix before calling this method.");
471 return A_->getDomainMap();
474 template <
class MatrixType>
475 Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
476 typename MatrixType::global_ordinal_type,
477 typename MatrixType::node_type>>
480 A_.
is_null(), std::runtime_error,
481 "Ifpack2::Relaxation::getRangeMap: "
482 "The input matrix A is null. Please call setMatrix() with a nonnull "
483 "input matrix before calling this method.");
484 return A_->getRangeMap();
487 template <
class MatrixType>
492 template <
class MatrixType>
494 return (NumInitialize_);
497 template <
class MatrixType>
499 return (NumCompute_);
502 template <
class MatrixType>
507 template <
class MatrixType>
509 return (InitializeTime_);
512 template <
class MatrixType>
514 return (ComputeTime_);
517 template <
class MatrixType>
522 template <
class MatrixType>
524 return (ComputeFlops_);
527 template <
class MatrixType>
529 return (ApplyFlops_);
532 template <
class MatrixType>
535 A_.
is_null(), std::runtime_error,
536 "Ifpack2::Relaxation::getNodeSmootherComplexity: "
537 "The input matrix A is null. Please call setMatrix() with a nonnull "
538 "input matrix, then call compute(), before calling this method.");
540 return A_->getLocalNumRows() + A_->getLocalNumEntries();
543 template <
class MatrixType>
545 apply(
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
546 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
553 using Teuchos::rcpFromRef;
558 A_.
is_null(), std::runtime_error,
559 "Ifpack2::Relaxation::apply: "
560 "The input matrix A is null. Please call setMatrix() with a nonnull "
561 "input matrix, then call compute(), before calling this method.");
565 "Ifpack2::Relaxation::apply: You must call compute() on this Ifpack2 "
566 "preconditioner instance before you may call apply(). You may call "
567 "isComputed() to find out if compute() has been called already.");
568 TEUCHOS_TEST_FOR_EXCEPTION(
569 X.getNumVectors() != Y.getNumVectors(),
571 "Ifpack2::Relaxation::apply: X and Y have different numbers of columns. "
573 << X.getNumVectors() <<
" columns, but Y has "
574 << Y.getNumVectors() <<
" columns.");
575 TEUCHOS_TEST_FOR_EXCEPTION(
576 beta != STS::zero(), std::logic_error,
577 "Ifpack2::Relaxation::apply: beta = " << beta <<
" != 0 case not "
581 const std::string timerName(
"Ifpack2::Relaxation::apply");
582 if (TimerForApply_) {
598 if (alpha == STS::zero()) {
600 Y.putScalar(STS::zero());
609 Xcopy = rcpFromRef(X);
617 case Ifpack2::Details::JACOBI:
618 ApplyInverseJacobi(*Xcopy, Y);
620 case Ifpack2::Details::GS:
621 ApplyInverseSerialGS(*Xcopy, Y, DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward);
623 case Ifpack2::Details::SGS:
624 ApplyInverseSerialGS(*Xcopy, Y, Tpetra::Symmetric);
626 case Ifpack2::Details::MTGS:
627 case Ifpack2::Details::GS2:
628 ApplyInverseMTGS_CrsMatrix(*Xcopy, Y, DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward);
630 case Ifpack2::Details::MTSGS:
631 case Ifpack2::Details::SGS2:
632 ApplyInverseMTGS_CrsMatrix(*Xcopy, Y, Tpetra::Symmetric);
634 case Ifpack2::Details::RICHARDSON:
635 ApplyInverseRichardson(*Xcopy, Y);
639 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
640 "Ifpack2::Relaxation::apply: Invalid preconditioner type enum value "
641 << PrecType_ <<
". Please report this bug to the Ifpack2 developers.");
643 if (alpha != STS::one()) {
645 const double numGlobalRows = as<double>(A_->getGlobalNumRows());
646 const double numVectors = as<double>(Y.getNumVectors());
647 ApplyFlops_ += numGlobalRows * numVectors;
651 ApplyTime_ += (time.
wallTime() - startTime);
655 template <
class MatrixType>
657 applyMat(
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
658 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
661 A_.
is_null(), std::runtime_error,
662 "Ifpack2::Relaxation::applyMat: "
663 "The input matrix A is null. Please call setMatrix() with a nonnull "
664 "input matrix, then call compute(), before calling this method.");
666 !isComputed(), std::runtime_error,
667 "Ifpack2::Relaxation::applyMat: "
668 "isComputed() must be true before you may call applyMat(). "
669 "Please call compute() before calling this method.");
670 TEUCHOS_TEST_FOR_EXCEPTION(
671 X.getNumVectors() != Y.getNumVectors(), std::invalid_argument,
672 "Ifpack2::Relaxation::applyMat: X.getNumVectors() = " << X.getNumVectors()
673 <<
" != Y.getNumVectors() = " << Y.getNumVectors() <<
".");
674 A_->apply(X, Y, mode);
677 template <
class MatrixType>
679 const char methodName[] =
"Ifpack2::Relaxation::initialize";
682 "input matrix A is null. Please call setMatrix() with "
683 "a nonnull input matrix before calling this method.");
688 double startTime = timer->wallTime();
692 isInitialized_ =
false;
695 auto rowMap = A_->getRowMap();
696 auto comm = rowMap.
is_null() ? Teuchos::null : rowMap->getComm();
697 IsParallel_ = !comm.is_null() && comm->getSize() > 1;
707 Importer_ = IsParallel_ ? A_->getGraph()->getImporter() : Teuchos::null;
711 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(A_);
712 hasBlockCrsMatrix_ = !A_bcrs.
is_null();
715 serialGaussSeidel_ = Teuchos::null;
717 if (PrecType_ == Details::MTGS || PrecType_ == Details::MTSGS ||
718 PrecType_ == Details::GS2 || PrecType_ == Details::SGS2) {
719 auto crsMat = Details::getCrsMatrix(A_);
721 "Multithreaded Gauss-Seidel methods currently only work "
722 "when the input matrix is a Tpetra::CrsMatrix.");
727 if (ifpack2_dump_matrix_) {
728 static int sequence_number = 0;
729 const std::string file_name =
"Ifpack2_MT_GS_" +
730 std::to_string(sequence_number++) +
".mtx";
731 using writer_type = Tpetra::MatrixMarket::Writer<crs_matrix_type>;
732 writer_type::writeSparseFile(file_name, crsMat);
735 this->mtKernelHandle_ =
Teuchos::rcp(
new mt_kernel_handle_type());
736 if (mtKernelHandle_->get_gs_handle() ==
nullptr) {
737 if (PrecType_ == Details::GS2 || PrecType_ == Details::SGS2)
738 mtKernelHandle_->create_gs_handle(KokkosSparse::GS_TWOSTAGE);
739 else if (this->clusterSize_ == 1) {
740 mtKernelHandle_->create_gs_handle(KokkosSparse::GS_DEFAULT, this->mtColoringAlgorithm_);
741 mtKernelHandle_->get_point_gs_handle()->set_long_row_threshold(longRowThreshold_);
743 mtKernelHandle_->create_gs_handle(KokkosSparse::CLUSTER_DEFAULT, this->clusterSize_, this->mtColoringAlgorithm_);
745 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice();
746 if (PrecType_ == Details::GS2 || PrecType_ == Details::SGS2) {
748 mtKernelHandle_->set_gs_set_num_inner_sweeps(NumInnerSweeps_);
749 mtKernelHandle_->set_gs_set_num_outer_sweeps(NumOuterSweeps_);
750 mtKernelHandle_->set_gs_set_inner_damp_factor(InnerDampingFactor_);
751 mtKernelHandle_->set_gs_twostage(!InnerSpTrsv_, A_->getLocalNumRows());
752 mtKernelHandle_->set_gs_twostage_compact_form(CompactForm_);
755 KokkosSparse::gauss_seidel_symbolic(
756 mtKernelHandle_.getRawPtr(),
757 A_->getLocalNumRows(),
758 A_->getLocalNumCols(),
761 is_matrix_structurally_symmetric_);
765 InitializeTime_ += (timer->wallTime() - startTime);
767 isInitialized_ =
true;
771 template <
typename BlockDiagView>
772 struct InvertDiagBlocks {
773 typedef typename BlockDiagView::size_type Size;
776 typedef Kokkos::MemoryTraits<Kokkos::Unmanaged> Unmanaged;
777 template <
typename View>
778 using UnmanagedView = Kokkos::View<
typename View::data_type,
typename View::array_layout,
779 typename View::device_type, Unmanaged>;
781 typedef typename BlockDiagView::non_const_value_type Scalar;
782 typedef typename BlockDiagView::device_type Device;
783 typedef Kokkos::View<Scalar**, Kokkos::LayoutRight, Device> RWrk;
784 typedef Kokkos::View<int**, Kokkos::LayoutRight, Device> IWrk;
786 UnmanagedView<BlockDiagView> block_diag_;
790 UnmanagedView<RWrk> rwrk_;
792 UnmanagedView<IWrk> iwrk_;
795 InvertDiagBlocks(BlockDiagView& block_diag)
796 : block_diag_(block_diag) {
797 const auto blksz = block_diag.extent(1);
798 Kokkos::resize(rwrk_buf_, block_diag_.extent(0), blksz);
800 Kokkos::resize(iwrk_buf_, block_diag_.extent(0), blksz);
804 KOKKOS_INLINE_FUNCTION
805 void operator()(
const Size i,
int& jinfo)
const {
806 auto D_cur = Kokkos::subview(block_diag_, i, Kokkos::ALL(), Kokkos::ALL());
807 auto ipiv = Kokkos::subview(iwrk_, i, Kokkos::ALL());
808 auto work = Kokkos::subview(rwrk_, i, Kokkos::ALL());
810 Tpetra::GETF2(D_cur, ipiv, info);
815 Tpetra::GETRI(D_cur, ipiv, work, info);
821 template <
class MatrixType>
822 void Relaxation<MatrixType>::computeBlockCrs() {
831 using Teuchos::rcp_dynamic_cast;
835 using Teuchos::reduceAll;
836 typedef local_ordinal_type LO;
838 const std::string timerName(
"Ifpack2::Relaxation::computeBlockCrs");
843 double startTime = timer->
wallTime();
848 A_.
is_null(), std::runtime_error,
849 "Ifpack2::Relaxation::"
850 "computeBlockCrs: The input matrix A is null. Please call setMatrix() "
851 "with a nonnull input matrix, then call initialize(), before calling "
853 auto blockCrsA = rcp_dynamic_cast<
const block_crs_matrix_type>(A_);
855 blockCrsA.is_null(), std::logic_error,
856 "Ifpack2::Relaxation::"
857 "computeBlockCrs: A_ is not a BlockCrsMatrix, but it should be if we "
858 "got this far. Please report this bug to the Ifpack2 developers.");
860 const scalar_type one = STS::one();
865 const LO lclNumMeshRows =
866 blockCrsA->getCrsGraph().getLocalNumRows();
867 const LO blockSize = blockCrsA->getBlockSize();
869 if (!savedDiagOffsets_) {
870 blockDiag_ = block_diag_type();
871 blockDiag_ = block_diag_type(
"Ifpack2::Relaxation::blockDiag_",
872 lclNumMeshRows, blockSize, blockSize);
873 if (Teuchos::as<LO>(diagOffsets_.extent(0)) < lclNumMeshRows) {
876 diagOffsets_ = Kokkos::View<size_t*, device_type>();
877 diagOffsets_ = Kokkos::View<size_t*, device_type>(
"offsets", lclNumMeshRows);
879 blockCrsA->getCrsGraph().getLocalDiagOffsets(diagOffsets_);
881 static_cast<size_t>(blockDiag_.extent(0)),
882 std::logic_error,
"diagOffsets_.extent(0) = " << diagOffsets_.extent(0) <<
" != blockDiag_.extent(0) = " << blockDiag_.extent(0) <<
". Please report this bug to the Ifpack2 developers.");
883 savedDiagOffsets_ =
true;
885 blockCrsA->getLocalDiagCopy(blockDiag_, diagOffsets_);
892 unmanaged_block_diag_type blockDiag = blockDiag_;
894 if (DoL1Method_ && IsParallel_) {
895 const scalar_type two = one + one;
896 const size_t maxLength = A_->getLocalMaxNumRowEntries();
897 nonconst_local_inds_host_view_type indices(
"indices", maxLength);
898 nonconst_values_host_view_type values_(
"values", maxLength * blockSize * blockSize);
899 size_t numEntries = 0;
901 for (LO i = 0; i < lclNumMeshRows; ++i) {
903 blockCrsA->getLocalRowCopy(i, indices, values_, numEntries);
904 scalar_type* values =
reinterpret_cast<scalar_type*
>(values_.data());
906 auto diagBlock = Kokkos::subview(blockDiag, i, ALL(), ALL());
907 for (LO subRow = 0; subRow < blockSize; ++subRow) {
908 magnitude_type diagonal_boost = STM::zero();
909 for (
size_t k = 0; k < numEntries; ++k) {
910 if (indices[k] >= lclNumMeshRows) {
911 const size_t offset = blockSize * blockSize * k + subRow * blockSize;
912 for (LO subCol = 0; subCol < blockSize; ++subCol) {
913 diagonal_boost += STS::magnitude(values[offset + subCol] / two);
917 if (STS::magnitude(diagBlock(subRow, subRow)) < L1Eta_ * diagonal_boost) {
918 diagBlock(subRow, subRow) += diagonal_boost;
926 Impl::InvertDiagBlocks<unmanaged_block_diag_type> idb(blockDiag);
927 typedef typename unmanaged_block_diag_type::execution_space exec_space;
928 typedef Kokkos::RangePolicy<exec_space, LO> range_type;
930 Kokkos::parallel_reduce(range_type(0, lclNumMeshRows), idb, info);
940 #ifdef HAVE_IFPACK2_DEBUG
941 const int numResults = 2;
943 int lclResults[2], gblResults[2];
944 lclResults[0] = info;
945 lclResults[1] = -info;
948 reduceAll<int, int>(*(A_->getGraph()->getComm()), REDUCE_MIN,
949 numResults, lclResults, gblResults);
951 "Ifpack2::Relaxation::compute: When processing the input "
952 "Tpetra::BlockCrsMatrix, one or more diagonal block LU factorizations "
953 "failed on one or more (MPI) processes.");
954 #endif // HAVE_IFPACK2_DEBUG
955 serialGaussSeidel_ =
rcp(
new SerialGaussSeidel(blockCrsA, blockDiag_, localSmoothingIndices_, DampingFactor_));
958 ComputeTime_ += (timer->
wallTime() - startTime);
963 template <
class MatrixType>
975 using Teuchos::reduceAll;
979 using IST =
typename vector_type::impl_scalar_type;
980 #if KOKKOS_VERSION >= 40799
981 using KAT = KokkosKernels::ArithTraits<IST>;
983 using KAT = Kokkos::ArithTraits<IST>;
986 const char methodName[] =
"Ifpack2::Relaxation::compute";
987 const scalar_type zero = STS::zero();
988 const scalar_type one = STS::one();
991 "The input matrix A is null. Please call setMatrix() with a nonnull "
992 "input matrix, then call initialize(), before calling this method.");
995 if (!isInitialized()) {
999 if (hasBlockCrsMatrix_) {
1006 double startTime = timer->
wallTime();
1018 IST oneOverMinDiagVal = KAT::one() /
static_cast<IST
>(SmallTraits<scalar_type>::eps());
1019 if (MinDiagonalValue_ != zero)
1020 oneOverMinDiagVal = KAT::one() /
static_cast<IST
>(MinDiagonalValue_);
1023 const magnitude_type minDiagValMag = STS::magnitude(MinDiagonalValue_);
1025 const LO numMyRows =
static_cast<LO
>(A_->getLocalNumRows());
1028 "Please report this bug to the Ifpack2 developers.");
1029 IsComputed_ =
false;
1031 if (Diagonal_.is_null()) {
1034 Diagonal_ =
rcp(
new vector_type(A_->getRowMap(),
false));
1037 if (checkDiagEntries_) {
1043 size_t numSmallDiagEntries = 0;
1044 size_t numZeroDiagEntries = 0;
1045 size_t numNegDiagEntries = 0;
1052 A_->getLocalDiagCopy(*Diagonal_);
1053 std::unique_ptr<vector_type> origDiag;
1054 origDiag = std::unique_ptr<vector_type>(
new vector_type(*Diagonal_,
Teuchos::Copy));
1055 auto diag2d = Diagonal_->getLocalViewHost(Tpetra::Access::ReadWrite);
1056 auto diag = Kokkos::subview(diag2d, Kokkos::ALL(), 0);
1063 if (numMyRows != 0) {
1065 minMagDiagEntryMag = d_0_mag;
1066 maxMagDiagEntryMag = d_0_mag;
1075 for (LO i = 0; i < numMyRows; ++i) {
1076 const IST d_i = diag(i);
1080 const auto d_i_real = getRealValue(d_i);
1084 if (d_i_real < STM::zero()) {
1085 ++numNegDiagEntries;
1087 if (d_i_mag < minMagDiagEntryMag) {
1088 minMagDiagEntryMag = d_i_mag;
1090 if (d_i_mag > maxMagDiagEntryMag) {
1091 maxMagDiagEntryMag = d_i_mag;
1094 if (fixTinyDiagEntries_) {
1096 if (d_i_mag <= minDiagValMag) {
1097 ++numSmallDiagEntries;
1098 if (d_i_mag == STM::zero()) {
1099 ++numZeroDiagEntries;
1101 diag(i) = oneOverMinDiagVal;
1103 diag(i) = KAT::one() / d_i;
1107 if (d_i_mag <= minDiagValMag) {
1108 ++numSmallDiagEntries;
1109 if (d_i_mag == STM::zero()) {
1110 ++numZeroDiagEntries;
1113 diag(i) = KAT::one() / d_i;
1131 const Comm<int>& comm = *(A_->getRowMap()->getComm());
1134 Array<magnitude_type> localVals(2);
1135 localVals[0] = minMagDiagEntryMag;
1137 localVals[1] = -maxMagDiagEntryMag;
1138 Array<magnitude_type> globalVals(2);
1139 reduceAll<int, magnitude_type>(comm, REDUCE_MIN, 2,
1140 localVals.getRawPtr(),
1141 globalVals.getRawPtr());
1142 globalMinMagDiagEntryMag_ = globalVals[0];
1143 globalMaxMagDiagEntryMag_ = -globalVals[1];
1145 Array<size_t> localCounts(3);
1146 localCounts[0] = numSmallDiagEntries;
1147 localCounts[1] = numZeroDiagEntries;
1148 localCounts[2] = numNegDiagEntries;
1149 Array<size_t> globalCounts(3);
1150 reduceAll<int, size_t>(comm, REDUCE_SUM, 3,
1151 localCounts.getRawPtr(),
1152 globalCounts.getRawPtr());
1153 globalNumSmallDiagEntries_ = globalCounts[0];
1154 globalNumZeroDiagEntries_ = globalCounts[1];
1155 globalNumNegDiagEntries_ = globalCounts[2];
1159 vector_type diff(A_->getRowMap());
1160 diff.reciprocal(*origDiag);
1161 diff.update(-one, *Diagonal_, one);
1162 globalDiagNormDiff_ = diff.norm2();
1169 bool debugAgainstSlowPath =
false;
1171 auto crsMat = Details::getCrsMatrix(A_);
1173 if (crsMat.get() && crsMat->isFillComplete()) {
1178 if (invDiagKernel_.is_null())
1181 invDiagKernel_->setMatrix(crsMat);
1182 invDiagKernel_->compute(*Diagonal_,
1183 DoL1Method_ && IsParallel_, L1Eta_,
1184 fixTinyDiagEntries_, minDiagValMag);
1185 savedDiagOffsets_ =
true;
1189 #ifdef HAVE_IFPACK2_DEBUG
1190 debugAgainstSlowPath =
true;
1194 if (crsMat.is_null() || !crsMat->isFillComplete() || debugAgainstSlowPath) {
1202 if (debugAgainstSlowPath)
1203 Diagonal =
rcp(
new vector_type(A_->getRowMap()));
1205 Diagonal = Diagonal_;
1207 A_->getLocalDiagCopy(*Diagonal);
1217 if (DoL1Method_ && IsParallel_) {
1219 auto diag = Diagonal->getLocalViewHost(Tpetra::Access::ReadWrite);
1221 const size_t maxLength = A_row.getLocalMaxNumRowEntries();
1222 nonconst_local_inds_host_view_type indices(
"indices", maxLength);
1223 nonconst_values_host_view_type values(
"values", maxLength);
1226 for (LO i = 0; i < numMyRows; ++i) {
1227 A_row.getLocalRowCopy(i, indices, values, numEntries);
1229 for (
size_t k = 0; k < numEntries; ++k) {
1230 if (indices[k] >= numMyRows) {
1231 diagonal_boost += STS::magnitude(values[k] / two);
1234 if (KAT::magnitude(diag(i, 0)) < L1Eta_ * diagonal_boost) {
1235 diag(i, 0) += diagonal_boost;
1243 if (fixTinyDiagEntries_) {
1247 auto localDiag = Diagonal->getLocalViewDevice(Tpetra::Access::ReadWrite);
1248 Kokkos::parallel_for(
1249 Kokkos::RangePolicy<MyExecSpace>(0, localDiag.extent(0)),
1250 KOKKOS_LAMBDA(local_ordinal_type i) {
1251 auto d_i = localDiag(i, 0);
1254 if (d_i_mag <= minDiagValMag) {
1255 d_i = oneOverMinDiagVal;
1259 d_i = IST(KAT::one() / d_i);
1261 localDiag(i, 0) = d_i;
1264 Diagonal->reciprocal(*Diagonal);
1267 if (debugAgainstSlowPath) {
1269 Diagonal->update(STS::one(), *Diagonal_, -STS::one());
1274 <<
"\"fast-path\" diagonal computation failed. "
1275 "\\|D1 - D2\\|_inf = "
1283 if (STS::isComplex) {
1284 ComputeFlops_ += 4.0 * numMyRows;
1286 ComputeFlops_ += numMyRows;
1289 if (PrecType_ == Ifpack2::Details::MTGS ||
1290 PrecType_ == Ifpack2::Details::MTSGS ||
1291 PrecType_ == Ifpack2::Details::GS2 ||
1292 PrecType_ == Ifpack2::Details::SGS2) {
1294 using scalar_view_t =
typename local_matrix_device_type::values_type;
1297 "Multithreaded Gauss-Seidel methods currently only work "
1298 "when the input matrix is a Tpetra::CrsMatrix.");
1299 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice();
1303 auto diagView_2d = Diagonal_->getLocalViewDevice(Tpetra::Access::ReadWrite);
1304 scalar_view_t diagView_1d = Kokkos::subview(diagView_2d, Kokkos::ALL(), 0);
1305 KokkosSparse::gauss_seidel_numeric(
1306 mtKernelHandle_.getRawPtr(),
1307 A_->getLocalNumRows(),
1308 A_->getLocalNumCols(),
1313 is_matrix_structurally_symmetric_);
1314 }
else if (PrecType_ == Ifpack2::Details::GS || PrecType_ == Ifpack2::Details::SGS) {
1316 serialGaussSeidel_ =
rcp(
new SerialGaussSeidel(crsMat, Diagonal_, localSmoothingIndices_, DampingFactor_));
1318 serialGaussSeidel_ =
rcp(
new SerialGaussSeidel(A_, Diagonal_, localSmoothingIndices_, DampingFactor_));
1322 ComputeTime_ += (timer->
wallTime() - startTime);
1327 template <
class MatrixType>
1329 ApplyInverseRichardson(
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
1330 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y)
const {
1332 const double numGlobalRows = as<double>(A_->getGlobalNumRows());
1333 const double numVectors = as<double>(X.getNumVectors());
1334 if (ZeroStartingSolution_) {
1338 Y.scale(DampingFactor_, X);
1344 double flopUpdate = 0.0;
1345 if (DampingFactor_ == STS::one()) {
1347 flopUpdate = numGlobalRows * numVectors;
1351 flopUpdate = numGlobalRows * numVectors;
1353 ApplyFlops_ += flopUpdate;
1354 if (NumSweeps_ == 1) {
1360 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1363 updateCachedMultiVector(Y.getMap(), as<size_t>(numVectors));
1365 for (
int j = startSweep; j < NumSweeps_; ++j) {
1367 Tpetra::Details::residual(*A_, Y, X, *cachedMV_);
1368 Y.update(DampingFactor_, *cachedMV_, STS::one());
1382 const double numGlobalNonzeros = as<double>(A_->getGlobalNumEntries());
1383 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1384 ApplyFlops_ += as<double>(NumSweeps_ - startSweep) * numVectors *
1385 (2.0 * numGlobalNonzeros + dampingFlops);
1388 template <
class MatrixType>
1389 void Relaxation<MatrixType>::
1390 ApplyInverseJacobi(
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
1391 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y)
const {
1393 if (hasBlockCrsMatrix_) {
1394 ApplyInverseJacobi_BlockCrsMatrix(X, Y);
1398 const double numGlobalRows = as<double>(A_->getGlobalNumRows());
1399 const double numVectors = as<double>(X.getNumVectors());
1400 if (ZeroStartingSolution_) {
1405 Y.elementWiseMultiply(DampingFactor_, *Diagonal_, X, STS::zero());
1411 double flopUpdate = 0.0;
1412 if (DampingFactor_ == STS::one()) {
1414 flopUpdate = numGlobalRows * numVectors;
1418 flopUpdate = 2.0 * numGlobalRows * numVectors;
1420 ApplyFlops_ += flopUpdate;
1421 if (NumSweeps_ == 1) {
1427 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1430 updateCachedMultiVector(Y.getMap(), as<size_t>(numVectors));
1432 for (
int j = startSweep; j < NumSweeps_; ++j) {
1434 Tpetra::Details::residual(*A_, Y, X, *cachedMV_);
1435 Y.elementWiseMultiply(DampingFactor_, *Diagonal_, *cachedMV_, STS::one());
1449 const double numGlobalNonzeros = as<double>(A_->getGlobalNumEntries());
1450 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1451 ApplyFlops_ += as<double>(NumSweeps_ - startSweep) * numVectors *
1452 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1455 template <
class MatrixType>
1456 void Relaxation<MatrixType>::
1457 ApplyInverseJacobi_BlockCrsMatrix(
const Tpetra::MultiVector<scalar_type,
1459 global_ordinal_type,
1461 Tpetra::MultiVector<scalar_type,
1463 global_ordinal_type,
1464 node_type>& Y)
const {
1465 using Tpetra::BlockMultiVector;
1466 using BMV = BlockMultiVector<scalar_type, local_ordinal_type,
1467 global_ordinal_type, node_type>;
1469 const block_crs_matrix_type* blockMatConst =
1470 dynamic_cast<const block_crs_matrix_type*
>(A_.
getRawPtr());
1472 "This method should "
1473 "never be called if the matrix A_ is not a BlockCrsMatrix. "
1474 "Please report this bug to the Ifpack2 developers.");
1479 block_crs_matrix_type* blockMat =
1480 const_cast<block_crs_matrix_type*
>(blockMatConst);
1482 auto meshRowMap = blockMat->getRowMap();
1483 auto meshColMap = blockMat->getColMap();
1484 const local_ordinal_type blockSize = blockMat->getBlockSize();
1486 BMV X_blk(X, *meshColMap, blockSize);
1487 BMV Y_blk(Y, *meshRowMap, blockSize);
1489 if (ZeroStartingSolution_) {
1493 Y_blk.blockWiseMultiply(DampingFactor_, blockDiag_, X_blk);
1494 if (NumSweeps_ == 1) {
1499 auto pointRowMap = Y.getMap();
1500 const size_t numVecs = X.getNumVectors();
1504 BMV A_times_Y(*meshRowMap, *pointRowMap, blockSize, numVecs);
1508 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1510 for (
int j = startSweep; j < NumSweeps_; ++j) {
1511 blockMat->applyBlock(Y_blk, A_times_Y);
1513 Y_blk.blockJacobiUpdate(DampingFactor_, blockDiag_,
1514 X_blk, A_times_Y, STS::one());
1518 template <
class MatrixType>
1519 void Relaxation<MatrixType>::
1520 ApplyInverseSerialGS(
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
1521 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
1522 Tpetra::ESweepDirection direction)
const {
1523 using this_type = Relaxation<MatrixType>;
1533 auto blockCrsMat = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(A_);
1534 auto crsMat = Details::getCrsMatrix(A_);
1535 if (blockCrsMat.get()) {
1536 const_cast<this_type&
>(*this).ApplyInverseSerialGS_BlockCrsMatrix(*blockCrsMat, X, Y, direction);
1537 }
else if (crsMat.get()) {
1538 ApplyInverseSerialGS_CrsMatrix(*crsMat, X, Y, direction);
1540 ApplyInverseSerialGS_RowMatrix(X, Y, direction);
1544 template <
class MatrixType>
1545 void Relaxation<MatrixType>::
1546 ApplyInverseSerialGS_RowMatrix(
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
1547 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
1548 Tpetra::ESweepDirection direction)
const {
1555 using Teuchos::rcpFromRef;
1556 typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> MV;
1561 if (ZeroStartingSolution_) {
1562 Y.putScalar(STS::zero());
1565 size_t NumVectors = X.getNumVectors();
1569 if (Importer_.is_null()) {
1570 updateCachedMultiVector(Y.getMap(), NumVectors);
1572 updateCachedMultiVector(Importer_->getTargetMap(), NumVectors);
1579 for (
int j = 0; j < NumSweeps_; ++j) {
1582 if (Importer_.is_null()) {
1588 Y2->doImport(Y, *Importer_, Tpetra::INSERT);
1591 serialGaussSeidel_->apply(*Y2, X, direction);
1595 Tpetra::deep_copy(Y, *Y2);
1600 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1601 const double numVectors = as<double>(X.getNumVectors());
1602 const double numGlobalRows = as<double>(A_->getGlobalNumRows());
1603 const double numGlobalNonzeros = as<double>(A_->getGlobalNumEntries());
1604 ApplyFlops_ += 2.0 * NumSweeps_ * numVectors *
1605 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1608 template <
class MatrixType>
1609 void Relaxation<MatrixType>::
1610 ApplyInverseSerialGS_CrsMatrix(
const crs_matrix_type& A,
1611 const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& B,
1612 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
1613 Tpetra::ESweepDirection direction)
const {
1614 using Teuchos::null;
1617 using Teuchos::rcp_const_cast;
1618 using Teuchos::rcpFromRef;
1619 typedef scalar_type Scalar;
1620 const char prefix[] =
"Ifpack2::Relaxation::SerialGS: ";
1624 !A.isFillComplete(), std::runtime_error,
1625 prefix <<
"The matrix is not fill complete.");
1627 NumSweeps_ < 0, std::invalid_argument,
1628 prefix <<
"The number of sweeps must be nonnegative, "
1629 "but you provided numSweeps = "
1630 << NumSweeps_ <<
" < 0.");
1632 if (NumSweeps_ == 0) {
1636 RCP<const import_type> importer = A.getGraph()->getImporter();
1638 RCP<const map_type> domainMap = A.getDomainMap();
1639 RCP<const map_type> rangeMap = A.getRangeMap();
1640 RCP<const map_type> rowMap = A.getGraph()->getRowMap();
1641 RCP<const map_type> colMap = A.getGraph()->getColMap();
1643 #ifdef HAVE_IFPACK2_DEBUG
1649 !X.getMap()->isSameAs(*domainMap), std::runtime_error,
1650 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1651 "multivector X be in the domain Map of the matrix.");
1653 !B.getMap()->isSameAs(*rangeMap), std::runtime_error,
1654 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1655 "B be in the range Map of the matrix.");
1657 !Diagonal_->getMap()->isSameAs(*rowMap), std::runtime_error,
1658 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1659 "D be in the row Map of the matrix.");
1661 !rowMap->isSameAs(*rangeMap), std::runtime_error,
1662 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the row Map and the "
1663 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
1665 !domainMap->isSameAs(*rangeMap), std::runtime_error,
1666 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the domain Map and "
1667 "the range Map of the matrix be the same.");
1679 RCP<multivector_type> X_colMap;
1680 RCP<multivector_type> X_domainMap;
1681 bool copyBackOutput =
false;
1682 if (importer.is_null()) {
1683 X_colMap = Teuchos::rcpFromRef(X);
1684 X_domainMap = Teuchos::rcpFromRef(X);
1690 if (ZeroStartingSolution_) {
1691 X_colMap->putScalar(ZERO);
1695 updateCachedMultiVector(colMap, X.getNumVectors());
1696 X_colMap = cachedMV_;
1697 X_domainMap = X_colMap->offsetViewNonConst(domainMap, 0);
1699 if (ZeroStartingSolution_) {
1701 X_colMap->putScalar(ZERO);
1711 X_colMap->doImport(X, *importer, Tpetra::INSERT);
1713 copyBackOutput =
true;
1716 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
1717 if (!importer.is_null() && sweep > 0) {
1720 X_colMap->doImport(*X_domainMap, *importer, Tpetra::INSERT);
1723 serialGaussSeidel_->apply(*X_colMap, B, direction);
1726 if (copyBackOutput) {
1728 deep_copy(X, *X_domainMap);
1729 }
catch (std::exception& e) {
1731 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) "
1732 "threw an exception: "
1748 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1749 const double numVectors = X.getNumVectors();
1750 const double numGlobalRows = A_->getGlobalNumRows();
1751 const double numGlobalNonzeros = A_->getGlobalNumEntries();
1752 ApplyFlops_ += NumSweeps_ * numVectors *
1753 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1756 template <
class MatrixType>
1757 void Relaxation<MatrixType>::
1758 ApplyInverseSerialGS_BlockCrsMatrix(
const block_crs_matrix_type& A,
1759 const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
1760 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
1761 Tpetra::ESweepDirection direction) {
1764 using Teuchos::rcpFromRef;
1765 using Tpetra::INSERT;
1774 block_multivector_type yBlock(Y, *(A.getGraph()->getDomainMap()), A.getBlockSize());
1775 const block_multivector_type xBlock(X, *(A.getColMap()), A.getBlockSize());
1777 bool performImport =
false;
1778 RCP<block_multivector_type> yBlockCol;
1779 if (Importer_.is_null()) {
1780 yBlockCol = rcpFromRef(yBlock);
1782 if (yBlockColumnPointMap_.is_null() ||
1783 yBlockColumnPointMap_->getNumVectors() != yBlock.getNumVectors() ||
1784 yBlockColumnPointMap_->getBlockSize() != yBlock.getBlockSize()) {
1785 yBlockColumnPointMap_ =
1786 rcp(
new block_multivector_type(*(A.getColMap()), A.getBlockSize(),
1787 static_cast<local_ordinal_type
>(yBlock.getNumVectors())));
1789 yBlockCol = yBlockColumnPointMap_;
1790 if (pointImporter_.is_null()) {
1791 auto srcMap =
rcp(
new map_type(yBlock.getPointMap()));
1792 auto tgtMap =
rcp(
new map_type(yBlockCol->getPointMap()));
1793 pointImporter_ =
rcp(
new import_type(srcMap, tgtMap));
1795 performImport =
true;
1798 multivector_type yBlock_mv;
1799 multivector_type yBlockCol_mv;
1800 RCP<const multivector_type> yBlockColPointDomain;
1801 if (performImport) {
1802 yBlock_mv = yBlock.getMultiVectorView();
1803 yBlockCol_mv = yBlockCol->getMultiVectorView();
1804 yBlockColPointDomain =
1805 yBlockCol_mv.offsetView(A.getDomainMap(), 0);
1808 if (ZeroStartingSolution_) {
1809 yBlockCol->putScalar(STS::zero());
1810 }
else if (performImport) {
1811 yBlockCol_mv.doImport(yBlock_mv, *pointImporter_, Tpetra::INSERT);
1814 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
1815 if (performImport && sweep > 0) {
1816 yBlockCol_mv.doImport(yBlock_mv, *pointImporter_, Tpetra::INSERT);
1818 serialGaussSeidel_->applyBlock(*yBlockCol, xBlock, direction);
1819 if (performImport) {
1820 Tpetra::deep_copy(Y, *yBlockColPointDomain);
1825 template <
class MatrixType>
1826 void Relaxation<MatrixType>::
1827 ApplyInverseMTGS_CrsMatrix(
1828 const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& B,
1829 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
1830 Tpetra::ESweepDirection direction)
const {
1832 using Teuchos::null;
1835 using Teuchos::rcp_const_cast;
1836 using Teuchos::rcpFromRef;
1838 typedef scalar_type Scalar;
1840 const char prefix[] =
"Ifpack2::Relaxation::(reordered)MTGaussSeidel: ";
1843 auto crsMat = Details::getCrsMatrix(A_);
1845 "Ifpack2::Relaxation::apply: "
1846 "Multithreaded Gauss-Seidel methods currently only work when the "
1847 "input matrix is a Tpetra::CrsMatrix.");
1851 "Ifpack2's implementation of Multithreaded Gauss-Seidel does not "
1852 "implement the case where the user supplies an iteration order. "
1853 "This error used to appear as \"MT GaussSeidel ignores the given "
1855 "I tried to add more explanation, but I didn't implement \"MT "
1856 "GaussSeidel\" [sic]. "
1857 "You'll have to ask the person who did.");
1860 "input CrsMatrix is not fill complete. Please call fillComplete "
1861 "on the matrix, then do setup again, before calling apply(). "
1862 "\"Do setup\" means that if only the matrix's values have changed "
1863 "since last setup, you need only call compute(). If the matrix's "
1864 "structure may also have changed, you must first call initialize(), "
1865 "then call compute(). If you have not set up this preconditioner "
1866 "for this matrix before, you must first call initialize(), then "
1868 TEUCHOS_TEST_FOR_EXCEPTION(NumSweeps_ < 0, std::logic_error, prefix <<
": NumSweeps_ = " << NumSweeps_ <<
" < 0. Please report this bug to the Ifpack2 "
1870 if (NumSweeps_ == 0) {
1874 RCP<const import_type> importer = crsMat->getGraph()->getImporter();
1876 !crsMat->getGraph()->getExporter().is_null(), std::runtime_error,
1877 "This method's implementation currently requires that the matrix's row, "
1878 "domain, and range Maps be the same. This cannot be the case, because "
1879 "the matrix has a nontrivial Export object.");
1881 RCP<const map_type> domainMap = crsMat->getDomainMap();
1882 RCP<const map_type> rangeMap = crsMat->getRangeMap();
1883 RCP<const map_type> rowMap = crsMat->getGraph()->getRowMap();
1884 RCP<const map_type> colMap = crsMat->getGraph()->getColMap();
1886 #ifdef HAVE_IFPACK2_DEBUG
1892 !X.getMap()->isSameAs(*domainMap), std::runtime_error,
1893 "Ifpack2::Relaxation::MTGaussSeidel requires that the input "
1894 "multivector X be in the domain Map of the matrix.");
1896 !B.getMap()->isSameAs(*rangeMap), std::runtime_error,
1897 "Ifpack2::Relaxation::MTGaussSeidel requires that the input "
1898 "B be in the range Map of the matrix.");
1900 !rowMap->isSameAs(*rangeMap), std::runtime_error,
1901 "Ifpack2::Relaxation::MTGaussSeidel requires that the row Map and the "
1902 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
1904 !domainMap->isSameAs(*rangeMap), std::runtime_error,
1905 "Ifpack2::Relaxation::MTGaussSeidel requires that the domain Map and "
1906 "the range Map of the matrix be the same.");
1912 #endif // HAVE_IFPACK2_DEBUG
1922 RCP<multivector_type> X_colMap;
1923 RCP<multivector_type> X_domainMap;
1924 bool copyBackOutput =
false;
1925 if (importer.is_null()) {
1926 if (X.isConstantStride()) {
1927 X_colMap = rcpFromRef(X);
1928 X_domainMap = rcpFromRef(X);
1935 if (ZeroStartingSolution_) {
1936 X_colMap->putScalar(ZERO);
1944 updateCachedMultiVector(colMap, as<size_t>(X.getNumVectors()));
1945 X_colMap = cachedMV_;
1949 X_domainMap = X_colMap;
1950 if (!ZeroStartingSolution_) {
1952 deep_copy(*X_domainMap, X);
1953 }
catch (std::exception& e) {
1954 std::ostringstream os;
1955 os <<
"Ifpack2::Relaxation::MTGaussSeidel: "
1956 "deep_copy(*X_domainMap, X) threw an exception: "
1961 copyBackOutput =
true;
1974 updateCachedMultiVector(colMap, as<size_t>(X.getNumVectors()));
1975 X_colMap = cachedMV_;
1977 X_domainMap = X_colMap->offsetViewNonConst(domainMap, 0);
1979 if (ZeroStartingSolution_) {
1981 X_colMap->putScalar(ZERO);
1991 X_colMap->doImport(X, *importer, Tpetra::CombineMode::INSERT);
1993 copyBackOutput =
true;
1999 RCP<const multivector_type> B_in;
2000 if (B.isConstantStride()) {
2001 B_in = rcpFromRef(B);
2006 RCP<multivector_type> B_in_nonconst =
rcp(
new multivector_type(rowMap, B.getNumVectors()));
2008 deep_copy(*B_in_nonconst, B);
2009 }
catch (std::exception& e) {
2010 std::ostringstream os;
2011 os <<
"Ifpack2::Relaxation::MTGaussSeidel: "
2012 "deep_copy(*B_in_nonconst, B) threw an exception: "
2016 B_in = rcp_const_cast<
const multivector_type>(B_in_nonconst);
2030 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice();
2032 bool update_y_vector =
true;
2034 bool zero_x_vector =
false;
2036 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
2037 if (!importer.is_null() && sweep > 0) {
2040 X_colMap->doImport(*X_domainMap, *importer, Tpetra::CombineMode::INSERT);
2043 if (direction == Tpetra::Symmetric) {
2044 KokkosSparse::symmetric_gauss_seidel_apply(mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2045 kcsr.graph.row_map, kcsr.graph.entries, kcsr.values,
2046 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2047 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2048 zero_x_vector, update_y_vector, DampingFactor_, 1);
2049 }
else if (direction == Tpetra::Forward) {
2050 KokkosSparse::forward_sweep_gauss_seidel_apply(mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2051 kcsr.graph.row_map, kcsr.graph.entries, kcsr.values,
2052 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2053 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2054 zero_x_vector, update_y_vector, DampingFactor_, 1);
2055 }
else if (direction == Tpetra::Backward) {
2056 KokkosSparse::backward_sweep_gauss_seidel_apply(mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2057 kcsr.graph.row_map, kcsr.graph.entries, kcsr.values,
2058 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2059 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2060 zero_x_vector, update_y_vector, DampingFactor_, 1);
2063 true, std::invalid_argument,
2064 prefix <<
"The 'direction' enum does not have any of its valid "
2065 "values: Forward, Backward, or Symmetric.");
2067 update_y_vector =
false;
2070 if (copyBackOutput) {
2072 deep_copy(X, *X_domainMap);
2073 }
catch (std::exception& e) {
2075 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) "
2076 "threw an exception: "
2081 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
2082 const double numVectors = as<double>(X.getNumVectors());
2083 const double numGlobalRows = as<double>(A_->getGlobalNumRows());
2084 const double numGlobalNonzeros = as<double>(A_->getGlobalNumEntries());
2085 double ApplyFlops = NumSweeps_ * numVectors *
2086 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
2087 if (direction == Tpetra::Symmetric)
2089 ApplyFlops_ += ApplyFlops;
2092 template <
class MatrixType>
2094 std::ostringstream os;
2099 os <<
"\"Ifpack2::Relaxation\": {";
2101 os <<
"Initialized: " << (isInitialized() ?
"true" :
"false") <<
", "
2102 <<
"Computed: " << (isComputed() ?
"true" :
"false") <<
", ";
2108 if (PrecType_ == Ifpack2::Details::JACOBI) {
2110 }
else if (PrecType_ == Ifpack2::Details::GS) {
2111 os <<
"Gauss-Seidel";
2112 }
else if (PrecType_ == Ifpack2::Details::SGS) {
2113 os <<
"Symmetric Gauss-Seidel";
2114 }
else if (PrecType_ == Ifpack2::Details::MTGS) {
2115 os <<
"MT Gauss-Seidel";
2116 }
else if (PrecType_ == Ifpack2::Details::MTSGS) {
2117 os <<
"MT Symmetric Gauss-Seidel";
2118 }
else if (PrecType_ == Ifpack2::Details::GS2) {
2119 os <<
"Two-stage Gauss-Seidel";
2120 }
else if (PrecType_ == Ifpack2::Details::SGS2) {
2121 os <<
"Two-stage Symmetric Gauss-Seidel";
2125 if (hasBlockCrsMatrix_)
2129 <<
"sweeps: " << NumSweeps_ <<
", "
2130 <<
"damping factor: " << DampingFactor_ <<
", ";
2132 if (PrecType_ == Ifpack2::Details::MTGS || PrecType_ == Ifpack2::Details::MTSGS) {
2133 os <<
"\"relaxation: mtgs cluster size\": " << clusterSize_ <<
", ";
2134 os <<
"\"relaxation: long row threshold\": " << longRowThreshold_ <<
", ";
2135 os <<
"\"relaxation: symmetric matrix structure\": " << (is_matrix_structurally_symmetric_ ?
"true" :
"false") <<
", ";
2136 os <<
"\"relaxation: relaxation: mtgs coloring algorithm\": ";
2137 switch (mtColoringAlgorithm_) {
2138 case KokkosGraph::COLORING_DEFAULT:
2141 case KokkosGraph::COLORING_SERIAL:
2144 case KokkosGraph::COLORING_VB:
2147 case KokkosGraph::COLORING_VBBIT:
2150 case KokkosGraph::COLORING_VBCS:
2153 case KokkosGraph::COLORING_VBD:
2156 case KokkosGraph::COLORING_VBDBIT:
2159 case KokkosGraph::COLORING_EB:
2168 if (PrecType_ == Ifpack2::Details::GS2 ||
2169 PrecType_ == Ifpack2::Details::SGS2) {
2170 os <<
"outer sweeps: " << NumOuterSweeps_ <<
", "
2171 <<
"inner sweeps: " << NumInnerSweeps_ <<
", "
2172 <<
"inner damping factor: " << InnerDampingFactor_ <<
", ";
2176 os <<
"use l1: " << DoL1Method_ <<
", "
2177 <<
"l1 eta: " << L1Eta_ <<
", ";
2181 os <<
"Matrix: null";
2183 os <<
"Global matrix dimensions: ["
2184 << A_->getGlobalNumRows() <<
", " << A_->getGlobalNumCols() <<
"]"
2185 <<
", Global nnz: " << A_->getGlobalNumEntries();
2192 template <
class MatrixType>
2209 const int myRank = this->getComm()->getRank();
2221 out <<
"\"Ifpack2::Relaxation\":" << endl;
2223 out <<
"MatrixType: \"" << TypeNameTraits<MatrixType>::name() <<
"\""
2225 if (this->getObjectLabel() !=
"") {
2226 out <<
"Label: " << this->getObjectLabel() << endl;
2228 out <<
"Initialized: " << (isInitialized() ?
"true" :
"false") << endl
2229 <<
"Computed: " << (isComputed() ?
"true" :
"false") << endl
2230 <<
"Parameters: " << endl;
2233 out <<
"\"relaxation: type\": ";
2234 if (PrecType_ == Ifpack2::Details::JACOBI) {
2236 }
else if (PrecType_ == Ifpack2::Details::GS) {
2237 out <<
"Gauss-Seidel";
2238 }
else if (PrecType_ == Ifpack2::Details::SGS) {
2239 out <<
"Symmetric Gauss-Seidel";
2240 }
else if (PrecType_ == Ifpack2::Details::MTGS) {
2241 out <<
"MT Gauss-Seidel";
2242 }
else if (PrecType_ == Ifpack2::Details::MTSGS) {
2243 out <<
"MT Symmetric Gauss-Seidel";
2244 }
else if (PrecType_ == Ifpack2::Details::GS2) {
2245 out <<
"Two-stage Gauss-Seidel";
2246 }
else if (PrecType_ == Ifpack2::Details::SGS2) {
2247 out <<
"Two-stage Symmetric Gauss-Seidel";
2254 <<
"\"relaxation: sweeps\": " << NumSweeps_ << endl
2255 <<
"\"relaxation: damping factor\": " << DampingFactor_ << endl
2256 <<
"\"relaxation: min diagonal value\": " << MinDiagonalValue_ << endl
2257 <<
"\"relaxation: zero starting solution\": " << ZeroStartingSolution_ << endl
2258 <<
"\"relaxation: backward mode\": " << DoBackwardGS_ << endl
2259 <<
"\"relaxation: use l1\": " << DoL1Method_ << endl
2260 <<
"\"relaxation: l1 eta\": " << L1Eta_ << endl;
2261 if (PrecType_ == Ifpack2::Details::MTGS || PrecType_ == Ifpack2::Details::MTSGS) {
2262 out <<
"\"relaxation: mtgs cluster size\": " << clusterSize_ << endl;
2263 out <<
"\"relaxation: long row threshold\": " << longRowThreshold_ << endl;
2264 out <<
"\"relaxation: symmetric matrix structure\": " << (is_matrix_structurally_symmetric_ ?
"true" :
"false") << endl;
2265 out <<
"\"relaxation: relaxation: mtgs coloring algorithm\": ";
2266 switch (mtColoringAlgorithm_) {
2267 case KokkosGraph::COLORING_DEFAULT:
2270 case KokkosGraph::COLORING_SERIAL:
2273 case KokkosGraph::COLORING_VB:
2276 case KokkosGraph::COLORING_VBBIT:
2279 case KokkosGraph::COLORING_VBCS:
2282 case KokkosGraph::COLORING_VBD:
2285 case KokkosGraph::COLORING_VBDBIT:
2288 case KokkosGraph::COLORING_EB:
2296 if (PrecType_ == Ifpack2::Details::GS2 || PrecType_ == Ifpack2::Details::SGS2) {
2297 out <<
"\"relaxation: inner damping factor\": " << InnerDampingFactor_ << endl;
2298 out <<
"\"relaxation: outer sweeps\" : " << NumOuterSweeps_ << endl;
2299 out <<
"\"relaxation: inner sweeps\" : " << NumInnerSweeps_ << endl;
2302 out <<
"Computed quantities:" << endl;
2305 out <<
"Global number of rows: " << A_->getGlobalNumRows() << endl
2306 <<
"Global number of columns: " << A_->getGlobalNumCols() << endl;
2308 if (checkDiagEntries_ && isComputed()) {
2309 out <<
"Properties of input diagonal entries:" << endl;
2312 out <<
"Magnitude of minimum-magnitude diagonal entry: "
2313 << globalMinMagDiagEntryMag_ << endl
2314 <<
"Magnitude of maximum-magnitude diagonal entry: "
2315 << globalMaxMagDiagEntryMag_ << endl
2316 <<
"Number of diagonal entries with small magnitude: "
2317 << globalNumSmallDiagEntries_ << endl
2318 <<
"Number of zero diagonal entries: "
2319 << globalNumZeroDiagEntries_ << endl
2320 <<
"Number of diagonal entries with negative real part: "
2321 << globalNumNegDiagEntries_ << endl
2322 <<
"Abs 2-norm diff between computed and actual inverse "
2323 <<
"diagonal: " << globalDiagNormDiff_ << endl;
2327 out <<
"Saved diagonal offsets: "
2328 << (savedDiagOffsets_ ?
"true" :
"false") << endl;
2330 out <<
"Call counts and total times (in seconds): " << endl;
2333 out <<
"initialize: " << endl;
2336 out <<
"Call count: " << NumInitialize_ << endl;
2337 out <<
"Total time: " << InitializeTime_ << endl;
2339 out <<
"compute: " << endl;
2342 out <<
"Call count: " << NumCompute_ << endl;
2343 out <<
"Total time: " << ComputeTime_ << endl;
2345 out <<
"apply: " << endl;
2348 out <<
"Call count: " << NumApply_ << endl;
2349 out <<
"Total time: " << ApplyTime_ << endl;
2357 #define IFPACK2_RELAXATION_INSTANT(S, LO, GO, N) \
2358 template class Ifpack2::Relaxation<Tpetra::RowMatrix<S, LO, GO, N>>;
2360 #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:488
double getComputeFlops() const
Total number of floating-point operations over all calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:523
double getApplyFlops() const
Total number of floating-point operations over all calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:528
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:2194
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:186
T & get(const std::string &name, T def_value)
void compute()
Compute the preconditioner ("numeric setup");.
Definition: Ifpack2_Relaxation_def.hpp:964
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:533
virtual const std::string getXMLTypeName() const =0
Diagonal preconditioner.
Definition: Ifpack2_Diagonal_decl.hpp:38
static std::ostream & printLines(std::ostream &os, const std::string &linePrefix, const std::string &lines)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix and vectors are distributed.
Definition: Ifpack2_Relaxation_def.hpp:446
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Relaxation_def.hpp:166
int getNumApply() const
Total number of calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:503
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:513
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:230
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:465
void setParameters(const Teuchos::ParameterList ¶ms)
Set the relaxation / preconditioner parameters.
Definition: Ifpack2_Relaxation_def.hpp:438
File for utility functions.
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_Relaxation_def.hpp:2093
int getNumInitialize() const
Total number of calls to initialize().
Definition: Ifpack2_Relaxation_def.hpp:493
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:478
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:227
double getInitializeTime() const
Total time in seconds spent in all calls to initialize().
Definition: Ifpack2_Relaxation_def.hpp:508
any & getAny(bool activeQry=true)
Compute scaled damped residual for Chebyshev.
Definition: Ifpack2_Details_InverseDiagonalKernel_decl.hpp:45
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:657
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix to be preconditioned.
Definition: Ifpack2_Relaxation_def.hpp:457
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:224
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:221
std::string typeName() const
const std::type_info & type() const
int getNumCompute() const
Total number of calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:498
Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices.
Definition: Ifpack2_Relaxation_decl.hpp:206
double getApplyTime() const
Total time in seconds spent in all calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:518
void getValidParameters(Teuchos::ParameterList ¶ms)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:18
void initialize()
Initialize the preconditioner ("symbolic setup").
Definition: Ifpack2_Relaxation_def.hpp:678
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:241
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:545
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a list of all the parameters that this class accepts.
Definition: Ifpack2_Relaxation_def.hpp:195
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_Relaxation_decl.hpp:236
virtual void validate(ParameterEntry const &entry, std::string const ¶mName, std::string const &sublistName) const =0