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" 59 struct TrisolverType {
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;
79 template<
class MatrixType>
80 class 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 Teuchos::ArrayRCP<const size_t> rowptr;
123 Teuchos::ArrayRCP<const local_ordinal_type> colidx;
124 Teuchos::ArrayRCP<const scalar_type> val;
125 T_in.getAllValues(rowptr, colidx, val);
128 Teuchos::RCP<HtsCrsMatrix> T_hts = Teuchos::rcpWithDealloc(
129 HTST::make_CrsMatrix(rowptr.size() - 1,
130 rowptr.getRawPtr(), colidx.getRawPtr(), val.getRawPtr(),
131 transpose_, conjugate_),
132 HtsCrsMatrixDeleter());
134 if (Teuchos::nonnull(Timpl_)) {
136 HTST::reprocess_numeric(Timpl_.get(), T_hts.get());
139 if (T_in.getCrsGraph().is_null()) {
140 if (Teuchos::nonnull(out))
141 *out <<
"HTS compute failed because T_in.getCrsGraph().is_null().\n";
144 if ( ! T_in.getCrsGraph()->isSorted()) {
145 if (Teuchos::nonnull(out))
146 *out <<
"HTS compute failed because ! T_in.getCrsGraph().isSorted().\n";
149 if ( ! T_in.isStorageOptimized()) {
150 if (Teuchos::nonnull(out))
151 *out <<
"HTS compute failed because ! T_in.isStorageOptimized().\n";
155 typename HTST::PreprocessArgs args;
156 args.T = T_hts.get();
159 args.nthreads = omp_get_max_threads();
163 args.save_for_reprocess =
true;
164 typename HTST::Options opts;
165 opts.levelset_block_size = levelset_block_size_;
166 args.options = &opts;
169 Timpl_ = Teuchos::rcpWithDealloc(HTST::preprocess(args), TImplDeleter());
170 }
catch (
const std::exception& e) {
171 if (Teuchos::nonnull(out))
172 *out <<
"HTS preprocess threw: " << e.what() <<
"\n";
181 #ifdef HAVE_IFPACK2_SHYLU_NODEHTS 182 return Teuchos::nonnull(Timpl_);
189 void localApply (
const MV& X, MV& Y,
190 const Teuchos::ETransp ,
196 #ifdef HAVE_IFPACK2_SHYLU_NODEHTS 197 const auto& X_view = X.getLocalViewHost (Tpetra::Access::ReadOnly);
198 const auto& Y_view = Y.getLocalViewHost (Tpetra::Access::ReadWrite);
201 HTST::reset_max_nrhs(Timpl_.get(), X_view.extent(1));
203 HTST::solve_omp(Timpl_.get(),
205 reinterpret_cast<const scalar_type*
>(X_view.data()),
214 #ifdef HAVE_IFPACK2_SHYLU_NODEHTS 215 typedef ::Experimental::HTS<local_ordinal_type, size_t, scalar_type> HTST;
216 typedef typename HTST::Impl TImpl;
217 typedef typename HTST::CrsMatrix HtsCrsMatrix;
219 struct TImplDeleter {
220 void free (TImpl* impl) {
221 HTST::delete_Impl(impl);
225 struct HtsCrsMatrixDeleter {
226 void free (HtsCrsMatrix* T) {
227 HTST::delete_CrsMatrix(T);
231 Teuchos::RCP<TImpl> Timpl_;
232 bool transpose_, conjugate_;
233 int levelset_block_size_;
237 template<
class MatrixType>
245 if (! A.is_null ()) {
246 Teuchos::RCP<const crs_matrix_type> A_crs =
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.");
256 template<
class MatrixType>
259 const Teuchos::RCP<Teuchos::FancyOStream>& out) :
264 if (! out_.is_null ()) {
265 *out_ <<
">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor" 270 if (! A.is_null ()) {
271 Teuchos::RCP<const crs_matrix_type> A_crs =
273 TEUCHOS_TEST_FOR_EXCEPTION
274 (A_crs.is_null (), std::invalid_argument,
275 "Ifpack2::LocalSparseTriangularSolver constructor: " 276 "The input matrix A is not a Tpetra::CrsMatrix.");
281 template<
class MatrixType>
288 template<
class MatrixType>
294 if (! out_.is_null ()) {
295 *out_ <<
">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor" 300 template<
class MatrixType>
303 isInitialized_ =
false;
305 reverseStorage_ =
false;
306 isInternallyChanged_ =
false;
310 initializeTime_ = 0.0;
313 isKokkosKernelsSptrsv_ =
false;
318 template<
class MatrixType>
322 if (Teuchos::nonnull (kh_))
324 kh_->destroy_sptrsv_handle();
328 template<
class MatrixType>
334 using Teuchos::ParameterList;
335 using Teuchos::Array;
337 Details::TrisolverType::Enum trisolverType = Details::TrisolverType::Internal;
339 static const char typeName[] =
"trisolver: type";
341 if ( ! pl.isType<std::string>(typeName))
break;
344 Array<std::string> trisolverTypeStrs;
345 Array<Details::TrisolverType::Enum> trisolverTypeEnums;
346 Details::TrisolverType::loadPLTypeOption (trisolverTypeStrs, trisolverTypeEnums);
347 Teuchos::StringToIntegralParameterEntryValidator<Details::TrisolverType::Enum>
348 s2i(trisolverTypeStrs (), trisolverTypeEnums (), typeName,
false);
350 trisolverType = s2i.getIntegralValue(pl.get<std::string>(typeName));
353 if (trisolverType == Details::TrisolverType::HTS) {
354 htsImpl_ = Teuchos::rcp (
new HtsImpl ());
355 htsImpl_->setParameters (pl);
358 if (trisolverType == Details::TrisolverType::KSPTRSV) {
359 kh_ = Teuchos::rcp (
new k_handle());
362 if (pl.isParameter(
"trisolver: reverse U"))
363 reverseStorage_ = pl.get<
bool>(
"trisolver: reverse U");
365 TEUCHOS_TEST_FOR_EXCEPTION
366 (reverseStorage_ && (trisolverType == Details::TrisolverType::HTS || trisolverType == Details::TrisolverType::KSPTRSV),
367 std::logic_error,
"Ifpack2::LocalSparseTriangularSolver::setParameters: " 368 "You are not allowed to enable both HTS or KSPTRSV and the \"trisolver: reverse U\" " 369 "options. See GitHub issue #2647.");
372 template<
class MatrixType>
377 using Tpetra::Details::determineLocalTriangularStructure;
380 using local_matrix_type =
typename crs_matrix_type::local_matrix_device_type;
383 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::initialize: ";
384 if (! out_.is_null ()) {
385 *out_ <<
">>> DEBUG " << prefix << std::endl;
388 TEUCHOS_TEST_FOR_EXCEPTION
389 (A_.is_null (), std::runtime_error, prefix <<
"You must call " 390 "setMatrix() with a nonnull input matrix before you may call " 391 "initialize() or compute().");
392 if (A_crs_.is_null ()) {
394 TEUCHOS_TEST_FOR_EXCEPTION
395 (A_crs.get () ==
nullptr, std::invalid_argument,
396 prefix <<
"The input matrix A is not a Tpetra::CrsMatrix.");
399 auto G = A_crs_->getGraph ();
400 TEUCHOS_TEST_FOR_EXCEPTION
401 (G.is_null (), std::logic_error, prefix <<
"A_ and A_crs_ are nonnull, " 402 "but A_crs_'s RowGraph G is null. " 403 "Please report this bug to the Ifpack2 developers.");
407 TEUCHOS_TEST_FOR_EXCEPTION
408 (! G->isFillComplete (), std::runtime_error,
"If you call this method, " 409 "the matrix's graph must be fill complete. It is not.");
412 constexpr
bool ignoreMapsForTriStructure =
true;
413 auto lclTriStructure = [&] {
414 auto lclMatrix = A_crs_->getLocalMatrixDevice ();
415 auto lclRowMap = A_crs_->getRowMap ()->getLocalMap ();
416 auto lclColMap = A_crs_->getColMap ()->getLocalMap ();
418 determineLocalTriangularStructure (lclMatrix.graph,
421 ignoreMapsForTriStructure);
422 const LO lclNumRows = lclRowMap.getNodeNumElements ();
423 this->diag_ = (lclTriStruct.diagCount < lclNumRows) ?
"U" :
"N";
424 this->uplo_ = lclTriStruct.couldBeLowerTriangular ?
"L" :
425 (lclTriStruct.couldBeUpperTriangular ?
"U" :
"N");
429 if (reverseStorage_ && lclTriStructure.couldBeUpperTriangular &&
430 htsImpl_.is_null ()) {
432 auto Alocal = A_crs_->getLocalMatrixDevice();
433 auto ptr = Alocal.graph.row_map;
434 auto ind = Alocal.graph.entries;
435 auto val = Alocal.values;
437 auto numRows = Alocal.numRows();
438 auto numCols = Alocal.numCols();
439 auto numNnz = Alocal.nnz();
441 typename decltype(ptr)::non_const_type newptr (
"ptr", ptr.extent (0));
442 typename decltype(ind)::non_const_type newind (
"ind", ind.extent (0));
443 decltype(val) newval (
"val", val.extent (0));
446 typename crs_matrix_type::execution_space().fence();
449 auto A_r = Alocal.row(numRows-1 - row);
451 auto numEnt = A_r.length;
453 newind(rowStart + k) = numCols-1 - A_r.colidx(numEnt-1 - k);
454 newval(rowStart + k) = A_r.value (numEnt-1 - k);
457 newptr(row+1) = rowStart;
459 typename crs_matrix_type::execution_space().fence();
462 using map_type =
typename crs_matrix_type::map_type;
463 Teuchos::RCP<map_type> newRowMap, newColMap;
466 auto rowMap = A_->getRowMap();
467 auto numElems = rowMap->getNodeNumElements();
468 auto rowElems = rowMap->getNodeElementList();
470 Teuchos::Array<global_ordinal_type> newRowElems(rowElems.size());
471 for (
size_t i = 0; i < numElems; i++)
472 newRowElems[i] = rowElems[numElems-1 - i];
474 newRowMap = Teuchos::rcp(
new map_type(rowMap->getGlobalNumElements(), newRowElems, rowMap->getIndexBase(), rowMap->getComm()));
478 auto colMap = A_->getColMap();
479 auto numElems = colMap->getNodeNumElements();
480 auto colElems = colMap->getNodeElementList();
482 Teuchos::Array<global_ordinal_type> newColElems(colElems.size());
483 for (
size_t i = 0; i < numElems; i++)
484 newColElems[i] = colElems[numElems-1 - i];
486 newColMap = Teuchos::rcp(
new map_type(colMap->getGlobalNumElements(), newColElems, colMap->getIndexBase(), colMap->getComm()));
490 local_matrix_type newLocalMatrix(
"Upermuted", numRows, numCols, numNnz, newval, newptr, newind);
492 A_crs_ = Teuchos::rcp(
new crs_matrix_type(newLocalMatrix, newRowMap, newColMap, A_crs_->getDomainMap(), A_crs_->getRangeMap()));
494 isInternallyChanged_ =
true;
499 auto newLclTriStructure =
500 determineLocalTriangularStructure (newLocalMatrix.graph,
501 newRowMap->getLocalMap (),
502 newColMap->getLocalMap (),
503 ignoreMapsForTriStructure);
504 const LO newLclNumRows = newRowMap->getNodeNumElements ();
505 this->diag_ = (newLclTriStructure.diagCount < newLclNumRows) ?
"U" :
"N";
506 this->uplo_ = newLclTriStructure.couldBeLowerTriangular ?
"L" :
507 (newLclTriStructure.couldBeUpperTriangular ?
"U" :
"N");
510 if (Teuchos::nonnull (htsImpl_))
512 htsImpl_->initialize (*A_crs_);
513 isInternallyChanged_ =
true;
516 const bool ksptrsv_valid_uplo = (this->uplo_ !=
"N");
517 if (Teuchos::nonnull(kh_) && ksptrsv_valid_uplo && this->diag_ !=
"U")
519 this->isKokkosKernelsSptrsv_ =
true;
523 this->isKokkosKernelsSptrsv_ =
false;
526 isInitialized_ =
true;
530 template<
class MatrixType>
535 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::compute: ";
536 if (! out_.is_null ()) {
537 *out_ <<
">>> DEBUG " << prefix << std::endl;
540 TEUCHOS_TEST_FOR_EXCEPTION
541 (A_.is_null (), std::runtime_error, prefix <<
"You must call " 542 "setMatrix() with a nonnull input matrix before you may call " 543 "initialize() or compute().");
544 TEUCHOS_TEST_FOR_EXCEPTION
545 (A_crs_.is_null (), std::logic_error, prefix <<
"A_ is nonnull, but " 546 "A_crs_ is null. Please report this bug to the Ifpack2 developers.");
548 TEUCHOS_TEST_FOR_EXCEPTION
549 (! A_crs_->isFillComplete (), std::runtime_error,
"If you call this " 550 "method, the matrix must be fill complete. It is not.");
552 if (! isInitialized_) {
555 TEUCHOS_TEST_FOR_EXCEPTION
556 (! isInitialized_, std::logic_error, prefix <<
"initialize() should have " 557 "been called by this point, but isInitialized_ is false. " 558 "Please report this bug to the Ifpack2 developers.");
560 if (Teuchos::nonnull (htsImpl_))
561 htsImpl_->compute (*A_crs_, out_);
563 if (Teuchos::nonnull(kh_) && this->isKokkosKernelsSptrsv_)
566 auto Alocal = A_crs->getLocalMatrixDevice();
567 auto ptr = Alocal.graph.row_map;
568 auto ind = Alocal.graph.entries;
569 auto val = Alocal.values;
571 auto numRows = Alocal.numRows();
572 const bool is_lower_tri = (this->uplo_ ==
"L") ?
true :
false;
575 kh_->destroy_sptrsv_handle();
576 #if defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) && defined(KOKKOS_ENABLE_CUDA) 578 if (std::is_same<Kokkos::Cuda, HandleExecSpace>::value && std::is_same<int,local_ordinal_type >::value)
580 kh_->create_sptrsv_handle(KokkosSparse::Experimental::SPTRSVAlgorithm::SPTRSV_CUSPARSE, numRows, is_lower_tri);
585 kh_->create_sptrsv_handle(KokkosSparse::Experimental::SPTRSVAlgorithm::SEQLVLSCHD_TP1, numRows, is_lower_tri);
587 KokkosSparse::Experimental::sptrsv_symbolic(kh_.getRawPtr(), ptr, ind, val);
594 template<
class MatrixType>
600 Teuchos::ETransp mode,
606 using Teuchos::rcpFromRef;
608 typedef Teuchos::ScalarTraits<ST> STS;
609 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::apply: ";
610 if (! out_.is_null ()) {
611 *out_ <<
">>> DEBUG " << prefix;
612 if (A_crs_.is_null ()) {
613 *out_ <<
"A_crs_ is null!" << std::endl;
616 Teuchos::RCP<const crs_matrix_type> A_crs =
618 const std::string uplo = this->uplo_;
619 const std::string trans = (mode == Teuchos::CONJ_TRANS) ?
"C" :
620 (mode == Teuchos::TRANS ?
"T" :
"N");
621 const std::string diag = this->diag_;
622 *out_ <<
"uplo=\"" << uplo
623 <<
"\", trans=\"" << trans
624 <<
"\", diag=\"" << diag <<
"\"" << std::endl;
628 TEUCHOS_TEST_FOR_EXCEPTION
629 (! isComputed (), std::runtime_error, prefix <<
"If compute() has not yet " 630 "been called, or if you have changed the matrix via setMatrix(), you must " 631 "call compute() before you may call this method.");
634 TEUCHOS_TEST_FOR_EXCEPTION
635 (A_.is_null (), std::logic_error, prefix <<
"A_ is null. " 636 "Please report this bug to the Ifpack2 developers.");
637 TEUCHOS_TEST_FOR_EXCEPTION
638 (A_crs_.is_null (), std::logic_error, prefix <<
"A_crs_ is null. " 639 "Please report this bug to the Ifpack2 developers.");
642 TEUCHOS_TEST_FOR_EXCEPTION
643 (! A_crs_->isFillComplete (), std::runtime_error,
"If you call this " 644 "method, the matrix must be fill complete. It is not. This means that " 645 " you must have called resumeFill() on the matrix before calling apply(). " 646 "This is NOT allowed. Note that this class may use the matrix's data in " 647 "place without copying it. Thus, you cannot change the matrix and expect " 648 "the solver to stay the same. If you have changed the matrix, first call " 649 "fillComplete() on it, then call compute() on this object, before you call" 650 " apply(). You do NOT need to call setMatrix, as long as the matrix " 651 "itself (that is, its address in memory) is the same.");
653 auto G = A_crs_->getGraph ();
654 TEUCHOS_TEST_FOR_EXCEPTION
655 (G.is_null (), std::logic_error, prefix <<
"A_ and A_crs_ are nonnull, " 656 "but A_crs_'s RowGraph G is null. " 657 "Please report this bug to the Ifpack2 developers.");
658 auto importer = G->getImporter ();
659 auto exporter = G->getExporter ();
661 if (! importer.is_null ()) {
662 if (X_colMap_.is_null () || X_colMap_->getNumVectors () != X.getNumVectors ()) {
663 X_colMap_ = rcp (
new MV (importer->getTargetMap (), X.getNumVectors ()));
666 X_colMap_->putScalar (STS::zero ());
671 X_colMap_->doImport (X, *importer, Tpetra::ZERO);
673 RCP<const MV> X_cur = importer.is_null () ? rcpFromRef (X) :
674 Teuchos::rcp_const_cast<
const MV> (X_colMap_);
676 if (! exporter.is_null ()) {
677 if (Y_rowMap_.is_null () || Y_rowMap_->getNumVectors () != Y.getNumVectors ()) {
678 Y_rowMap_ = rcp (
new MV (exporter->getSourceMap (), Y.getNumVectors ()));
681 Y_rowMap_->putScalar (STS::zero ());
683 Y_rowMap_->doExport (Y, *importer, Tpetra::ADD);
685 RCP<MV> Y_cur = exporter.is_null () ? rcpFromRef (Y) : Y_rowMap_;
687 localApply (*X_cur, *Y_cur, mode, alpha, beta);
689 if (! exporter.is_null ()) {
690 Y.putScalar (STS::zero ());
691 Y.doExport (*Y_cur, *exporter, Tpetra::ADD);
697 template<
class MatrixType>
702 const Teuchos::ETransp mode)
const 704 using Teuchos::CONJ_TRANS;
705 using Teuchos::NO_TRANS;
706 using Teuchos::TRANS;
707 const char tfecfFuncName[] =
"localTriangularSolve: ";
709 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
710 (! A_crs_->isFillComplete (), std::runtime_error,
711 "The matrix is not fill complete.");
712 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
713 (! X.isConstantStride () || ! Y.isConstantStride (), std::invalid_argument,
714 "X and Y must be constant stride.");
715 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
716 ( A_crs_->getNodeNumRows() > 0 && this->uplo_ ==
"N", std::runtime_error,
717 "The matrix is neither upper triangular or lower triangular. " 718 "You may only call this method if the matrix is triangular. " 719 "Remember that this is a local (per MPI process) property, and that " 720 "Tpetra only knows how to do a local (per process) triangular solve.");
721 using STS = Teuchos::ScalarTraits<scalar_type>;
722 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
723 (STS::isComplex && mode == TRANS, std::logic_error,
"This method does " 724 "not currently support non-conjugated transposed solve (mode == " 725 "Teuchos::TRANS) for complex scalar types.");
727 const std::string uplo = this->uplo_;
728 const std::string trans = (mode == Teuchos::CONJ_TRANS) ?
"C" :
729 (mode == Teuchos::TRANS ?
"T" :
"N");
731 if (Teuchos::nonnull(kh_) && this->isKokkosKernelsSptrsv_ && trans ==
"N")
733 auto A_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (this->A_);
734 auto A_lclk = A_crs->getLocalMatrixDevice ();
735 auto ptr = A_lclk.graph.row_map;
736 auto ind = A_lclk.graph.entries;
737 auto val = A_lclk.values;
739 const size_t numVecs = std::min (X.getNumVectors (), Y.getNumVectors ());
741 for (
size_t j = 0; j < numVecs; ++j) {
742 auto X_j = X.getVectorNonConst (j);
743 auto Y_j = Y.getVector (j);
744 auto X_lcl = X_j->getLocalViewDevice (Tpetra::Access::ReadWrite);
745 auto Y_lcl = Y_j->getLocalViewDevice (Tpetra::Access::ReadOnly);
746 auto X_lcl_1d = Kokkos::subview (X_lcl, Kokkos::ALL (), 0);
747 auto Y_lcl_1d = Kokkos::subview (Y_lcl, Kokkos::ALL (), 0);
748 KokkosSparse::Experimental::sptrsv_solve(kh_.getRawPtr(), ptr, ind, val, Y_lcl_1d, X_lcl_1d);
750 typename k_handle::HandleExecSpace().fence();
755 const std::string diag = this->diag_;
760 auto A_lcl = this->A_crs_->getLocalMatrixHost ();
762 if (X.isConstantStride () && Y.isConstantStride ()) {
763 auto X_lcl = X.getLocalViewHost (Tpetra::Access::ReadWrite);
764 auto Y_lcl = Y.getLocalViewHost (Tpetra::Access::ReadOnly);
765 KokkosSparse::trsv (uplo.c_str (), trans.c_str (), diag.c_str (),
766 A_lcl, Y_lcl, X_lcl);
769 const size_t numVecs =
770 std::min (X.getNumVectors (), Y.getNumVectors ());
771 for (
size_t j = 0; j < numVecs; ++j) {
772 auto X_j = X.getVectorNonConst (j);
773 auto Y_j = Y.getVector (j);
774 auto X_lcl = X_j->getLocalViewHost (Tpetra::Access::ReadWrite);
775 auto Y_lcl = Y_j->getLocalViewHost (Tpetra::Access::ReadOnly);
776 KokkosSparse::trsv (uplo.c_str (), trans.c_str (),
777 diag.c_str (), A_lcl, Y_lcl, X_lcl);
783 template<
class MatrixType>
785 LocalSparseTriangularSolver<MatrixType>::
786 localApply (
const MV& X,
788 const Teuchos::ETransp mode,
789 const scalar_type& alpha,
790 const scalar_type& beta)
const 792 if (mode == Teuchos::NO_TRANS && Teuchos::nonnull (htsImpl_) &&
793 htsImpl_->isComputed ()) {
794 htsImpl_->localApply (X, Y, mode, alpha, beta);
799 typedef scalar_type ST;
800 typedef Teuchos::ScalarTraits<ST> STS;
802 if (beta == STS::zero ()) {
803 if (alpha == STS::zero ()) {
804 Y.putScalar (STS::zero ());
807 this->localTriangularSolve (X, Y, mode);
808 if (alpha != STS::one ()) {
814 if (alpha == STS::zero ()) {
818 MV Y_tmp (Y, Teuchos::Copy);
819 this->localTriangularSolve (X, Y_tmp, mode);
820 Y.update (alpha, Y_tmp, beta);
826 template <
class MatrixType>
830 return numInitialize_;
833 template <
class MatrixType>
840 template <
class MatrixType>
847 template <
class MatrixType>
851 return initializeTime_;
854 template<
class MatrixType>
861 template<
class MatrixType>
868 template <
class MatrixType>
873 std::ostringstream os;
878 os <<
"\"Ifpack2::LocalSparseTriangularSolver\": {";
879 if (this->getObjectLabel () !=
"") {
880 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
882 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", " 883 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
886 os <<
"Matrix: null";
889 os <<
"Matrix: not null" 890 <<
", Global matrix dimensions: [" 891 << A_->getGlobalNumRows () <<
", " 892 << A_->getGlobalNumCols () <<
"]";
895 if (Teuchos::nonnull (htsImpl_))
896 os <<
", HTS computed: " << (htsImpl_->isComputed () ?
"true" :
"false");
902 template <
class MatrixType>
905 const Teuchos::EVerbosityLevel verbLevel)
const 910 const Teuchos::EVerbosityLevel vl
911 = (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
913 if (vl != Teuchos::VERB_NONE) {
918 auto comm = A_.is_null () ?
919 Tpetra::getDefaultComm () :
924 if (! comm.is_null () && comm->getRank () == 0) {
926 Teuchos::OSTab tab0 (out);
929 out <<
"\"Ifpack2::LocalSparseTriangularSolver\":" << endl;
930 Teuchos::OSTab tab1 (out);
931 out <<
"Scalar: " << Teuchos::TypeNameTraits<scalar_type>::name () << endl
932 <<
"LocalOrdinal: " << Teuchos::TypeNameTraits<local_ordinal_type>::name () << endl
933 <<
"GlobalOrdinal: " << Teuchos::TypeNameTraits<global_ordinal_type>::name () << endl
934 <<
"Node: " << Teuchos::TypeNameTraits<node_type>::name () << endl;
939 template <
class MatrixType>
940 Teuchos::RCP<const typename LocalSparseTriangularSolver<MatrixType>::map_type>
944 TEUCHOS_TEST_FOR_EXCEPTION
945 (A_.is_null (), std::runtime_error,
946 "Ifpack2::LocalSparseTriangularSolver::getDomainMap: " 947 "The matrix is null. Please call setMatrix() with a nonnull input " 948 "before calling this method.");
949 return A_->getDomainMap ();
952 template <
class MatrixType>
953 Teuchos::RCP<const typename LocalSparseTriangularSolver<MatrixType>::map_type>
957 TEUCHOS_TEST_FOR_EXCEPTION
958 (A_.is_null (), std::runtime_error,
959 "Ifpack2::LocalSparseTriangularSolver::getRangeMap: " 960 "The matrix is null. Please call setMatrix() with a nonnull input " 961 "before calling this method.");
962 return A_->getRangeMap ();
965 template<
class MatrixType>
967 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
969 const char prefix[] =
"Ifpack2::LocalSparseTriangularSolver::setMatrix: ";
975 if (A.getRawPtr () != A_.getRawPtr () || isInternallyChanged_) {
977 TEUCHOS_TEST_FOR_EXCEPTION
978 (! A.is_null () && A->getComm ()->getSize () == 1 &&
979 A->getNodeNumRows () != A->getNodeNumCols (),
980 std::runtime_error, prefix <<
"If A's communicator only contains one " 981 "process, then A must be square. Instead, you provided a matrix A with " 982 << A->getNodeNumRows () <<
" rows and " << A->getNodeNumCols ()
988 isInitialized_ =
false;
994 A_crs_ = Teuchos::null;
998 Teuchos::RCP<const crs_matrix_type> A_crs =
1000 TEUCHOS_TEST_FOR_EXCEPTION
1001 (A_crs.is_null (), std::invalid_argument, prefix <<
1002 "The input matrix A is not a Tpetra::CrsMatrix.");
1007 if (Teuchos::nonnull (htsImpl_))
1014 #define IFPACK2_LOCALSPARSETRIANGULARSOLVER_INSTANT(S,LO,GO,N) \ 1015 template class Ifpack2::LocalSparseTriangularSolver< Tpetra::RowMatrix<S, LO, GO, N> >; 1017 #endif // IFPACK2_LOCALSPARSETRIANGULARSOLVER_DEF_HPP 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:596
int getNumInitialize() const
Return the number of calls to initialize().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:829
Teuchos::RCP< const map_type > getDomainMap() const
The domain of this operator.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:942
Teuchos::RCP< const map_type > getRangeMap() const
The range of this operator.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:955
int getNumCompute() const
Return the number of calls to compute().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:836
double getApplyTime() const
Return the time spent in apply().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:864
void compute()
"Numeric" phase of setup
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:533
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:904
MatrixType::global_ordinal_type global_ordinal_type
Type of the global indices of the input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:95
Ifpack2 implementation details.
MatrixType::node_type node_type
Node type of the input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:97
bool isComputed() const
Return true if compute() has been called.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:222
double getComputeTime() const
Return the time spent in compute().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:857
int getNumApply() const
Return the number of calls to apply().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:843
"Preconditioner" that solves local sparse triangular systems.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:79
void initialize()
"Symbolic" phase of setup
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:375
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
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set this preconditioner's matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:967
double getInitializeTime() const
Return the time spent in initialize().
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:850
MatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:93
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
void setParameters(const Teuchos::ParameterList ¶ms)
Set this object's parameters.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:331
MatrixType::scalar_type scalar_type
Type of the entries of the input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:91
std::string description() const
A one-line description of this object.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:871
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
virtual ~LocalSparseTriangularSolver()
Destructor (virtual for memory safety).
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:320
LocalSparseTriangularSolver()
Constructor that takes no input matrix.
Definition: Ifpack2_LocalSparseTriangularSolver_def.hpp:283