43#ifndef IFPACK2_LOCALSPARSETRIANGULARSOLVER_DEF_HPP
44#define IFPACK2_LOCALSPARSETRIANGULARSOLVER_DEF_HPP
46#include "Tpetra_CrsMatrix.hpp"
47#include "Tpetra_Core.hpp"
48#include "Teuchos_StandardParameterEntryValidators.hpp"
49#include "Tpetra_Details_determineLocalTriangularStructure.hpp"
50#include "KokkosSparse_trsv.hpp"
52#ifdef HAVE_IFPACK2_SHYLU_NODEHTS
53# include "shylu_hts.hpp"
66 static void loadPLTypeOption (Teuchos::Array<std::string>& type_strs, Teuchos::Array<Enum>& type_enums) {
68 type_strs[0] =
"Internal";
70 type_strs[2] =
"KSPTRSV";
72 type_enums[0] = Internal;
74 type_enums[2] = KSPTRSV;
79template<
class MatrixType>
80class LocalSparseTriangularSolver<MatrixType>::HtsImpl {
82 typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type>
crs_matrix_type;
85#ifdef HAVE_IFPACK2_SHYLU_NODEHTS
86 Timpl_ = Teuchos::null;
87 levelset_block_size_ = 1;
93#ifdef HAVE_IFPACK2_SHYLU_NODEHTS
94 const char* block_size_s =
"trisolver: block size";
95 if (pl.isParameter(block_size_s)) {
96 TEUCHOS_TEST_FOR_EXCEPT_MSG( ! pl.isType<
int>(block_size_s),
97 "The parameter \"" << block_size_s <<
"\" must be of type int.");
98 levelset_block_size_ = pl.get<
int>(block_size_s);
100 if (levelset_block_size_ < 1)
101 levelset_block_size_ = 1;
110#ifdef HAVE_IFPACK2_SHYLU_NODEHTS
112 transpose_ = conjugate_ =
false;
119#ifdef HAVE_IFPACK2_SHYLU_NODEHTS
120 using Teuchos::ArrayRCP;
122 auto rowptr = T_in.getLocalRowPtrsHost();
123 auto colidx = T_in.getLocalIndicesHost();
124 auto val = T_in.getLocalValuesHost(Tpetra::Access::ReadOnly);
127 Teuchos::RCP<HtsCrsMatrix> T_hts = Teuchos::rcpWithDealloc(
128 HTST::make_CrsMatrix(rowptr.size() - 1,
129 rowptr.data(), colidx.data(),
132 transpose_, conjugate_),
133 HtsCrsMatrixDeleter());
135 if (Teuchos::nonnull(Timpl_)) {
137 HTST::reprocess_numeric(Timpl_.get(), T_hts.get());
140 if (T_in.getCrsGraph().is_null()) {
141 if (Teuchos::nonnull(out))
142 *out <<
"HTS compute failed because T_in.getCrsGraph().is_null().\n";
145 if ( ! T_in.getCrsGraph()->isSorted()) {
146 if (Teuchos::nonnull(out))
147 *out <<
"HTS compute failed because ! T_in.getCrsGraph().isSorted().\n";
150 if ( ! T_in.isStorageOptimized()) {
151 if (Teuchos::nonnull(out))
152 *out <<
"HTS compute failed because ! T_in.isStorageOptimized().\n";
156 typename HTST::PreprocessArgs args;
157 args.T = T_hts.get();
160 args.nthreads = omp_get_max_threads();
164 args.save_for_reprocess =
true;
165 typename HTST::Options opts;
166 opts.levelset_block_size = levelset_block_size_;
167 args.options = &opts;
170 Timpl_ = Teuchos::rcpWithDealloc(HTST::preprocess(args), TImplDeleter());
171 }
catch (
const std::exception& e) {
172 if (Teuchos::nonnull(out))
173 *out <<
"HTS preprocess threw: " << e.what() <<
"\n";
182#ifdef HAVE_IFPACK2_SHYLU_NODEHTS
183 return Teuchos::nonnull(Timpl_);
190 void localApply (
const MV& X, MV& Y,
191 const Teuchos::ETransp ,
197#ifdef HAVE_IFPACK2_SHYLU_NODEHTS
198 const auto& X_view = X.getLocalViewHost (Tpetra::Access::ReadOnly);
199 const auto& Y_view = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
202 HTST::reset_max_nrhs(Timpl_.get(), X_view.extent(1));
204 HTST::solve_omp(Timpl_.get(),
206 reinterpret_cast<const scalar_type*
>(X_view.data()),
215#ifdef HAVE_IFPACK2_SHYLU_NODEHTS
216 typedef ::Experimental::HTS<local_ordinal_type, size_t, scalar_type> HTST;
217 typedef typename HTST::Impl TImpl;
218 typedef typename HTST::CrsMatrix HtsCrsMatrix;
220 struct TImplDeleter {
221 void free (TImpl* impl) {
222 HTST::delete_Impl(impl);
226 struct HtsCrsMatrixDeleter {
227 void free (HtsCrsMatrix* T) {
228 HTST::delete_CrsMatrix(T);
232 Teuchos::RCP<TImpl> Timpl_;
233 bool transpose_, conjugate_;
234 int levelset_block_size_;
238template<
class MatrixType>
245 if (! A.is_null ()) {
246 Teuchos::RCP<const crs_matrix_type> A_crs =
247 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A);
248 TEUCHOS_TEST_FOR_EXCEPTION
249 (A_crs.is_null (), std::invalid_argument,
250 "Ifpack2::LocalSparseTriangularSolver constructor: "
251 "The input matrix A is not a Tpetra::CrsMatrix.");
256template<
class MatrixType>
259 const Teuchos::RCP<Teuchos::FancyOStream>& out) :
264 if (! out_.is_null ()) {
265 *out_ <<
">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor"
269 if (! A.is_null ()) {
270 Teuchos::RCP<const crs_matrix_type> A_crs =
271 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A);
272 TEUCHOS_TEST_FOR_EXCEPTION
273 (A_crs.is_null (), std::invalid_argument,
274 "Ifpack2::LocalSparseTriangularSolver constructor: "
275 "The input matrix A is not a Tpetra::CrsMatrix.");
280template<
class MatrixType>
287template<
class MatrixType>
293 if (! out_.is_null ()) {
294 *out_ <<
">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor"
299template<
class MatrixType>
302 isInitialized_ =
false;
304 reverseStorage_ =
false;
305 isInternallyChanged_ =
false;
309 initializeTime_ = 0.0;
312 isKokkosKernelsSptrsv_ =
false;
317template<
class MatrixType>
321 if (Teuchos::nonnull (kh_))
323 kh_->destroy_sptrsv_handle();
327template<
class MatrixType>
333 using Teuchos::ParameterList;
334 using Teuchos::Array;
336 Details::TrisolverType::Enum trisolverType = Details::TrisolverType::Internal;
338 static const char typeName[] =
"trisolver: type";
340 if ( ! pl.isType<std::string>(typeName))
break;
343 Array<std::string> trisolverTypeStrs;
344 Array<Details::TrisolverType::Enum> trisolverTypeEnums;
345 Details::TrisolverType::loadPLTypeOption (trisolverTypeStrs, trisolverTypeEnums);
346 Teuchos::StringToIntegralParameterEntryValidator<Details::TrisolverType::Enum>
347 s2i(trisolverTypeStrs (), trisolverTypeEnums (), typeName,
false);
349 trisolverType = s2i.getIntegralValue(pl.get<std::string>(typeName));
352 if (trisolverType == Details::TrisolverType::HTS) {
353 htsImpl_ = Teuchos::rcp (
new HtsImpl ());
354 htsImpl_->setParameters (pl);
357 if (trisolverType == Details::TrisolverType::KSPTRSV) {
358 kh_ = Teuchos::rcp (
new k_handle());
361 if (pl.isParameter(
"trisolver: reverse U"))
362 reverseStorage_ = pl.get<
bool>(
"trisolver: reverse U");
364 TEUCHOS_TEST_FOR_EXCEPTION
365 (reverseStorage_ && (trisolverType == Details::TrisolverType::HTS || trisolverType == Details::TrisolverType::KSPTRSV),
366 std::logic_error,
"Ifpack2::LocalSparseTriangularSolver::setParameters: "
367 "You are not allowed to enable both HTS or KSPTRSV and the \"trisolver: reverse U\" "
368 "options. See GitHub issue #2647.");
371template<
class MatrixType>
376 using Tpetra::Details::determineLocalTriangularStructure;
378 using local_matrix_type =
typename crs_matrix_type::local_matrix_device_type;
381 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::initialize: ";
382 if (! out_.is_null ()) {
383 *out_ <<
">>> DEBUG " << prefix << std::endl;
386 TEUCHOS_TEST_FOR_EXCEPTION
387 (A_.is_null (), std::runtime_error, prefix <<
"You must call "
388 "setMatrix() with a nonnull input matrix before you may call "
389 "initialize() or compute().");
390 if (A_crs_.is_null ()) {
391 auto A_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
392 TEUCHOS_TEST_FOR_EXCEPTION
393 (A_crs.get () ==
nullptr, std::invalid_argument,
394 prefix <<
"The input matrix A is not a Tpetra::CrsMatrix.");
397 auto G = A_crs_->getGraph ();
398 TEUCHOS_TEST_FOR_EXCEPTION
399 (G.is_null (), std::logic_error, prefix <<
"A_ and A_crs_ are nonnull, "
400 "but A_crs_'s RowGraph G is null. "
401 "Please report this bug to the Ifpack2 developers.");
405 TEUCHOS_TEST_FOR_EXCEPTION
406 (! G->isFillComplete (), std::runtime_error,
"If you call this method, "
407 "the matrix's graph must be fill complete. It is not.");
410 constexpr bool ignoreMapsForTriStructure =
true;
411 auto lclTriStructure = [&] {
412 auto lclMatrix = A_crs_->getLocalMatrixDevice ();
413 auto lclRowMap = A_crs_->getRowMap ()->getLocalMap ();
414 auto lclColMap = A_crs_->getColMap ()->getLocalMap ();
416 determineLocalTriangularStructure (lclMatrix.graph,
419 ignoreMapsForTriStructure);
420 const LO lclNumRows = lclRowMap.getLocalNumElements ();
421 this->diag_ = (lclTriStruct.diagCount < lclNumRows) ?
"U" :
"N";
422 this->uplo_ = lclTriStruct.couldBeLowerTriangular ?
"L" :
423 (lclTriStruct.couldBeUpperTriangular ?
"U" :
"N");
427 if (reverseStorage_ && lclTriStructure.couldBeUpperTriangular &&
428 htsImpl_.is_null ()) {
430 auto Alocal = A_crs_->getLocalMatrixDevice();
431 auto ptr = Alocal.graph.row_map;
432 auto ind = Alocal.graph.entries;
433 auto val = Alocal.values;
435 auto numRows = Alocal.numRows();
436 auto numCols = Alocal.numCols();
437 auto numNnz = Alocal.nnz();
439 typename decltype(ptr)::non_const_type newptr (
"ptr", ptr.extent (0));
440 typename decltype(ind)::non_const_type newind (
"ind", ind.extent (0));
441 decltype(val) newval (
"val", val.extent (0));
444 typename crs_matrix_type::execution_space().fence();
447 auto A_r = Alocal.row(numRows-1 - row);
449 auto numEnt = A_r.length;
451 newind(rowStart + k) = numCols-1 - A_r.colidx(numEnt-1 - k);
452 newval(rowStart + k) = A_r.value (numEnt-1 - k);
455 newptr(row+1) = rowStart;
457 typename crs_matrix_type::execution_space().fence();
460 Teuchos::RCP<map_type> newRowMap, newColMap;
463 auto rowMap = A_->getRowMap();
464 auto numElems = rowMap->getLocalNumElements();
465 auto rowElems = rowMap->getLocalElementList();
467 Teuchos::Array<global_ordinal_type> newRowElems(rowElems.size());
468 for (
size_t i = 0; i < numElems; i++)
469 newRowElems[i] = rowElems[numElems-1 - i];
471 newRowMap = Teuchos::rcp(
new map_type(rowMap->getGlobalNumElements(), newRowElems, rowMap->getIndexBase(), rowMap->getComm()));
475 auto colMap = A_->getColMap();
476 auto numElems = colMap->getLocalNumElements();
477 auto colElems = colMap->getLocalElementList();
479 Teuchos::Array<global_ordinal_type> newColElems(colElems.size());
480 for (
size_t i = 0; i < numElems; i++)
481 newColElems[i] = colElems[numElems-1 - i];
483 newColMap = Teuchos::rcp(
new map_type(colMap->getGlobalNumElements(), newColElems, colMap->getIndexBase(), colMap->getComm()));
487 local_matrix_type newLocalMatrix(
"Upermuted", numRows, numCols, numNnz, newval, newptr, newind);
489 A_crs_ = Teuchos::rcp(
new crs_matrix_type(newLocalMatrix, newRowMap, newColMap, A_crs_->getDomainMap(), A_crs_->getRangeMap()));
491 isInternallyChanged_ =
true;
496 auto newLclTriStructure =
497 determineLocalTriangularStructure (newLocalMatrix.graph,
498 newRowMap->getLocalMap (),
499 newColMap->getLocalMap (),
500 ignoreMapsForTriStructure);
501 const LO newLclNumRows = newRowMap->getLocalNumElements ();
502 this->diag_ = (newLclTriStructure.diagCount < newLclNumRows) ?
"U" :
"N";
503 this->uplo_ = newLclTriStructure.couldBeLowerTriangular ?
"L" :
504 (newLclTriStructure.couldBeUpperTriangular ?
"U" :
"N");
507 if (Teuchos::nonnull (htsImpl_))
509 htsImpl_->initialize (*A_crs_);
510 isInternallyChanged_ =
true;
513 const bool ksptrsv_valid_uplo = (this->uplo_ !=
"N");
514 if (Teuchos::nonnull(kh_) && ksptrsv_valid_uplo && this->diag_ !=
"U")
516 this->isKokkosKernelsSptrsv_ =
true;
520 this->isKokkosKernelsSptrsv_ =
false;
523 isInitialized_ =
true;
527template<
class MatrixType>
532 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::compute: ";
533 if (! out_.is_null ()) {
534 *out_ <<
">>> DEBUG " << prefix << std::endl;
537 TEUCHOS_TEST_FOR_EXCEPTION
538 (A_.is_null (), std::runtime_error, prefix <<
"You must call "
539 "setMatrix() with a nonnull input matrix before you may call "
540 "initialize() or compute().");
541 TEUCHOS_TEST_FOR_EXCEPTION
542 (A_crs_.is_null (), std::logic_error, prefix <<
"A_ is nonnull, but "
543 "A_crs_ is null. Please report this bug to the Ifpack2 developers.");
545 TEUCHOS_TEST_FOR_EXCEPTION
546 (! A_crs_->isFillComplete (), std::runtime_error,
"If you call this "
547 "method, the matrix must be fill complete. It is not.");
549 if (! isInitialized_) {
552 TEUCHOS_TEST_FOR_EXCEPTION
553 (! isInitialized_, std::logic_error, prefix <<
"initialize() should have "
554 "been called by this point, but isInitialized_ is false. "
555 "Please report this bug to the Ifpack2 developers.");
558 if (Teuchos::nonnull (htsImpl_))
559 htsImpl_->compute (*A_crs_, out_);
561 if (Teuchos::nonnull(kh_) && this->isKokkosKernelsSptrsv_)
563 auto A_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
564 auto Alocal = A_crs->getLocalMatrixDevice();
565 auto ptr = Alocal.graph.row_map;
566 auto ind = Alocal.graph.entries;
567 auto val = Alocal.values;
569 auto numRows = Alocal.numRows();
570 const bool is_lower_tri = (this->uplo_ ==
"L") ?
true :
false;
573 kh_->destroy_sptrsv_handle();
574#if defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) && defined(KOKKOS_ENABLE_CUDA)
577 if (std::is_same<Kokkos::Cuda, HandleExecSpace>::value &&
578 std::is_same<int, local_ordinal_type>::value &&
579 (std::is_same<scalar_type, float>::value ||
580 std::is_same<scalar_type, double>::value ||
581 std::is_same<
scalar_type, Kokkos::complex<float>>::value ||
582 std::is_same<
scalar_type, Kokkos::complex<double>>::value))
584 kh_->create_sptrsv_handle(KokkosSparse::Experimental::SPTRSVAlgorithm::SPTRSV_CUSPARSE, numRows, is_lower_tri);
589 kh_->create_sptrsv_handle(KokkosSparse::Experimental::SPTRSVAlgorithm::SEQLVLSCHD_TP1, numRows, is_lower_tri);
591 KokkosSparse::Experimental::sptrsv_symbolic(kh_.getRawPtr(), ptr, ind, val);
599template<
class MatrixType>
605 Teuchos::ETransp mode,
611 using Teuchos::rcpFromRef;
613 typedef Teuchos::ScalarTraits<ST> STS;
614 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::apply: ";
615 if (! out_.is_null ()) {
616 *out_ <<
">>> DEBUG " << prefix;
617 if (A_crs_.is_null ()) {
618 *out_ <<
"A_crs_ is null!" << std::endl;
621 Teuchos::RCP<const crs_matrix_type> A_crs =
622 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
623 const std::string uplo = this->uplo_;
624 const std::string trans = (mode == Teuchos::CONJ_TRANS) ?
"C" :
625 (mode == Teuchos::TRANS ?
"T" :
"N");
626 const std::string diag = this->diag_;
627 *out_ <<
"uplo=\"" << uplo
628 <<
"\", trans=\"" << trans
629 <<
"\", diag=\"" << diag <<
"\"" << std::endl;
633 TEUCHOS_TEST_FOR_EXCEPTION
634 (! isComputed (), std::runtime_error, prefix <<
"If compute() has not yet "
635 "been called, or if you have changed the matrix via setMatrix(), you must "
636 "call compute() before you may call this method.");
639 TEUCHOS_TEST_FOR_EXCEPTION
640 (A_.is_null (), std::logic_error, prefix <<
"A_ is null. "
641 "Please report this bug to the Ifpack2 developers.");
642 TEUCHOS_TEST_FOR_EXCEPTION
643 (A_crs_.is_null (), std::logic_error, prefix <<
"A_crs_ is null. "
644 "Please report this bug to the Ifpack2 developers.");
647 TEUCHOS_TEST_FOR_EXCEPTION
648 (! A_crs_->isFillComplete (), std::runtime_error,
"If you call this "
649 "method, the matrix must be fill complete. It is not. This means that "
650 " you must have called resumeFill() on the matrix before calling apply(). "
651 "This is NOT allowed. Note that this class may use the matrix's data in "
652 "place without copying it. Thus, you cannot change the matrix and expect "
653 "the solver to stay the same. If you have changed the matrix, first call "
654 "fillComplete() on it, then call compute() on this object, before you call"
655 " apply(). You do NOT need to call setMatrix, as long as the matrix "
656 "itself (that is, its address in memory) is the same.");
658 auto G = A_crs_->getGraph ();
659 TEUCHOS_TEST_FOR_EXCEPTION
660 (G.is_null (), std::logic_error, prefix <<
"A_ and A_crs_ are nonnull, "
661 "but A_crs_'s RowGraph G is null. "
662 "Please report this bug to the Ifpack2 developers.");
663 auto importer = G->getImporter ();
664 auto exporter = G->getExporter ();
666 if (! importer.is_null ()) {
667 if (X_colMap_.is_null () || X_colMap_->getNumVectors () != X.getNumVectors ()) {
668 X_colMap_ = rcp (
new MV (importer->getTargetMap (), X.getNumVectors ()));
671 X_colMap_->putScalar (STS::zero ());
676 X_colMap_->doImport (X, *importer, Tpetra::ZERO);
678 RCP<const MV> X_cur = importer.is_null () ? rcpFromRef (X) :
679 Teuchos::rcp_const_cast<const MV> (X_colMap_);
681 if (! exporter.is_null ()) {
682 if (Y_rowMap_.is_null () || Y_rowMap_->getNumVectors () != Y.getNumVectors ()) {
683 Y_rowMap_ = rcp (
new MV (exporter->getSourceMap (), Y.getNumVectors ()));
686 Y_rowMap_->putScalar (STS::zero ());
688 Y_rowMap_->doExport (Y, *importer, Tpetra::ADD);
690 RCP<MV> Y_cur = exporter.is_null () ? rcpFromRef (Y) : Y_rowMap_;
692 localApply (*X_cur, *Y_cur, mode, alpha, beta);
694 if (! exporter.is_null ()) {
695 Y.putScalar (STS::zero ());
696 Y.doExport (*Y_cur, *exporter, Tpetra::ADD);
702template<
class MatrixType>
707 const Teuchos::ETransp mode)
const
709 using Teuchos::CONJ_TRANS;
710 using Teuchos::NO_TRANS;
711 using Teuchos::TRANS;
712 const char tfecfFuncName[] =
"localTriangularSolve: ";
714 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
715 (! A_crs_->isFillComplete (), std::runtime_error,
716 "The matrix is not fill complete.");
717 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
718 (! X.isConstantStride () || ! Y.isConstantStride (), std::invalid_argument,
719 "X and Y must be constant stride.");
720 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
721 ( A_crs_->getLocalNumRows() > 0 && this->uplo_ ==
"N", std::runtime_error,
722 "The matrix is neither upper triangular or lower triangular. "
723 "You may only call this method if the matrix is triangular. "
724 "Remember that this is a local (per MPI process) property, and that "
725 "Tpetra only knows how to do a local (per process) triangular solve.");
726 using STS = Teuchos::ScalarTraits<scalar_type>;
727 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
728 (STS::isComplex && mode == TRANS, std::logic_error,
"This method does "
729 "not currently support non-conjugated transposed solve (mode == "
730 "Teuchos::TRANS) for complex scalar types.");
732 const std::string uplo = this->uplo_;
733 const std::string trans = (mode == Teuchos::CONJ_TRANS) ?
"C" :
734 (mode == Teuchos::
TRANS ?
"T" :
"N");
736 if (Teuchos::nonnull(kh_) && this->isKokkosKernelsSptrsv_ && trans ==
"N")
738 auto A_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type> (this->A_);
739 auto A_lclk = A_crs->getLocalMatrixDevice ();
740 auto ptr = A_lclk.graph.row_map;
741 auto ind = A_lclk.graph.entries;
742 auto val = A_lclk.values;
744 const size_t numVecs = std::min (X.getNumVectors (), Y.getNumVectors ());
746 for (
size_t j = 0; j < numVecs; ++j) {
747 auto X_j = X.getVectorNonConst (j);
748 auto Y_j = Y.getVector (j);
749 auto X_lcl = X_j->getLocalViewDevice (Tpetra::Access::ReadWrite);
750 auto Y_lcl = Y_j->getLocalViewDevice (Tpetra::Access::ReadOnly);
751 auto X_lcl_1d = Kokkos::subview (X_lcl, Kokkos::ALL (), 0);
752 auto Y_lcl_1d = Kokkos::subview (Y_lcl, Kokkos::ALL (), 0);
753 KokkosSparse::Experimental::sptrsv_solve(kh_.getRawPtr(), ptr, ind, val, Y_lcl_1d, X_lcl_1d);
755 typename k_handle::HandleExecSpace().fence();
760 const std::string diag = this->diag_;
765 auto A_lcl = this->A_crs_->getLocalMatrixHost ();
767 if (X.isConstantStride () && Y.isConstantStride ()) {
768 auto X_lcl = X.getLocalViewHost (Tpetra::Access::ReadWrite);
769 auto Y_lcl = Y.getLocalViewHost (Tpetra::Access::ReadOnly);
770 KokkosSparse::trsv (uplo.c_str (), trans.c_str (), diag.c_str (),
771 A_lcl, Y_lcl, X_lcl);
774 const size_t numVecs =
775 std::min (X.getNumVectors (), Y.getNumVectors ());
776 for (
size_t j = 0; j < numVecs; ++j) {
777 auto X_j = X.getVectorNonConst (j);
778 auto Y_j = Y.getVector (j);
779 auto X_lcl = X_j->getLocalViewHost (Tpetra::Access::ReadWrite);
780 auto Y_lcl = Y_j->getLocalViewHost (Tpetra::Access::ReadOnly);
781 KokkosSparse::trsv (uplo.c_str (), trans.c_str (),
782 diag.c_str (), A_lcl, Y_lcl, X_lcl);
788template<
class MatrixType>
790LocalSparseTriangularSolver<MatrixType>::
791localApply (
const MV& X,
793 const Teuchos::ETransp mode,
794 const scalar_type& alpha,
795 const scalar_type& beta)
const
797 if (mode == Teuchos::NO_TRANS && Teuchos::nonnull (htsImpl_) &&
798 htsImpl_->isComputed ()) {
799 htsImpl_->localApply (X, Y, mode, alpha, beta);
804 typedef scalar_type ST;
805 typedef Teuchos::ScalarTraits<ST> STS;
807 if (beta == STS::zero ()) {
808 if (alpha == STS::zero ()) {
809 Y.putScalar (STS::zero ());
812 this->localTriangularSolve (X, Y, mode);
813 if (alpha != STS::one ()) {
819 if (alpha == STS::zero ()) {
823 MV Y_tmp (Y, Teuchos::Copy);
824 this->localTriangularSolve (X, Y_tmp, mode);
825 Y.update (alpha, Y_tmp, beta);
831template <
class MatrixType>
835 return numInitialize_;
838template <
class MatrixType>
845template <
class MatrixType>
852template <
class MatrixType>
856 return initializeTime_;
859template<
class MatrixType>
866template<
class MatrixType>
873template <
class MatrixType>
878 std::ostringstream os;
883 os <<
"\"Ifpack2::LocalSparseTriangularSolver\": {";
884 if (this->getObjectLabel () !=
"") {
885 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
887 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", "
888 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
891 os <<
"Matrix: null";
894 os <<
"Matrix dimensions: ["
895 << A_->getGlobalNumRows () <<
", "
896 << A_->getGlobalNumCols () <<
"]"
897 <<
", Number of nonzeros: " << A_->getGlobalNumEntries();
900 if (Teuchos::nonnull (htsImpl_))
901 os <<
", HTS computed: " << (htsImpl_->isComputed () ?
"true" :
"false");
906template <
class MatrixType>
908describe (Teuchos::FancyOStream& out,
909 const Teuchos::EVerbosityLevel verbLevel)
const
914 const Teuchos::EVerbosityLevel vl
915 = (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
917 if (vl != Teuchos::VERB_NONE) {
922 auto comm = A_.is_null () ?
923 Tpetra::getDefaultComm () :
928 if (! comm.is_null () && comm->getRank () == 0) {
930 Teuchos::OSTab tab0 (out);
933 out <<
"\"Ifpack2::LocalSparseTriangularSolver\":" << endl;
934 Teuchos::OSTab tab1 (out);
935 out <<
"Scalar: " << Teuchos::TypeNameTraits<scalar_type>::name () << endl
936 <<
"LocalOrdinal: " << Teuchos::TypeNameTraits<local_ordinal_type>::name () << endl
937 <<
"GlobalOrdinal: " << Teuchos::TypeNameTraits<global_ordinal_type>::name () << endl
938 <<
"Node: " << Teuchos::TypeNameTraits<node_type>::name () << endl;
943template <
class MatrixType>
944Teuchos::RCP<const typename LocalSparseTriangularSolver<MatrixType>::map_type>
948 TEUCHOS_TEST_FOR_EXCEPTION
949 (A_.is_null (), std::runtime_error,
950 "Ifpack2::LocalSparseTriangularSolver::getDomainMap: "
951 "The matrix is null. Please call setMatrix() with a nonnull input "
952 "before calling this method.");
953 return A_->getDomainMap ();
956template <
class MatrixType>
957Teuchos::RCP<const typename LocalSparseTriangularSolver<MatrixType>::map_type>
961 TEUCHOS_TEST_FOR_EXCEPTION
962 (A_.is_null (), std::runtime_error,
963 "Ifpack2::LocalSparseTriangularSolver::getRangeMap: "
964 "The matrix is null. Please call setMatrix() with a nonnull input "
965 "before calling this method.");
966 return A_->getRangeMap ();
969template<
class MatrixType>
971setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
973 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::setMatrix: ";
979 if (A.getRawPtr () != A_.getRawPtr () || isInternallyChanged_) {
981 TEUCHOS_TEST_FOR_EXCEPTION
982 (! A.is_null () && A->getComm ()->getSize () == 1 &&
983 A->getLocalNumRows () != A->getLocalNumCols (),
984 std::runtime_error, prefix <<
"If A's communicator only contains one "
985 "process, then A must be square. Instead, you provided a matrix A with "
986 << A->getLocalNumRows () <<
" rows and " << A->getLocalNumCols ()
992 isInitialized_ =
false;
996 A_crs_ = Teuchos::null;
1000 Teuchos::RCP<const crs_matrix_type> A_crs =
1001 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A);
1002 TEUCHOS_TEST_FOR_EXCEPTION
1003 (A_crs.is_null (), std::invalid_argument, prefix <<
1004 "The input matrix A is not a Tpetra::CrsMatrix.");
1009 if (Teuchos::nonnull (htsImpl_))
1017#if defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) && defined(KOKKOS_ENABLE_CUDA) && (CUDA_VERSION >= 11030)
1018 isComputed_ =
false;
1024#define IFPACK2_LOCALSPARSETRIANGULARSOLVER_INSTANT(S,LO,GO,N) \
1025 template class Ifpack2::LocalSparseTriangularSolver< Tpetra::RowMatrix<S, LO, GO, N> >;
"Preconditioner" that solves local sparse triangular systems.
Definition Ifpack2_LocalSparseTriangularSolver_decl.hpp:88
bool isComputed() const
Return true if compute() has been called.
Definition Ifpack2_LocalSparseTriangularSolver_decl.hpp:222
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print this object with given verbosity to the given output stream.
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:908
std::string description() const
A one-line description of this object.
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:876
Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Specialization of Tpetra::Map used by this class.
Definition Ifpack2_LocalSparseTriangularSolver_decl.hpp:102
LocalSparseTriangularSolver()
Constructor that takes no input matrix.
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:282
MatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the input matrix.
Definition Ifpack2_LocalSparseTriangularSolver_decl.hpp:93
void setParameters(const Teuchos::ParameterList ¶ms)
Set this object's parameters.
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:330
double getComputeTime() const
Return the time spent in compute().
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:862
void initialize()
"Symbolic" phase of setup
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:374
MatrixType::global_ordinal_type global_ordinal_type
Type of the global indices of the input matrix.
Definition Ifpack2_LocalSparseTriangularSolver_decl.hpp:95
int getNumCompute() const
Return the number of calls to compute().
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:841
double getInitializeTime() const
Return the time spent in initialize().
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:855
Teuchos::RCP< const map_type > getDomainMap() const
The domain of this operator.
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:946
MatrixType::scalar_type scalar_type
Type of the entries of the input matrix.
Definition Ifpack2_LocalSparseTriangularSolver_decl.hpp:91
MatrixType::node_type node_type
Node type of the input matrix.
Definition Ifpack2_LocalSparseTriangularSolver_decl.hpp:97
int getNumInitialize() const
Return the number of calls to initialize().
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:834
Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > crs_matrix_type
Specialization of Tpetra::CrsMatrix used by this class.
Definition Ifpack2_LocalSparseTriangularSolver_decl.hpp:108
Teuchos::RCP< const map_type > getRangeMap() const
The range of this operator.
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:959
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set this preconditioner's matrix.
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:971
virtual ~LocalSparseTriangularSolver()
Destructor (virtual for memory safety).
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:319
double getApplyTime() const
Return the time spent in apply().
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:869
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, and put the result in Y.
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:601
int getNumApply() const
Return the number of calls to apply().
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:848
void compute()
"Numeric" phase of setup
Definition Ifpack2_LocalSparseTriangularSolver_def.hpp:530
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74