MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_DroppingCommon.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_DROPPINGCOMMON_HPP
11 #define MUELU_DROPPINGCOMMON_HPP
12 
13 #include "Kokkos_Core.hpp"
14 #include "Kokkos_ArithTraits.hpp"
15 #include "Xpetra_Access.hpp"
16 #include "Xpetra_Matrix.hpp"
17 
18 namespace MueLu {
19 
24 enum DecisionType : char {
25  UNDECIDED = 0, // no decision has been taken yet, used for initialization
26  KEEP = 1, // keeep the entry
27  DROP = 2, // drop it
28  BOUNDARY = 3 // entry is a boundary
29 };
30 
31 namespace Misc {
32 
33 template <class local_ordinal_type>
34 class NoOpFunctor {
35  public:
37 
38  KOKKOS_FORCEINLINE_FUNCTION
39  void operator()(local_ordinal_type rlid) const {
40  }
41 };
42 
47 template <class local_matrix_type>
49  private:
50  using scalar_type = typename local_matrix_type::value_type;
51  using local_ordinal_type = typename local_matrix_type::ordinal_type;
52  using memory_space = typename local_matrix_type::memory_space;
53  using results_view = Kokkos::View<DecisionType*, memory_space>;
54  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
55 
56  local_matrix_type A;
59 
60  public:
61  PointwiseDropBoundaryFunctor(local_matrix_type& A_, boundary_nodes_view boundaryNodes_, results_view& results_)
62  : A(A_)
63  , boundaryNodes(boundaryNodes_)
64  , results(results_) {}
65 
66  KOKKOS_FORCEINLINE_FUNCTION
67  void operator()(local_ordinal_type rlid) const {
68  auto row = A.rowConst(rlid);
69  const size_t offset = A.graph.row_map(rlid);
70  const bool isBoundaryRow = boundaryNodes(rlid);
71  if (isBoundaryRow) {
72  for (local_ordinal_type k = 0; k < row.length; ++k) {
73  auto clid = row.colidx(k);
74  results(offset + k) = Kokkos::max(rlid == clid ? KEEP : DROP,
75  results(offset + k));
76  }
77  }
78  }
79 };
80 
85 template <class local_matrix_type>
87  private:
88  using scalar_type = typename local_matrix_type::value_type;
89  using local_ordinal_type = typename local_matrix_type::ordinal_type;
90  using memory_space = typename local_matrix_type::memory_space;
91  using results_view = Kokkos::View<DecisionType*, memory_space>;
92  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
93  using block_indices_view_type = Kokkos::View<local_ordinal_type*, memory_space>;
94 
95  local_matrix_type A;
99 
100  public:
101  VectorDropBoundaryFunctor(local_matrix_type& A_, block_indices_view_type point_to_block_, boundary_nodes_view boundaryNodes_, results_view& results_)
102  : A(A_)
103  , point_to_block(point_to_block_)
104  , boundaryNodes(boundaryNodes_)
105  , results(results_) {}
106 
107  KOKKOS_FORCEINLINE_FUNCTION
108  void operator()(local_ordinal_type rlid) const {
109  auto row = A.rowConst(rlid);
110  const size_t offset = A.graph.row_map(rlid);
111  const bool isBoundaryRow = boundaryNodes(point_to_block(rlid));
112  if (isBoundaryRow) {
113  for (local_ordinal_type k = 0; k < row.length; ++k) {
114  auto clid = row.colidx(k);
115  results(offset + k) = Kokkos::max(rlid == clid ? KEEP : DROP,
116  results(offset + k));
117  }
118  }
119  }
120 };
121 
126 template <class local_matrix_type>
128  private:
129  using scalar_type = typename local_matrix_type::value_type;
130  using local_ordinal_type = typename local_matrix_type::ordinal_type;
131  using memory_space = typename local_matrix_type::memory_space;
132  using results_view = Kokkos::View<DecisionType*, memory_space>;
133 
134  local_matrix_type A;
136 
137  public:
138  KeepDiagonalFunctor(local_matrix_type& A_, results_view& results_)
139  : A(A_)
140  , results(results_) {}
141 
142  KOKKOS_FORCEINLINE_FUNCTION
143  void operator()(local_ordinal_type rlid) const {
144  auto row = A.rowConst(rlid);
145  const size_t offset = A.graph.row_map(rlid);
146  for (local_ordinal_type k = 0; k < row.length; ++k) {
147  auto clid = row.colidx(k);
148  if ((rlid == clid) && (results(offset + k) != BOUNDARY)) {
149  results(offset + k) = KEEP;
150  break;
151  }
152  }
153  }
154 };
155 
160 template <class local_matrix_type>
162  private:
163  using scalar_type = typename local_matrix_type::value_type;
164  using local_ordinal_type = typename local_matrix_type::ordinal_type;
165  using memory_space = typename local_matrix_type::memory_space;
166  using results_view = Kokkos::View<DecisionType*, memory_space>;
167 
168  local_matrix_type A;
170 
171  public:
172  DropOffRankFunctor(local_matrix_type& A_, results_view& results_)
173  : A(A_)
174  , results(results_) {}
175 
176  KOKKOS_FORCEINLINE_FUNCTION
177  void operator()(local_ordinal_type rlid) const {
178  auto row = A.rowConst(rlid);
179  const size_t offset = A.graph.row_map(rlid);
180  for (local_ordinal_type k = 0; k < row.length; ++k) {
181  auto clid = row.colidx(k);
182  if (clid >= A.numRows()) {
183  results(offset + k) = Kokkos::max(DROP, results(offset + k));
184  }
185  }
186  }
187 };
188 
193 template <class local_matrix_type>
195  private:
196  using scalar_type = typename local_matrix_type::value_type;
197  using local_ordinal_type = typename local_matrix_type::ordinal_type;
198  using memory_space = typename local_matrix_type::memory_space;
199  using results_view = Kokkos::View<DecisionType*, memory_space>;
200 
201  using boundary_nodes_view = Kokkos::View<bool*, memory_space>;
202 
203  local_matrix_type A;
206 
207  public:
208  MarkSingletonFunctor(local_matrix_type& A_, boundary_nodes_view boundaryNodes_, results_view& results_)
209  : A(A_)
210  , boundaryNodes(boundaryNodes_)
211  , results(results_) {}
212 
213  KOKKOS_FORCEINLINE_FUNCTION
214  void operator()(local_ordinal_type rlid) const {
215  auto row = A.rowConst(rlid);
216  const size_t offset = A.graph.row_map(rlid);
217  for (local_ordinal_type k = 0; k < row.length; ++k) {
218  auto clid = row.colidx(k);
219  if ((results(offset + k) == KEEP) && (rlid != clid))
220  return;
221  }
222  boundaryNodes(rlid) = true;
223  for (local_ordinal_type k = 0; k < row.length; ++k) {
224  auto clid = row.colidx(k);
225  if (rlid == clid)
226  results(offset + k) = KEEP;
227  else
228  results(offset + k) = BOUNDARY;
229  }
230  }
231 };
232 
237 template <class local_matrix_type>
239  private:
240  using scalar_type = typename local_matrix_type::value_type;
241  using local_ordinal_type = typename local_matrix_type::ordinal_type;
242  using memory_space = typename local_matrix_type::memory_space;
243  using results_view = Kokkos::View<DecisionType*, memory_space>;
244  using block_indices_view_type = Kokkos::View<local_ordinal_type*, memory_space>;
245 
246  using boundary_nodes_view = Kokkos::View<bool*, memory_space>;
247 
248  local_matrix_type A;
252 
253  public:
254  MarkSingletonVectorFunctor(local_matrix_type& A_, block_indices_view_type point_to_block_, boundary_nodes_view boundaryNodes_, results_view& results_)
255  : A(A_)
256  , point_to_block(point_to_block_)
257  , boundaryNodes(boundaryNodes_)
258  , results(results_) {}
259 
260  KOKKOS_FORCEINLINE_FUNCTION
261  void operator()(local_ordinal_type rlid) const {
262  auto row = A.rowConst(rlid);
263  const size_t offset = A.graph.row_map(rlid);
264  for (local_ordinal_type k = 0; k < row.length; ++k) {
265  auto clid = row.colidx(k);
266  if ((results(offset + k) == KEEP) && (rlid != clid))
267  return;
268  }
269  auto brlid = point_to_block(rlid);
270  boundaryNodes(brlid) = true;
271  for (local_ordinal_type k = 0; k < row.length; ++k) {
272  auto clid = row.colidx(k);
273  if (rlid == clid)
274  results(offset + k) = KEEP;
275  else
276  results(offset + k) = BOUNDARY;
277  }
278  }
279 };
280 
285 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
287  private:
289  using local_matrix_type = typename matrix_type::local_matrix_type;
290 
291  using scalar_type = typename local_matrix_type::value_type;
292  using local_ordinal_type = typename local_matrix_type::ordinal_type;
293  using memory_space = typename local_matrix_type::memory_space;
294  using results_view = Kokkos::View<DecisionType*, memory_space>;
295 
297  using local_block_indices_view_type = typename block_indices_type::dual_view_type_const::t_dev;
298 
303 
304  public:
306  : A(A_.getLocalMatrixDevice())
307  , point_to_block(point_to_block_.getLocalViewDevice(Xpetra::Access::ReadOnly))
308  , results(results_) {
309  auto importer = A_.getCrsGraph()->getImporter();
310 
311  if (!importer.is_null()) {
312  auto ghosted_point_to_blockMV = Xpetra::VectorFactory<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>::Build(importer->getTargetMap());
313  ghosted_point_to_blockMV->doImport(point_to_block_, *importer, Xpetra::INSERT);
314  ghosted_point_to_block = ghosted_point_to_blockMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
315  } else
317  }
318 
319  KOKKOS_FORCEINLINE_FUNCTION
320  void operator()(local_ordinal_type rlid) const {
321  auto row = A.rowConst(rlid);
322  const size_t offset = A.graph.row_map(rlid);
323  for (local_ordinal_type k = 0; k < row.length; ++k) {
324  auto clid = row.colidx(k);
325  if (point_to_block(rlid, 0) == ghosted_point_to_block(clid, 0)) {
326  results(offset + k) = Kokkos::max(KEEP, results(offset + k));
327  } else {
328  results(offset + k) = Kokkos::max(DROP, results(offset + k));
329  }
330  }
331  }
332 };
333 
338 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
340  private:
342  using local_matrix_type = typename matrix_type::local_matrix_type;
343 
344  using scalar_type = typename local_matrix_type::value_type;
345  using local_ordinal_type = typename local_matrix_type::ordinal_type;
346  using memory_space = typename local_matrix_type::memory_space;
347  using results_view = Kokkos::View<DecisionType*, memory_space>;
348 
350  using local_block_indices_view_type = typename block_indices_type::dual_view_type_const::t_dev;
351  using id_translation_type = Kokkos::View<local_ordinal_type*, memory_space>;
353 
361 
362  public:
363  BlockDiagonalizeVectorFunctor(matrix_type& A_, block_indices_type& point_to_block_, Teuchos::RCP<const map_type> non_unique_map_, results_view& results_, id_translation_type row_translation_, id_translation_type col_translation_)
364  : A(A_.getLocalMatrixDevice())
365  , point_to_block(point_to_block_.getLocalViewDevice(Xpetra::Access::ReadOnly))
366  , results(results_)
367  , row_translation(row_translation_)
368  , col_translation(col_translation_) {
369  auto importer = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(point_to_block_.getMap(), non_unique_map_);
370 
371  if (!importer.is_null()) {
373  ghosted_point_to_blockMV->doImport(point_to_block_, *importer, Xpetra::INSERT);
374  ghosted_point_to_block = ghosted_point_to_blockMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
375  } else
377  }
378 
379  KOKKOS_FORCEINLINE_FUNCTION
380  void operator()(local_ordinal_type rlid) const {
381  auto row = A.rowConst(rlid);
382  const size_t offset = A.graph.row_map(rlid);
383  auto brlid = row_translation(rlid);
384  for (local_ordinal_type k = 0; k < row.length; ++k) {
385  auto clid = row.colidx(k);
386  auto bclid = col_translation(clid);
387  if (point_to_block(brlid, 0) == ghosted_point_to_block(bclid, 0)) {
388  results(offset + k) = Kokkos::max(KEEP, results(offset + k));
389  } else {
390  results(offset + k) = Kokkos::max(DROP, results(offset + k));
391  }
392  }
393  }
394 };
395 
400 template <class local_matrix_type>
402  private:
403  using scalar_type = typename local_matrix_type::value_type;
404  using local_ordinal_type = typename local_matrix_type::ordinal_type;
405  using memory_space = typename local_matrix_type::memory_space;
406  using results_view = Kokkos::View<DecisionType*, memory_space>;
407 
408  using boundary_nodes_view = Kokkos::View<bool*, memory_space>;
409 
410  local_matrix_type A;
412 
413  public:
414  DebugFunctor(local_matrix_type& A_, results_view& results_)
415  : A(A_)
416  , results(results_) {}
417 
418  KOKKOS_FORCEINLINE_FUNCTION
419  void operator()(local_ordinal_type rlid) const {
420  auto row = A.rowConst(rlid);
421  const size_t offset = A.graph.row_map(rlid);
422  for (local_ordinal_type k = 0; k < row.length; ++k) {
423  if (results(offset + k) == UNDECIDED) {
424  Kokkos::printf("No dropping decision was taken for entry (%d, %d)\n", rlid, row.colidx(k));
425  assert(false);
426  }
427  }
428  }
429 };
430 
435 template <class local_matrix_type>
437  private:
438  using scalar_type = typename local_matrix_type::value_type;
439  using local_ordinal_type = typename local_matrix_type::ordinal_type;
440  using memory_space = typename local_matrix_type::memory_space;
441  using results_view = Kokkos::View<DecisionType*, memory_space>;
442 
443  local_matrix_type A;
445 
446  public:
447  SymmetrizeFunctor(local_matrix_type& A_, results_view& results_)
448  : A(A_)
449  , results(results_) {}
450 
451  KOKKOS_FORCEINLINE_FUNCTION
452  void operator()(local_ordinal_type rlid) const {
453  auto row = A.rowConst(rlid);
454  const size_t offset = A.graph.row_map(rlid);
455  for (local_ordinal_type k = 0; k < row.length; ++k) {
456  if (results(offset + k) == KEEP) {
457  auto clid = row.colidx(k);
458  if (clid >= A.numRows())
459  continue;
460  auto row2 = A.rowConst(clid);
461  const size_t offset2 = A.graph.row_map(clid);
462  for (local_ordinal_type k2 = 0; k2 < row2.length; ++k2) {
463  auto clid2 = row2.colidx(k2);
464  if (clid2 == rlid) {
465  if (results(offset2 + k2) == DROP)
466  results(offset2 + k2) = KEEP;
467  break;
468  }
469  }
470  }
471  }
472  }
473 };
474 
475 template <class view_type, class comparator_type>
476 KOKKOS_INLINE_FUNCTION void serialHeapSort(view_type& v, comparator_type comparator) {
477  auto N = v.extent(0);
478  size_t start = N / 2;
479  size_t end = N;
480  while (end > 1) {
481  if (start > 0)
482  start = start - 1;
483  else {
484  end = end - 1;
485  auto temp = v(0);
486  v(0) = v(end);
487  v(end) = temp;
488  }
489  size_t root = start;
490  while (2 * root + 1 < end) {
491  size_t child = 2 * root + 1;
492  if ((child + 1 < end) and (comparator(v(child), v(child + 1))))
493  ++child;
494 
495  if (comparator(v(root), v(child))) {
496  auto temp = v(root);
497  v(root) = v(child);
498  v(child) = temp;
499  root = child;
500  } else
501  break;
502  }
503  }
504 }
505 
508 enum StrengthMeasure : int {
509  /*
510  \f[
511  \frac{|A_{ij}|^2}{|A_{ii}| |A_{jj}|} \le \theta^2
512  \f]
513  */
515  /*
516  \f[
517  \frac{-\operatorname{Re}A_{ij}}{| max_j -A_{ij}|} \le \theta
518  \f]
519  */
521 
522  /*
523  \f[
524  \frac{-\operatorname{sign}(A_{ij}) |A_{ij}|^2}{|A_{ii}| |A_{jj}|} \le \theta^2
525  \f]
526  */
528 };
529 
530 } // namespace Misc
531 
532 } // namespace MueLu
533 
534 #endif
typename matrix_type::local_matrix_type local_matrix_type
Kokkos::View< local_ordinal_type *, memory_space > block_indices_view_type
typename block_indices_type::dual_view_type_const::t_dev local_block_indices_view_type
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
typename local_matrix_type::value_type scalar_type
DebugFunctor(local_matrix_type &A_, results_view &results_)
typename local_matrix_type::ordinal_type local_ordinal_type
typename local_matrix_type::memory_space memory_space
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
typename local_matrix_type::memory_space memory_space
typename local_matrix_type::ordinal_type local_ordinal_type
typename local_matrix_type::value_type scalar_type
typename local_matrix_type::value_type scalar_type
Functor that drops boundary nodes for a blockSize &gt; 1 problem.
Kokkos::View< DecisionType *, memory_space > results_view
Functor that marks singletons (all off-diagonal entries in a row are dropped) as boundary.
typename local_matrix_type::ordinal_type local_ordinal_type
typename local_matrix_type::memory_space memory_space
BlockDiagonalizeFunctor(matrix_type &A_, block_indices_type &point_to_block_, results_view &results_)
typename local_matrix_type::value_type scalar_type
typename local_matrix_type::memory_space memory_space
typename local_matrix_type::ordinal_type local_ordinal_type
local_block_indices_view_type ghosted_point_to_block
typename local_matrix_type::ordinal_type local_ordinal_type
Kokkos::View< DecisionType *, memory_space > results_view
local_block_indices_view_type point_to_block
Functor that checks that all entries have been marked.
typename local_matrix_type::memory_space memory_space
Kokkos::View< bool *, memory_space > boundary_nodes_view
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
Kokkos::View< DecisionType *, memory_space > results_view
typename local_matrix_type::memory_space memory_space
Functor that marks singletons (all off-diagonal entries in a row are dropped) as boundary.
typename local_matrix_type::value_type scalar_type
Kokkos::View< DecisionType *, memory_space > results_view
typename local_matrix_type::ordinal_type local_ordinal_type
Kokkos::View< DecisionType *, memory_space > results_view
typename local_matrix_type::ordinal_type local_ordinal_type
VectorDropBoundaryFunctor(local_matrix_type &A_, block_indices_view_type point_to_block_, boundary_nodes_view boundaryNodes_, results_view &results_)
DropOffRankFunctor(local_matrix_type &A_, results_view &results_)
typename block_indices_type::dual_view_type_const::t_dev local_block_indices_view_type
KOKKOS_INLINE_FUNCTION void serialHeapSort(view_type &v, comparator_type comparator)
typename local_matrix_type::value_type scalar_type
Teuchos::RCP< block_indices_type > ghosted_point_to_blockMV
Functor that symmetrizes the dropping decisions.
Functor that drops boundary nodes for a blockSize == 1 problem.
Functor that drops off-rank entries.
typename local_matrix_type::memory_space memory_space
typename local_matrix_type::memory_space memory_space
Kokkos::View< const bool *, memory_space > boundary_nodes_view
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
typename local_matrix_type::memory_space memory_space
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
typename local_matrix_type::ordinal_type local_ordinal_type
local_block_indices_view_type ghosted_point_to_block
Kokkos::View< bool *, memory_space > boundary_nodes_view
typename matrix_type::local_matrix_type local_matrix_type
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
Functor that drops all entries that are not on the block diagonal.
Functor that marks diagonal as kept, unless the are already marked as boundary.
Kokkos::View< DecisionType *, memory_space > results_view
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
void start()
MarkSingletonFunctor(local_matrix_type &A_, boundary_nodes_view boundaryNodes_, results_view &results_)
typename local_matrix_type::value_type scalar_type
virtual RCP< const CrsGraph > getCrsGraph() const =0
BlockDiagonalizeVectorFunctor(matrix_type &A_, block_indices_type &point_to_block_, Teuchos::RCP< const map_type > non_unique_map_, results_view &results_, id_translation_type row_translation_, id_translation_type col_translation_)
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
typename local_matrix_type::memory_space memory_space
Kokkos::View< local_ordinal_type *, memory_space > id_translation_type
typename local_matrix_type::ordinal_type local_ordinal_type
typename local_matrix_type::value_type scalar_type
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
Kokkos::View< bool *, memory_space > boundary_nodes_view
SymmetrizeFunctor(local_matrix_type &A_, results_view &results_)
Kokkos::View< DecisionType *, memory_space > results_view
MarkSingletonVectorFunctor(local_matrix_type &A_, block_indices_view_type point_to_block_, boundary_nodes_view boundaryNodes_, results_view &results_)
Kokkos::View< const bool *, memory_space > boundary_nodes_view
PointwiseDropBoundaryFunctor(local_matrix_type &A_, boundary_nodes_view boundaryNodes_, results_view &results_)
Kokkos::View< DecisionType *, memory_space > results_view
Functor that drops all entries that are not on the block diagonal.
Kokkos::View< DecisionType *, memory_space > results_view
typename local_matrix_type::ordinal_type local_ordinal_type
Kokkos::View< local_ordinal_type *, memory_space > block_indices_view_type
TransListIter end
typename local_matrix_type::value_type scalar_type
KeepDiagonalFunctor(local_matrix_type &A_, results_view &results_)
Kokkos::View< DecisionType *, memory_space > results_view
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)
typename local_matrix_type::value_type scalar_type