41 #ifndef IFPACK2_MDF_DEF_HPP
42 #define IFPACK2_MDF_DEF_HPP
44 #include "Ifpack2_LocalFilter.hpp"
46 #include "Tpetra_CrsMatrix.hpp"
47 #include "Teuchos_StandardParameterEntryValidators.hpp"
48 #include "Ifpack2_LocalSparseTriangularSolver.hpp"
49 #include "Ifpack2_Details_getParamTryingTypes.hpp"
50 #include "Kokkos_Core.hpp"
51 #include "Kokkos_Sort.hpp"
52 #include "KokkosKernels_Sorting.hpp"
54 #include <type_traits>
63 template<
class dev_view_t>
64 auto copy_view(
const dev_view_t & vals)
66 using Kokkos::view_alloc;
67 using Kokkos::WithoutInitializing;
68 typename dev_view_t::non_const_type newvals (view_alloc (vals.label(), WithoutInitializing), vals.extent (0));
69 Kokkos::deep_copy(newvals,vals);
73 template<
class array_t,
class dev_view_t>
74 void copy_dev_view_to_host_array(array_t & array,
const dev_view_t & dev_view)
76 using host_view_t =
typename dev_view_t::HostMirror;
79 const auto ext = dev_view.extent(0);
82 ext !=
size_t(array.size()), std::logic_error,
"Ifpack2::MDF::copy_dev_view_to_host_array: "
83 "Size of permuations on host and device do not match. "
84 "Please report this bug to the Ifpack2 developers.");
87 Kokkos::deep_copy(host_view_t(array.get(),ext),dev_view);
90 template<
class scalar_type,
class local_ordinal_type,
class global_ordinal_type,
class node_type>
91 void applyReorderingPermutations(
92 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
93 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
97 "Ifpack2::MDF::applyReorderingPermuations ERROR: X.getNumVectors() != Y.getNumVectors().");
102 for(
size_t k=0; k < X.getNumVectors(); k++)
103 for(local_ordinal_type i=0; (size_t)i< X.getLocalLength(); i++)
104 y_ptr[k][perm[i]] = x_ptr[k][i];
108 template<
class scalar_type,
class local_ordinal_type,
class global_ordinal_type,
class node_type>
109 auto get_local_crs_row_matrix(
110 Teuchos::RCP<
const Tpetra::RowMatrix<scalar_type,local_ordinal_type,global_ordinal_type,node_type>> A_local)
115 using Teuchos::rcp_const_cast;
116 using Teuchos::rcp_dynamic_cast;
118 using crs_matrix_type = Tpetra::CrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type>;
120 using nonconst_local_inds_host_view_type =
typename crs_matrix_type::nonconst_local_inds_host_view_type;
121 using nonconst_values_host_view_type =
typename crs_matrix_type::nonconst_values_host_view_type;
123 RCP<const crs_matrix_type> A_local_crs = rcp_dynamic_cast<
const crs_matrix_type>(A_local);
124 if (A_local_crs.is_null ()) {
125 local_ordinal_type numRows = A_local->getLocalNumRows();
126 Array<size_t> entriesPerRow(numRows);
127 for(local_ordinal_type i = 0; i < numRows; i++) {
128 entriesPerRow[i] = A_local->getNumEntriesInLocalRow(i);
130 RCP<crs_matrix_type> A_local_crs_nc =
131 rcp (
new crs_matrix_type (A_local->getRowMap (),
132 A_local->getColMap (),
135 nonconst_local_inds_host_view_type indices(
"indices",A_local->getLocalMaxNumRowEntries());
136 nonconst_values_host_view_type values(
"values",A_local->getLocalMaxNumRowEntries());
137 for(local_ordinal_type i = 0; i < numRows; i++) {
138 size_t numEntries = 0;
139 A_local->getLocalRowCopy(i, indices, values, numEntries);
140 A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()), indices.data());
142 A_local_crs_nc->fillComplete (A_local->getDomainMap (), A_local->getRangeMap ());
143 A_local_crs = rcp_const_cast<
const crs_matrix_type> (A_local_crs_nc);
155 template<
class MatrixType>
161 isAllocated_ (false),
162 isInitialized_ (false),
167 initializeTime_ (0.0),
172 allocatePermutations();
176 template<
class MatrixType>
182 isAllocated_ (false),
183 isInitialized_ (false),
188 initializeTime_ (0.0),
193 allocatePermutations();
196 template<
class MatrixType>
202 if (force || permutations_.is_null() || A_->getLocalNumRows() != size_t(permutations_.size()))
204 permutations_ = Teuchos::null;
205 reversePermutations_ = Teuchos::null;
206 permutations_ = permutations_type(A_->getLocalNumRows());
207 reversePermutations_ = permutations_type(A_->getLocalNumRows());
211 template<
class MatrixType>
212 void MDF<MatrixType>::allocateSolvers ()
214 L_solver_ = Teuchos::null;
215 U_solver_ = Teuchos::null;
216 L_solver_ =
Teuchos::rcp (
new LocalSparseTriangularSolver<row_matrix_type> ());
217 L_solver_->setObjectLabel(
"lower");
218 U_solver_ =
Teuchos::rcp (
new LocalSparseTriangularSolver<row_matrix_type> ());
219 U_solver_->setObjectLabel(
"upper");
222 template<
class MatrixType>
231 isAllocated_ =
false;
232 isInitialized_ =
false;
234 A_local_ = Teuchos::null;
235 MDF_handle_ = Teuchos::null;
242 if (! L_solver_.is_null ()) {
243 L_solver_->setMatrix (Teuchos::null);
245 if (! U_solver_.is_null ()) {
246 U_solver_->setMatrix (Teuchos::null);
253 allocatePermutations(
true);
259 template<
class MatrixType>
264 L_.is_null (), std::runtime_error,
"Ifpack2::MDF::getL: The L factor "
265 "is null. Please call initialize() and compute() "
266 "before calling this method. If the input matrix has not yet been set, "
267 "you must first call setMatrix() with a nonnull input matrix before you "
268 "may call initialize() or compute().");
272 template<
class MatrixType>
277 permutations_.is_null (), std::runtime_error,
"Ifpack2::MDF::getPermutations: "
278 "The permulations are null. Please call initialize() and compute() "
279 "before calling this method. If the input matrix has not yet been set, "
280 "you must first call setMatrix() with a nonnull input matrix before you "
281 "may call initialize() or compute().");
284 template<
class MatrixType>
289 reversePermutations_.is_null (), std::runtime_error,
"Ifpack2::MDF::getReversePermutations: "
290 "The permulations are null. Please call initialize() and compute() "
291 "before calling this method. If the input matrix has not yet been set, "
292 "you must first call setMatrix() with a nonnull input matrix before you "
293 "may call initialize() or compute().");
297 template<
class MatrixType>
302 U_.is_null (), std::runtime_error,
"Ifpack2::MDF::getU: The U factor "
303 "is null. Please call initialize() and compute() "
304 "before calling this method. If the input matrix has not yet been set, "
305 "you must first call setMatrix() with a nonnull input matrix before you "
306 "may call initialize() or compute().");
311 template<
class MatrixType>
314 A_.
is_null (), std::runtime_error,
"Ifpack2::MDF::getNodeSmootherComplexity: "
315 "The input matrix A is null. Please call setMatrix() with a nonnull "
316 "input matrix, then call compute(), before calling this method.");
318 if(!L_.is_null() && !U_.is_null())
319 return A_->getLocalNumEntries() + L_->getLocalNumEntries() + U_->getLocalNumEntries();
325 template<
class MatrixType>
331 A_.
is_null (), std::runtime_error,
"Ifpack2::MDF::getDomainMap: "
332 "The matrix is null. Please call setMatrix() with a nonnull input "
333 "before calling this method.");
337 L_.is_null (), std::runtime_error,
"Ifpack2::MDF::getDomainMap: "
338 "The computed graph is null. Please call initialize() and compute() before calling "
340 return L_->getDomainMap ();
344 template<
class MatrixType>
350 A_.
is_null (), std::runtime_error,
"Ifpack2::MDF::getRangeMap: "
351 "The matrix is null. Please call setMatrix() with a nonnull input "
352 "before calling this method.");
356 L_.is_null (), std::runtime_error,
"Ifpack2::MDF::getRangeMap: "
357 "The computed graph is null. Please call initialize() abd compute() before calling "
359 return L_->getRangeMap ();
362 template<
class MatrixType>
370 using Details::getParamTryingTypes;
371 const char prefix[] =
"Ifpack2::MDF: ";
375 double overalloc = 2.;
387 const std::string paramName (
"fact: mdf level-of-fill");
388 getParamTryingTypes<int, int, global_ordinal_type, double, float>
389 (fillLevel, params, paramName, prefix);
392 (fillLevel != 0, std::runtime_error, prefix <<
"MDF with level of fill != 0 is not yet implemented.");
395 const std::string paramName (
"Verbosity");
396 getParamTryingTypes<int, int, global_ordinal_type, double, float>
397 (verbosity, params, paramName, prefix);
400 const std::string paramName (
"fact: mdf overalloc");
401 getParamTryingTypes<double, double>
402 (overalloc, params, paramName, prefix);
406 L_solver_->setParameters(params);
407 U_solver_->setParameters(params);
413 LevelOfFill_ = fillLevel;
414 Overalloc_ = overalloc;
415 Verbosity_ = verbosity;
419 template<
class MatrixType>
426 template<
class MatrixType>
433 template<
class MatrixType>
439 using Teuchos::rcp_dynamic_cast;
440 using Teuchos::rcp_implicit_cast;
445 if (A->getRowMap ()->getComm ()->getSize () == 1 ||
446 A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
453 RCP<const LocalFilter<row_matrix_type> > A_lf_r =
454 rcp_dynamic_cast<
const LocalFilter<row_matrix_type> > (A);
455 if (! A_lf_r.is_null ()) {
456 return rcp_implicit_cast<
const row_matrix_type> (A_lf_r);
462 return rcp (
new LocalFilter<row_matrix_type> (A));
467 template<
class MatrixType>
472 using Teuchos::rcp_const_cast;
473 using Teuchos::rcp_dynamic_cast;
474 using Teuchos::rcp_implicit_cast;
477 const char prefix[] =
"Ifpack2::MDF::initialize: ";
480 (A_.
is_null (), std::runtime_error, prefix <<
"The matrix is null. Please "
481 "call setMatrix() with a nonnull input before calling this method.");
483 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not "
484 "fill complete. You may not invoke initialize() or compute() with this "
485 "matrix until the matrix is fill complete. If your matrix is a "
486 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
487 "range Maps, if appropriate) before calling this method.");
490 double startTime = timer.wallTime();
501 isInitialized_ =
false;
502 isAllocated_ =
false;
504 MDF_handle_ = Teuchos::null;
506 A_local_ = makeLocalFilter (A_);
508 A_local_.is_null (), std::logic_error,
"Ifpack2::MDF::initialize: "
509 "makeLocalFilter returned null; it failed to compute A_local. "
510 "Please report this bug to the Ifpack2 developers.");
518 RCP<const crs_matrix_type> A_local_crs = Details::MDFImpl::get_local_crs_row_matrix(A_local_);
520 auto A_local_device = A_local_crs->getLocalMatrixDevice();
521 MDF_handle_ =
rcp(
new MDF_handle_device_type(A_local_device) );
522 MDF_handle_->set_verbosity(Verbosity_);
524 KokkosSparse::Experimental::mdf_symbolic(A_local_device,*MDF_handle_);
529 checkOrderingConsistency (*A_local_);
532 isInitialized_ =
true;
534 initializeTime_ += (timer.wallTime() - startTime);
537 template<
class MatrixType>
547 bool gidsAreConsistentlyOrdered=
true;
548 global_ordinal_type indexOfInconsistentGID=0;
549 for (global_ordinal_type i=0; i<rowGIDs.
size(); ++i) {
550 if (rowGIDs[i] != colGIDs[i]) {
551 gidsAreConsistentlyOrdered=
false;
552 indexOfInconsistentGID=i;
557 "The ordering of the local GIDs in the row and column maps is not the same"
558 << std::endl <<
"at index " << indexOfInconsistentGID
559 <<
". Consistency is required, as all calculations are done with"
560 << std::endl <<
"local indexing.");
563 template<
class MatrixType>
568 using Teuchos::rcp_const_cast;
569 using Teuchos::rcp_dynamic_cast;
572 const char prefix[] =
"Ifpack2::MDF::compute: ";
578 (A_.
is_null (), std::runtime_error, prefix <<
"The matrix is null. Please "
579 "call setMatrix() with a nonnull input before calling this method.");
581 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not "
582 "fill complete. You may not invoke initialize() or compute() with this "
583 "matrix until the matrix is fill complete. If your matrix is a "
584 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
585 "range Maps, if appropriate) before calling this method.");
587 if (! isInitialized ()) {
595 double startTime = timer.
wallTime();
600 RCP<const crs_matrix_type> A_local_crs = Details::MDFImpl::get_local_crs_row_matrix(A_local_);
603 auto A_local_device = A_local_crs->getLocalMatrixDevice();
605 KokkosSparse::Experimental::mdf_numeric(A_local_device,*MDF_handle_);
609 Details::MDFImpl::copy_dev_view_to_host_array(reversePermutations_, MDF_handle_->permutation);
610 Details::MDFImpl::copy_dev_view_to_host_array(permutations_, MDF_handle_->permutation_inv);
615 auto L_mdf = MDF_handle_->getL();
617 A_local_->getRowMap (),
618 A_local_->getColMap (),
619 Details::MDFImpl::copy_view(L_mdf.graph.row_map),
620 Details::MDFImpl::copy_view(L_mdf.graph.entries),
621 Details::MDFImpl::copy_view(L_mdf.values)
625 auto U_mdf = MDF_handle_->getU();
627 A_local_->getRowMap (),
628 A_local_->getColMap (),
629 Details::MDFImpl::copy_view(U_mdf.graph.row_map),
630 Details::MDFImpl::copy_view(U_mdf.graph.entries),
631 Details::MDFImpl::copy_view(U_mdf.values)
636 L_solver_->setMatrix (L_);
637 L_solver_->initialize ();
638 L_solver_->compute ();
639 U_solver_->setMatrix (U_);
640 U_solver_->initialize ();
641 U_solver_->compute ();
645 computeTime_ += (timer.
wallTime() - startTime);
648 template<
class MatrixType>
651 apply_impl (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
652 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
655 scalar_type beta)
const
657 const scalar_type one = STS::one ();
658 const scalar_type zero = STS::zero ();
660 if (alpha == one && beta == zero) {
661 MV tmp (Y.getMap (), Y.getNumVectors ());
662 Details::MDFImpl::applyReorderingPermutations(X,tmp,permutations_);
665 L_solver_->apply (tmp, Y, mode);
666 U_solver_->apply (Y, tmp, mode);
670 U_solver_->apply (tmp, Y, mode);
671 L_solver_->apply (Y, tmp, mode);
673 Details::MDFImpl::applyReorderingPermutations(tmp,Y,reversePermutations_);
686 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
687 apply_impl (X, Y_tmp, mode);
688 Y.update (alpha, Y_tmp, beta);
693 template<
class MatrixType>
696 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
697 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
703 using Teuchos::rcpFromRef;
706 A_.
is_null (), std::runtime_error,
"Ifpack2::MDF::apply: The matrix is "
707 "null. Please call setMatrix() with a nonnull input, then initialize() "
708 "and compute(), before calling this method.");
710 ! isComputed (), std::runtime_error,
711 "Ifpack2::MDF::apply: If you have not yet called compute(), "
712 "you must call compute() before calling this method.");
713 TEUCHOS_TEST_FOR_EXCEPTION(
714 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
715 "Ifpack2::MDF::apply: X and Y do not have the same number of columns. "
716 "X.getNumVectors() = " << X.getNumVectors ()
717 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
718 TEUCHOS_TEST_FOR_EXCEPTION(
720 "Ifpack2::MDF::apply: mode = Teuchos::CONJ_TRANS is not implemented for "
721 "complex Scalar type. Please talk to the Ifpack2 developers to get this "
722 "fixed. There is a FIXME in this file about this very issue.");
723 #ifdef HAVE_IFPACK2_DEBUG
728 for (
size_t j = 0; j < X.getNumVectors (); ++j) {
729 if (STM::isnaninf (norms[j])) {
734 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::MDF::apply: The 1-norm of the input X is NaN or Inf.");
736 #endif // HAVE_IFPACK2_DEBUG
739 double startTime = timer.
wallTime();
742 apply_impl(X,Y,mode,alpha,beta);
745 #ifdef HAVE_IFPACK2_DEBUG
750 for (
size_t j = 0; j < Y.getNumVectors (); ++j) {
751 if (STM::isnaninf (norms[j])) {
756 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::MDF::apply: The 1-norm of the output Y is NaN or Inf.");
758 #endif // HAVE_IFPACK2_DEBUG
761 applyTime_ += (timer.
wallTime() - startTime);
764 template<
class MatrixType>
767 std::ostringstream os;
772 os <<
"\"Ifpack2::MDF\": {";
773 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", "
774 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
776 os <<
"Level-of-fill: " << getLevelOfFill() <<
", ";
779 os <<
"Matrix: null";
782 os <<
"Global matrix dimensions: ["
783 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]"
784 <<
", Global nnz: " << A_->getGlobalNumEntries();
787 if (! L_solver_.is_null ()) os <<
", " << L_solver_->description ();
788 if (! U_solver_.is_null ()) os <<
", " << U_solver_->description ();
796 #define IFPACK2_MDF_INSTANT(S,LO,GO,N) \
797 template class Ifpack2::MDF< Tpetra::RowMatrix<S, LO, GO, N> >;
Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > crs_matrix_type
Tpetra::CrsMatrix specialization used by this class for representing L and U.
Definition: Ifpack2_MDF_decl.hpp:122
Ifpack2::ScalingType enumerable type.
void compute()
Compute the (numeric) incomplete factorization.
Definition: Ifpack2_MDF_def.hpp:564
MDF (incomplete LU factorization with minimum discarded fill reordering) of a Tpetra sparse matrix...
Definition: Ifpack2_MDF_decl.hpp:83
permutations_type & getReversePermutations() const
Return the reverse permutations of the MDF factorization.
Definition: Ifpack2_MDF_def.hpp:286
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_MDF_decl.hpp:104
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::string description() const
A one-line description of this object.
Definition: Ifpack2_MDF_def.hpp:765
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_MDF_def.hpp:348
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_MDF_decl.hpp:101
Teuchos::RCP< const crs_matrix_type > getCrsMatrix() const
Return the input matrix A as a Tpetra::CrsMatrix, if possible; else throws.
Definition: Ifpack2_MDF_def.hpp:428
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_MDF_def.hpp:312
permutations_type & getPermutations() const
Return the permutations of the MDF factorization.
Definition: Ifpack2_MDF_def.hpp:274
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_MDF_def.hpp:224
void setParameters(const Teuchos::ParameterList ¶ms)
Definition: Ifpack2_MDF_def.hpp:365
const crs_matrix_type & getL() const
Return the L factor of the MDF factorization.
Definition: Ifpack2_MDF_def.hpp:261
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_MDF_def.hpp:329
Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > row_matrix_type
Tpetra::RowMatrix specialization used by this class.
Definition: Ifpack2_MDF_decl.hpp:119
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_MDF_decl.hpp:95
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the input matrix.
Definition: Ifpack2_MDF_def.hpp:421
const crs_matrix_type & getU() const
Return the U factor of the MDF factorization.
Definition: Ifpack2_MDF_def.hpp:299
void initialize()
Initialize by computing the symbolic incomplete factorization.
Definition: Ifpack2_MDF_def.hpp:468
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 (inverse of the) incomplete factorization to X, resulting in Y.
Definition: Ifpack2_MDF_def.hpp:696