Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_LocalSparseTriangularSolver_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK2_LOCALSPARSETRIANGULARSOLVER_DEF_HPP
44 #define IFPACK2_LOCALSPARSETRIANGULARSOLVER_DEF_HPP
45 
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"
51 
52 #ifdef HAVE_IFPACK2_SHYLU_NODEHTS
53 # include "shylu_hts.hpp"
54 #endif
55 
56 namespace Ifpack2 {
57 
58 namespace Details {
59 struct TrisolverType {
60  enum Enum {
61  Internal,
62  HTS,
63  KSPTRSV
64  };
65 
66  static void loadPLTypeOption (Teuchos::Array<std::string>& type_strs, Teuchos::Array<Enum>& type_enums) {
67  type_strs.resize(3);
68  type_strs[0] = "Internal";
69  type_strs[1] = "HTS";
70  type_strs[2] = "KSPTRSV";
71  type_enums.resize(3);
72  type_enums[0] = Internal;
73  type_enums[1] = HTS;
74  type_enums[2] = KSPTRSV;
75  }
76 };
77 }
78 
79 template<class MatrixType>
80 class LocalSparseTriangularSolver<MatrixType>::HtsImpl {
81 public:
82  typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> crs_matrix_type;
83 
84  void reset () {
85 #ifdef HAVE_IFPACK2_SHYLU_NODEHTS
86  Timpl_ = Teuchos::null;
87  levelset_block_size_ = 1;
88 #endif
89  }
90 
91  void setParameters (const Teuchos::ParameterList& pl) {
92  (void)pl;
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);
99  }
100  if (levelset_block_size_ < 1)
101  levelset_block_size_ = 1;
102 #endif
103  }
104 
105  // HTS has the phases symbolic+numeric, numeric, and apply. Hence the first
106  // call to compute() will trigger the symbolic+numeric phase, and subsequent
107  // calls (with the same Timpl_) will trigger the numeric phase. In the call to
108  // initialize(), essentially nothing happens.
109  void initialize (const crs_matrix_type& /* unused */) {
110 #ifdef HAVE_IFPACK2_SHYLU_NODEHTS
111  reset();
112  transpose_ = conjugate_ = false;
113 #endif
114  }
115 
116  void compute (const crs_matrix_type& T_in, const Teuchos::RCP<Teuchos::FancyOStream>& out) {
117  (void)T_in;
118  (void)out;
119 #ifdef HAVE_IFPACK2_SHYLU_NODEHTS
120  using Teuchos::ArrayRCP;
121 
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);
126  Kokkos::fence();
127 
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());
133 
134  if (Teuchos::nonnull(Timpl_)) {
135  // Reuse the nonzero pattern.
136  HTST::reprocess_numeric(Timpl_.get(), T_hts.get());
137  } else {
138  // Build from scratch.
139  if (T_in.getCrsGraph().is_null()) {
140  if (Teuchos::nonnull(out))
141  *out << "HTS compute failed because T_in.getCrsGraph().is_null().\n";
142  return;
143  }
144  if ( ! T_in.getCrsGraph()->isSorted()) {
145  if (Teuchos::nonnull(out))
146  *out << "HTS compute failed because ! T_in.getCrsGraph().isSorted().\n";
147  return;
148  }
149  if ( ! T_in.isStorageOptimized()) {
150  if (Teuchos::nonnull(out))
151  *out << "HTS compute failed because ! T_in.isStorageOptimized().\n";
152  return;
153  }
154 
155  typename HTST::PreprocessArgs args;
156  args.T = T_hts.get();
157  args.max_nrhs = 1;
158 #ifdef _OPENMP
159  args.nthreads = omp_get_max_threads();
160 #else
161  args.nthreads = 1;
162 #endif
163  args.save_for_reprocess = true;
164  typename HTST::Options opts;
165  opts.levelset_block_size = levelset_block_size_;
166  args.options = &opts;
167 
168  try {
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";
173  }
174  }
175 #endif
176  }
177 
178  // HTS may not be able to handle a matrix, so query whether compute()
179  // succeeded.
180  bool isComputed () {
181 #ifdef HAVE_IFPACK2_SHYLU_NODEHTS
182  return Teuchos::nonnull(Timpl_);
183 #else
184  return false;
185 #endif
186  }
187 
188  // Y := beta * Y + alpha * (M * X)
189  void localApply (const MV& X, MV& Y,
190  const Teuchos::ETransp /* mode */,
191  const scalar_type& alpha, const scalar_type& beta) const {
192  (void)X;
193  (void)Y;
194  (void)alpha;
195  (void)beta;
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);
199 
200  // Only does something if #rhs > current capacity.
201  HTST::reset_max_nrhs(Timpl_.get(), X_view.extent(1));
202  // Switch alpha and beta because of HTS's opposite convention.
203  HTST::solve_omp(Timpl_.get(),
204  // For std/Kokkos::complex.
205  reinterpret_cast<const scalar_type*>(X_view.data()),
206  X_view.extent(1),
207  // For std/Kokkos::complex.
208  reinterpret_cast<scalar_type*>(Y_view.data()),
209  beta, alpha);
210 #endif
211  }
212 
213 private:
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;
218 
219  struct TImplDeleter {
220  void free (TImpl* impl) {
221  HTST::delete_Impl(impl);
222  }
223  };
224 
225  struct HtsCrsMatrixDeleter {
226  void free (HtsCrsMatrix* T) {
227  HTST::delete_CrsMatrix(T);
228  }
229  };
230 
231  Teuchos::RCP<TImpl> Timpl_;
232  bool transpose_, conjugate_;
233  int levelset_block_size_;
234 #endif
235 };
236 
237 template<class MatrixType>
239 LocalSparseTriangularSolver (const Teuchos::RCP<const row_matrix_type>& A) :
240  A_ (A)
241 {
242  initializeState();
243  typedef typename Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
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.");
252  A_crs_ = A_crs;
253  }
254 }
255 
256 template<class MatrixType>
258 LocalSparseTriangularSolver (const Teuchos::RCP<const row_matrix_type>& A,
259  const Teuchos::RCP<Teuchos::FancyOStream>& out) :
260  A_ (A),
261  out_ (out)
262 {
263  initializeState();
264  if (! out_.is_null ()) {
265  *out_ << ">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor"
266  << std::endl;
267  }
268  typedef typename Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
270  if (! A.is_null ()) {
271  Teuchos::RCP<const crs_matrix_type> A_crs =
272  Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A);
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.");
277  A_crs_ = A_crs;
278  }
279 }
280 
281 template<class MatrixType>
284 {
285  initializeState();
286 }
287 
288 template<class MatrixType>
290 LocalSparseTriangularSolver (const bool /* unused */, const Teuchos::RCP<Teuchos::FancyOStream>& out) :
291  out_ (out)
292 {
293  initializeState();
294  if (! out_.is_null ()) {
295  *out_ << ">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor"
296  << std::endl;
297  }
298 }
299 
300 template<class MatrixType>
302 {
303  isInitialized_ = false;
304  isComputed_ = false;
305  reverseStorage_ = false;
306  isInternallyChanged_ = false;
307  numInitialize_ = 0;
308  numCompute_ = 0;
309  numApply_ = 0;
310  initializeTime_ = 0.0;
311  computeTime_ = 0.0;
312  applyTime_ = 0.0;
313  isKokkosKernelsSptrsv_ = false;
314  uplo_ = "N";
315  diag_ = "N";
316 }
317 
318 template<class MatrixType>
321 {
322  if (Teuchos::nonnull (kh_))
323  {
324  kh_->destroy_sptrsv_handle();
325  }
326 }
327 
328 template<class MatrixType>
329 void
331 setParameters (const Teuchos::ParameterList& pl)
332 {
333  using Teuchos::RCP;
334  using Teuchos::ParameterList;
335  using Teuchos::Array;
336 
337  Details::TrisolverType::Enum trisolverType = Details::TrisolverType::Internal;
338  do {
339  static const char typeName[] = "trisolver: type";
340 
341  if ( ! pl.isType<std::string>(typeName)) break;
342 
343  // Map std::string <-> TrisolverType::Enum.
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);
349 
350  trisolverType = s2i.getIntegralValue(pl.get<std::string>(typeName));
351  } while (0);
352 
353  if (trisolverType == Details::TrisolverType::HTS) {
354  htsImpl_ = Teuchos::rcp (new HtsImpl ());
355  htsImpl_->setParameters (pl);
356  }
357 
358  if (trisolverType == Details::TrisolverType::KSPTRSV) {
359  kh_ = Teuchos::rcp (new k_handle());
360  }
361 
362  if (pl.isParameter("trisolver: reverse U"))
363  reverseStorage_ = pl.get<bool>("trisolver: reverse U");
364 
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.");
370 }
371 
372 template<class MatrixType>
373 void
376 {
377  using Tpetra::Details::determineLocalTriangularStructure;
378  using crs_matrix_type = Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
380  using local_matrix_type = typename crs_matrix_type::local_matrix_device_type;
381  using LO = local_ordinal_type;
382 
383  const char prefix[] = "Ifpack2::LocalSparseTriangularSolver::initialize: ";
384  if (! out_.is_null ()) {
385  *out_ << ">>> DEBUG " << prefix << std::endl;
386  }
387 
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 ()) {
393  auto A_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
394  TEUCHOS_TEST_FOR_EXCEPTION
395  (A_crs.get () == nullptr, std::invalid_argument,
396  prefix << "The input matrix A is not a Tpetra::CrsMatrix.");
397  A_crs_ = A_crs;
398  }
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.");
404  // At this point, the graph MUST be fillComplete. The "initialize"
405  // (symbolic) part of setup only depends on the graph structure, so
406  // the matrix itself need not be fill complete.
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.");
410 
411  // mfh 30 Apr 2018: See GitHub Issue #2658.
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 ();
417  auto lclTriStruct =
418  determineLocalTriangularStructure (lclMatrix.graph,
419  lclRowMap,
420  lclColMap,
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");
426  return lclTriStruct;
427  } ();
428 
429  if (reverseStorage_ && lclTriStructure.couldBeUpperTriangular &&
430  htsImpl_.is_null ()) {
431  // Reverse the storage for an upper triangular matrix
432  auto Alocal = A_crs_->getLocalMatrixDevice();
433  auto ptr = Alocal.graph.row_map;
434  auto ind = Alocal.graph.entries;
435  auto val = Alocal.values;
436 
437  auto numRows = Alocal.numRows();
438  auto numCols = Alocal.numCols();
439  auto numNnz = Alocal.nnz();
440 
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));
444 
445  // FIXME: The code below assumes UVM
446  typename crs_matrix_type::execution_space().fence();
447  newptr(0) = 0;
448  for (local_ordinal_type row = 0, rowStart = 0; row < numRows; ++row) {
449  auto A_r = Alocal.row(numRows-1 - row);
450 
451  auto numEnt = A_r.length;
452  for (local_ordinal_type k = 0; k < numEnt; ++k) {
453  newind(rowStart + k) = numCols-1 - A_r.colidx(numEnt-1 - k);
454  newval(rowStart + k) = A_r.value (numEnt-1 - k);
455  }
456  rowStart += numEnt;
457  newptr(row+1) = rowStart;
458  }
459  typename crs_matrix_type::execution_space().fence();
460 
461  // Reverse maps
462  using map_type = typename crs_matrix_type::map_type;
463  Teuchos::RCP<map_type> newRowMap, newColMap;
464  {
465  // Reverse row map
466  auto rowMap = A_->getRowMap();
467  auto numElems = rowMap->getNodeNumElements();
468  auto rowElems = rowMap->getNodeElementList();
469 
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];
473 
474  newRowMap = Teuchos::rcp(new map_type(rowMap->getGlobalNumElements(), newRowElems, rowMap->getIndexBase(), rowMap->getComm()));
475  }
476  {
477  // Reverse column map
478  auto colMap = A_->getColMap();
479  auto numElems = colMap->getNodeNumElements();
480  auto colElems = colMap->getNodeElementList();
481 
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];
485 
486  newColMap = Teuchos::rcp(new map_type(colMap->getGlobalNumElements(), newColElems, colMap->getIndexBase(), colMap->getComm()));
487  }
488 
489  // Construct new matrix
490  local_matrix_type newLocalMatrix("Upermuted", numRows, numCols, numNnz, newval, newptr, newind);
491 
492  A_crs_ = Teuchos::rcp(new crs_matrix_type(newLocalMatrix, newRowMap, newColMap, A_crs_->getDomainMap(), A_crs_->getRangeMap()));
493 
494  isInternallyChanged_ = true;
495 
496  // FIXME (mfh 18 Apr 2019) Recomputing this is unnecessary, but I
497  // didn't want to break any invariants, especially considering
498  // that this branch is likely poorly tested.
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");
508  }
509 
510  if (Teuchos::nonnull (htsImpl_))
511  {
512  htsImpl_->initialize (*A_crs_);
513  isInternallyChanged_ = true;
514  }
515 
516  const bool ksptrsv_valid_uplo = (this->uplo_ != "N");
517  if (Teuchos::nonnull(kh_) && ksptrsv_valid_uplo && this->diag_ != "U")
518  {
519  this->isKokkosKernelsSptrsv_ = true;
520  }
521  else
522  {
523  this->isKokkosKernelsSptrsv_ = false;
524  }
525 
526  isInitialized_ = true;
527  ++numInitialize_;
528 }
529 
530 template<class MatrixType>
531 void
534 {
535  const char prefix[] = "Ifpack2::LocalSparseTriangularSolver::compute: ";
536  if (! out_.is_null ()) {
537  *out_ << ">>> DEBUG " << prefix << std::endl;
538  }
539 
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.");
547  // At this point, the matrix MUST be fillComplete.
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.");
551 
552  if (! isInitialized_) {
553  initialize ();
554  }
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.");
559 
560  if (Teuchos::nonnull (htsImpl_))
561  htsImpl_->compute (*A_crs_, out_);
562 
563  if (Teuchos::nonnull(kh_) && this->isKokkosKernelsSptrsv_)
564  {
565  auto A_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
566  auto Alocal = A_crs->getLocalMatrixDevice();
567  auto ptr = Alocal.graph.row_map;
568  auto ind = Alocal.graph.entries;
569  auto val = Alocal.values;
570 
571  auto numRows = Alocal.numRows();
572  const bool is_lower_tri = (this->uplo_ == "L") ? true : false;
573 
574  // Destroy existing handle and recreate in case new matrix provided - requires rerunning symbolic analysis
575  kh_->destroy_sptrsv_handle();
576 #if defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) && defined(KOKKOS_ENABLE_CUDA)
577  // CuSparse only supports int type ordinals
578  if (std::is_same<Kokkos::Cuda, HandleExecSpace>::value && std::is_same<int,local_ordinal_type >::value)
579  {
580  kh_->create_sptrsv_handle(KokkosSparse::Experimental::SPTRSVAlgorithm::SPTRSV_CUSPARSE, numRows, is_lower_tri);
581  }
582  else
583 #endif
584  {
585  kh_->create_sptrsv_handle(KokkosSparse::Experimental::SPTRSVAlgorithm::SEQLVLSCHD_TP1, numRows, is_lower_tri);
586  }
587  KokkosSparse::Experimental::sptrsv_symbolic(kh_.getRawPtr(), ptr, ind, val);
588  }
589 
590  isComputed_ = true;
591  ++numCompute_;
592 }
593 
594 template<class MatrixType>
596 apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type,
598  Tpetra::MultiVector<scalar_type, local_ordinal_type,
600  Teuchos::ETransp mode,
601  scalar_type alpha,
602  scalar_type beta) const
603 {
604  using Teuchos::RCP;
605  using Teuchos::rcp;
606  using Teuchos::rcpFromRef;
607  typedef scalar_type ST;
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;
614  }
615  else {
616  Teuchos::RCP<const crs_matrix_type> A_crs =
617  Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
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;
625  }
626  }
627 
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.");
632  // If isComputed() is true, it's impossible for the matrix to be
633  // null, or for it not to be a Tpetra::CrsMatrix.
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.");
640  // However, it _is_ possible that the user called resumeFill() on
641  // the matrix, after calling compute(). This is NOT allowed.
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.");
652 
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 ();
660 
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 ()));
664  }
665  else {
666  X_colMap_->putScalar (STS::zero ());
667  }
668  // See discussion of Github Issue #672 for why the Import needs to
669  // use the ZERO CombineMode. The case where the Export is
670  // nontrivial is likely never exercised.
671  X_colMap_->doImport (X, *importer, Tpetra::ZERO);
672  }
673  RCP<const MV> X_cur = importer.is_null () ? rcpFromRef (X) :
674  Teuchos::rcp_const_cast<const MV> (X_colMap_);
675 
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 ()));
679  }
680  else {
681  Y_rowMap_->putScalar (STS::zero ());
682  }
683  Y_rowMap_->doExport (Y, *importer, Tpetra::ADD);
684  }
685  RCP<MV> Y_cur = exporter.is_null () ? rcpFromRef (Y) : Y_rowMap_;
686 
687  localApply (*X_cur, *Y_cur, mode, alpha, beta);
688 
689  if (! exporter.is_null ()) {
690  Y.putScalar (STS::zero ());
691  Y.doExport (*Y_cur, *exporter, Tpetra::ADD);
692  }
693 
694  ++numApply_;
695 }
696 
697 template<class MatrixType>
698 void
700 localTriangularSolve (const MV& Y,
701  MV& X,
702  const Teuchos::ETransp mode) const
703 {
704  using Teuchos::CONJ_TRANS;
705  using Teuchos::NO_TRANS;
706  using Teuchos::TRANS;
707  const char tfecfFuncName[] = "localTriangularSolve: ";
708 
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.");
726 
727  const std::string uplo = this->uplo_;
728  const std::string trans = (mode == Teuchos::CONJ_TRANS) ? "C" :
729  (mode == Teuchos::TRANS ? "T" : "N");
730 
731  if (Teuchos::nonnull(kh_) && this->isKokkosKernelsSptrsv_ && trans == "N")
732  {
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;
738 
739  const size_t numVecs = std::min (X.getNumVectors (), Y.getNumVectors ());
740 
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);
749  // TODO is this fence needed...
750  typename k_handle::HandleExecSpace().fence();
751  }
752  }
753  else
754  {
755  const std::string diag = this->diag_;
756  // NOTE (mfh 20 Aug 2017): KokkosSparse::trsv currently is a
757  // sequential, host-only code. See
758  // https://github.com/kokkos/kokkos-kernels/issues/48.
759 
760  auto A_lcl = this->A_crs_->getLocalMatrixHost ();
761 
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);
767  }
768  else {
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);
778  }
779  }
780  }
781 }
782 
783 template<class MatrixType>
784 void
785 LocalSparseTriangularSolver<MatrixType>::
786 localApply (const MV& X,
787  MV& Y,
788  const Teuchos::ETransp mode,
789  const scalar_type& alpha,
790  const scalar_type& beta) const
791 {
792  if (mode == Teuchos::NO_TRANS && Teuchos::nonnull (htsImpl_) &&
793  htsImpl_->isComputed ()) {
794  htsImpl_->localApply (X, Y, mode, alpha, beta);
795  return;
796  }
797 
798  using Teuchos::RCP;
799  typedef scalar_type ST;
800  typedef Teuchos::ScalarTraits<ST> STS;
801 
802  if (beta == STS::zero ()) {
803  if (alpha == STS::zero ()) {
804  Y.putScalar (STS::zero ()); // Y := 0 * Y (ignore contents of Y)
805  }
806  else { // alpha != 0
807  this->localTriangularSolve (X, Y, mode);
808  if (alpha != STS::one ()) {
809  Y.scale (alpha);
810  }
811  }
812  }
813  else { // beta != 0
814  if (alpha == STS::zero ()) {
815  Y.scale (beta); // Y := beta * Y
816  }
817  else { // alpha != 0
818  MV Y_tmp (Y, Teuchos::Copy);
819  this->localTriangularSolve (X, Y_tmp, mode); // Y_tmp := M * X
820  Y.update (alpha, Y_tmp, beta); // Y := beta * Y + alpha * Y_tmp
821  }
822  }
823 }
824 
825 
826 template <class MatrixType>
827 int
830  return numInitialize_;
831 }
832 
833 template <class MatrixType>
834 int
836 getNumCompute () const {
837  return numCompute_;
838 }
839 
840 template <class MatrixType>
841 int
843 getNumApply () const {
844  return numApply_;
845 }
846 
847 template <class MatrixType>
848 double
851  return initializeTime_;
852 }
853 
854 template<class MatrixType>
855 double
857 getComputeTime () const {
858  return computeTime_;
859 }
860 
861 template<class MatrixType>
862 double
864 getApplyTime() const {
865  return applyTime_;
866 }
867 
868 template <class MatrixType>
869 std::string
871 description () const
872 {
873  std::ostringstream os;
874 
875  // Output is a valid YAML dictionary in flow style. If you don't
876  // like everything on a single line, you should call describe()
877  // instead.
878  os << "\"Ifpack2::LocalSparseTriangularSolver\": {";
879  if (this->getObjectLabel () != "") {
880  os << "Label: \"" << this->getObjectLabel () << "\", ";
881  }
882  os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
883  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
884 
885  if (A_.is_null ()) {
886  os << "Matrix: null";
887  }
888  else {
889  os << "Matrix: not null"
890  << ", Global matrix dimensions: ["
891  << A_->getGlobalNumRows () << ", "
892  << A_->getGlobalNumCols () << "]";
893  }
894 
895  if (Teuchos::nonnull (htsImpl_))
896  os << ", HTS computed: " << (htsImpl_->isComputed () ? "true" : "false");
897 
898  os << "}";
899  return os.str ();
900 }
901 
902 template <class MatrixType>
904 describe (Teuchos::FancyOStream& out,
905  const Teuchos::EVerbosityLevel verbLevel) const
906 {
907  using std::endl;
908  // Default verbosity level is VERB_LOW, which prints only on Process
909  // 0 of the matrix's communicator.
910  const Teuchos::EVerbosityLevel vl
911  = (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
912 
913  if (vl != Teuchos::VERB_NONE) {
914  // Print only on Process 0 in the matrix's communicator. If the
915  // matrix is null, though, we have to get the communicator from
916  // somewhere, so we ask Tpetra for its default communicator. If
917  // MPI is enabled, this wraps MPI_COMM_WORLD or a clone thereof.
918  auto comm = A_.is_null () ?
919  Tpetra::getDefaultComm () :
920  A_->getComm ();
921 
922  // Users aren't supposed to do anything with the matrix on
923  // processes where its communicator is null.
924  if (! comm.is_null () && comm->getRank () == 0) {
925  // By convention, describe() should always begin with a tab.
926  Teuchos::OSTab tab0 (out);
927  // Output is in YAML format. We have to escape the class name,
928  // because it has a colon.
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;
935  }
936  }
937 }
938 
939 template <class MatrixType>
940 Teuchos::RCP<const typename LocalSparseTriangularSolver<MatrixType>::map_type>
942 getDomainMap () const
943 {
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 ();
950 }
951 
952 template <class MatrixType>
953 Teuchos::RCP<const typename LocalSparseTriangularSolver<MatrixType>::map_type>
955 getRangeMap () const
956 {
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 ();
963 }
964 
965 template<class MatrixType>
967 setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
968 {
969  const char prefix[] = "Ifpack2::LocalSparseTriangularSolver::setMatrix: ";
970 
971  // If the pointer didn't change, do nothing. This is reasonable
972  // because users are supposed to call this method with the same
973  // object over all participating processes, and pointer identity
974  // implies object identity.
975  if (A.getRawPtr () != A_.getRawPtr () || isInternallyChanged_) {
976  // Check in serial or one-process mode if the matrix is square.
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 ()
983  << " columns.");
984 
985  // It's legal for A to be null; in that case, you may not call
986  // initialize() until calling setMatrix() with a nonnull input.
987  // Regardless, setting the matrix invalidates the preconditioner.
988  isInitialized_ = false;
989  isComputed_ = false;
990 
991  typedef typename Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
993  if (A.is_null ()) {
994  A_crs_ = Teuchos::null;
995  A_ = Teuchos::null;
996  }
997  else { // A is not null
998  Teuchos::RCP<const crs_matrix_type> A_crs =
999  Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A);
1000  TEUCHOS_TEST_FOR_EXCEPTION
1001  (A_crs.is_null (), std::invalid_argument, prefix <<
1002  "The input matrix A is not a Tpetra::CrsMatrix.");
1003  A_crs_ = A_crs;
1004  A_ = A;
1005  }
1006 
1007  if (Teuchos::nonnull (htsImpl_))
1008  htsImpl_->reset ();
1009  } // pointers are not the same
1010 }
1011 
1012 } // namespace Ifpack2
1013 
1014 #define IFPACK2_LOCALSPARSETRIANGULARSOLVER_INSTANT(S,LO,GO,N) \
1015  template class Ifpack2::LocalSparseTriangularSolver< Tpetra::RowMatrix<S, LO, GO, N> >;
1016 
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&#39;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 &params)
Set this object&#39;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