MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_CoalesceDropFactory_kokkos_def.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_COALESCEDROPFACTORY_KOKKOS_DEF_HPP
11 #define MUELU_COALESCEDROPFACTORY_KOKKOS_DEF_HPP
12 
13 #include <Kokkos_Core.hpp>
14 #include <KokkosSparse_CrsMatrix.hpp>
15 #include <sstream>
16 #include <tuple>
17 
18 #include "Xpetra_Matrix.hpp"
19 
20 #include "MueLu_CoalesceDropFactory_kokkos_decl.hpp"
21 
22 #include "MueLu_AmalgamationInfo.hpp"
23 #include "MueLu_Exceptions.hpp"
24 #include "MueLu_Level.hpp"
25 #include "MueLu_LWGraph_kokkos.hpp"
26 #include "MueLu_MasterList.hpp"
27 #include "MueLu_Monitor.hpp"
28 
29 // #define MUELU_COALESCE_DROP_DEBUG 1
30 
31 #include "MueLu_BoundaryDetection.hpp"
32 #include "MueLu_ClassicalDropping.hpp"
33 #include "MueLu_CutDrop.hpp"
34 #include "MueLu_DroppingCommon.hpp"
35 #include "MueLu_DistanceLaplacianDropping.hpp"
36 #include "MueLu_MatrixConstruction.hpp"
37 
38 namespace MueLu {
39 
40 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
41 RCP<const ParameterList> CoalesceDropFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
42  RCP<ParameterList> validParamList = rcp(new ParameterList());
43 
44 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
45  SET_VALID_ENTRY("aggregation: drop tol");
46  SET_VALID_ENTRY("aggregation: use ml scaling of drop tol");
47  SET_VALID_ENTRY("aggregation: Dirichlet threshold");
48  SET_VALID_ENTRY("aggregation: greedy Dirichlet");
49  SET_VALID_ENTRY("aggregation: row sum drop tol");
50  SET_VALID_ENTRY("aggregation: drop scheme");
51  SET_VALID_ENTRY("aggregation: block diagonal: interleaved blocksize");
52  SET_VALID_ENTRY("aggregation: distance laplacian metric");
53  SET_VALID_ENTRY("aggregation: distance laplacian directional weights");
54  SET_VALID_ENTRY("aggregation: dropping may create Dirichlet");
55  SET_VALID_ENTRY("aggregation: distance laplacian algo");
56  SET_VALID_ENTRY("aggregation: classical algo");
57  SET_VALID_ENTRY("aggregation: coloring: localize color graph");
58 
59  SET_VALID_ENTRY("filtered matrix: use lumping");
60  SET_VALID_ENTRY("filtered matrix: reuse graph");
61  SET_VALID_ENTRY("filtered matrix: reuse eigenvalue");
62 
63  SET_VALID_ENTRY("filtered matrix: use root stencil");
64  SET_VALID_ENTRY("filtered matrix: use spread lumping");
65  SET_VALID_ENTRY("filtered matrix: spread lumping diag dom growth factor");
66  SET_VALID_ENTRY("filtered matrix: spread lumping diag dom cap");
67  SET_VALID_ENTRY("filtered matrix: Dirichlet threshold");
68 
69 #undef SET_VALID_ENTRY
70  validParamList->set<bool>("lightweight wrap", true, "Experimental option for lightweight graph access");
71 
72  // "signed classical" is the Ruge-Stuben style (relative to max off-diagonal), "sign classical sa" is the signed version of the sa criterion (relative to the diagonal values)
73  validParamList->getEntry("aggregation: drop scheme").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("signed classical sa", "classical", "distance laplacian", "signed classical", "block diagonal", "block diagonal classical", "block diagonal distance laplacian", "block diagonal signed classical", "block diagonal colored signed classical", "signed classical distance laplacian", "signed classical sa distance laplacian"))));
74  validParamList->getEntry("aggregation: classical algo").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("default", "unscaled cut", "scaled cut", "scaled cut symmetric"))));
75  validParamList->getEntry("aggregation: distance laplacian algo").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("default", "unscaled cut", "scaled cut", "scaled cut symmetric"))));
76  validParamList->getEntry("aggregation: distance laplacian metric").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("unweighted", "material"))));
77 
78  validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory of the matrix A");
79  validParamList->set<RCP<const FactoryBase>>("UnAmalgamationInfo", Teuchos::null, "Generating factory for UnAmalgamationInfo");
80  validParamList->set<RCP<const FactoryBase>>("Coordinates", Teuchos::null, "Generating factory for Coordinates");
81  validParamList->set<RCP<const FactoryBase>>("BlockNumber", Teuchos::null, "Generating factory for BlockNumber");
82  validParamList->set<RCP<const FactoryBase>>("Material", Teuchos::null, "Generating factory for Material");
83 
84  return validParamList;
85 }
86 
87 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
88 void CoalesceDropFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
89  Input(currentLevel, "A");
90  Input(currentLevel, "UnAmalgamationInfo");
91 
92  const ParameterList& pL = GetParameterList();
93  std::string algo = pL.get<std::string>("aggregation: drop scheme");
94  std::string distLaplMetric = pL.get<std::string>("aggregation: distance laplacian metric");
95  if (algo == "distance laplacian" || algo == "block diagonal distance laplacian") {
96  Input(currentLevel, "Coordinates");
97  if (distLaplMetric == "material")
98  Input(currentLevel, "Material");
99  }
100  if (algo == "signed classical sa")
101  ;
102  else if (algo.find("block diagonal") != std::string::npos) {
103  Input(currentLevel, "BlockNumber");
104  }
105 }
106 
107 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
108 void CoalesceDropFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
109  Build(Level& currentLevel) const {
110  auto A = Get<RCP<Matrix>>(currentLevel, "A");
111  TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0, Exceptions::RuntimeError, "A->GetFixedBlockSize() needs to be a multiple of A->GetStorageBlockSize()");
112  LO blkSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
113 
114  std::tuple<GlobalOrdinal, boundary_nodes_type> results;
115  if (blkSize == 1)
116  results = BuildScalar(currentLevel);
117  else
118  results = BuildVector(currentLevel);
119 
120  if (GetVerbLevel() & Statistics1) {
121  GlobalOrdinal numDropped = std::get<0>(results);
122  auto boundaryNodes = std::get<1>(results);
123 
124  GO numLocalBoundaryNodes = 0;
125  GO numGlobalBoundaryNodes = 0;
126 
127  Kokkos::parallel_reduce(
128  "MueLu:CoalesceDropF:Build:bnd", range_type(0, boundaryNodes.extent(0)),
129  KOKKOS_LAMBDA(const LO i, GO& n) {
130  if (boundaryNodes(i))
131  n++;
132  },
133  numLocalBoundaryNodes);
134 
135  auto comm = A->getRowMap()->getComm();
136  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
137 
138  GO numGlobalTotal = A->getGlobalNumEntries();
139  GO numGlobalDropped;
140  MueLu_sumAll(comm, numDropped, numGlobalDropped);
141 
142  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
143  if (numGlobalTotal != 0) {
144  GetOStream(Statistics1) << "Number of dropped entries: "
145  << numGlobalDropped << "/" << numGlobalTotal
146  << " (" << 100 * Teuchos::as<double>(numGlobalDropped) / Teuchos::as<double>(numGlobalTotal) << "%)" << std::endl;
147  }
148  }
149 }
150 
151 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
152 std::tuple<GlobalOrdinal, typename MueLu::LWGraph_kokkos<LocalOrdinal, GlobalOrdinal, Node>::boundary_nodes_type> CoalesceDropFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
153  BuildScalar(Level& currentLevel) const {
154  FactoryMonitor m(*this, "Build", currentLevel);
155 
156  using MatrixType = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
157  using GraphType = Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>;
158  using local_matrix_type = typename MatrixType::local_matrix_type;
159  using local_graph_type = typename GraphType::local_graph_type;
160  using rowptr_type = typename local_graph_type::row_map_type::non_const_type;
161  using entries_type = typename local_graph_type::entries_type::non_const_type;
162  using values_type = typename local_matrix_type::values_type::non_const_type;
163  using device_type = typename Node::device_type;
164  using memory_space = typename device_type::memory_space;
165 
166  typedef Teuchos::ScalarTraits<SC> STS;
167  typedef typename STS::magnitudeType MT;
168  const MT zero = Teuchos::ScalarTraits<MT>::zero();
169 
170  auto A = Get<RCP<Matrix>>(currentLevel, "A");
171 
173  // Process parameterlist
174  const ParameterList& pL = GetParameterList();
175 
176  // Boundary detection
177  const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<double>("aggregation: Dirichlet threshold")));
178  const typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<double>("aggregation: row sum drop tol"));
179  const LocalOrdinal dirichletNonzeroThreshold = 1;
180 
181  // Dropping
182  const std::string algo = pL.get<std::string>("aggregation: drop scheme");
183  std::string classicalAlgoStr = pL.get<std::string>("aggregation: classical algo");
184  std::string distanceLaplacianAlgoStr = pL.get<std::string>("aggregation: distance laplacian algo");
185  std::string distanceLaplacianMetric = pL.get<std::string>("aggregation: distance laplacian metric");
186  MT threshold;
187  // If we're doing the ML-style halving of the drop tol at each level, we do that here.
188  if (pL.get<bool>("aggregation: use ml scaling of drop tol"))
189  threshold = pL.get<double>("aggregation: drop tol") / pow(2.0, currentLevel.GetLevelID());
190  else
191  threshold = as<MT>(pL.get<double>("aggregation: drop tol"));
192  bool aggregationMayCreateDirichlet = pL.get<bool>("aggregation: dropping may create Dirichlet");
193 
194  // Fill
195  const bool lumping = pL.get<bool>("filtered matrix: use lumping");
196  const bool reuseGraph = pL.get<bool>("filtered matrix: reuse graph");
197  const bool reuseEigenvalue = pL.get<bool>("filtered matrix: reuse eigenvalue");
198 
199  const bool useRootStencil = pL.get<bool>("filtered matrix: use root stencil");
200  const bool useSpreadLumping = pL.get<bool>("filtered matrix: use spread lumping");
201 
202  const MT filteringDirichletThreshold = as<MT>(pL.get<double>("filtered matrix: Dirichlet threshold"));
203  TEUCHOS_ASSERT(!useRootStencil);
204  TEUCHOS_ASSERT(!useSpreadLumping);
205 
206  if (algo == "classical")
207  GetOStream(Runtime0) << "algorithm = \"" << algo << "\" classical algorithm = \"" << classicalAlgoStr << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
208  else if (algo == "distance laplacian")
209  GetOStream(Runtime0) << "algorithm = \"" << algo << "\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr << "\" distance laplacian metric = \"" << distanceLaplacianMetric << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
210  else
211  GetOStream(Runtime0) << "algorithm = \"" << algo << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
212 
213  if (((algo == "classical") && (classicalAlgoStr.find("scaled") != std::string::npos)) || ((algo == "distance laplacian") && (distanceLaplacianAlgoStr.find("scaled") != std::string::npos)))
214  TEUCHOS_TEST_FOR_EXCEPTION(threshold > 1.0, Exceptions::RuntimeError, "For cut-drop algorithms, \"aggregation: drop tol\" = " << threshold << ", needs to be <= 1.0");
215 
217  // We perform four sweeps over the rows of A:
218  // Pass 1: detection of boundary nodes
219  // Pass 2: diagonal extraction
220  // Pass 3: drop decision for each entry and construction of the rowptr of the filtered matrix
221  // Pass 4: fill of the filtered matrix
222  //
223  // Pass 1 and 3 apply a sequence of criteria to each row of the matrix.
224 
225  // TODO: We could merge pass 1 and 2.
226 
227  auto crsA = toCrsMatrix(A);
228  auto lclA = crsA->getLocalMatrixDevice();
229  auto range = range_type(0, lclA.numRows());
230 
232  // Pass 1: Detect boundary nodes
233  //
234  // The following criteria are available:
235  // - BoundaryDetection::PointDirichletFunctor
236  // Marks rows as Dirichlet based on value threshold and number of off-diagonal entries
237  // - BoundaryDetection::RowSumFunctor
238  // Marks rows as Dirichlet bases on row-sum criterion
239 
240  // Dirichlet nodes
241  auto boundaryNodes = boundary_nodes_type("boundaryNodes", lclA.numRows()); // initialized to false
242  {
243  SubFactoryMonitor mBoundary(*this, "Boundary detection", currentLevel);
244 
245  // macro that applies boundary detection functors
246 #define MueLu_runBoundaryFunctors(...) \
247  { \
248  auto boundaries = BoundaryDetection::BoundaryFunctor(lclA, __VA_ARGS__); \
249  Kokkos::parallel_for("CoalesceDrop::BoundaryDetection", range, boundaries); \
250  }
251 
252  auto dirichlet_detection = BoundaryDetection::PointDirichletFunctor(lclA, boundaryNodes, dirichletThreshold, dirichletNonzeroThreshold);
253 
254  if (rowSumTol <= 0.) {
255  MueLu_runBoundaryFunctors(dirichlet_detection);
256  } else {
257  auto apply_rowsum = BoundaryDetection::RowSumFunctor(lclA, boundaryNodes, rowSumTol);
258  MueLu_runBoundaryFunctors(dirichlet_detection,
259  apply_rowsum);
260  }
261 #undef MueLu_runBoundaryFunctors
262  }
263  // In what follows, boundaryNodes can still still get modified if aggregationMayCreateDirichlet == true.
264  // Otherwise we're now done with it now.
265 
267  // Pass 2 & 3: Diagonal extraction and determine dropping and construct
268  // rowptr of filtered matrix
269  //
270  // The following criteria are available:
271  // - Misc::PointwiseDropBoundaryFunctor
272  // Drop all rows that have been marked as Dirichlet
273  // - Misc::DropOffRankFunctor
274  // Drop all entries that are off-rank
275  // - ClassicalDropping::DropFunctor
276  // Classical dropping
277  // - DistanceLaplacian::DropFunctor
278  // Distance Laplacian dropping
279  // - Misc::KeepDiagonalFunctor
280  // Mark diagonal as KEEP
281  // - Misc::MarkSingletonFunctor
282  // Mark singletons after dropping as Dirichlet
283  // - Misc::BlockDiagonalizeFunctor
284  // Drop coupling between blocks
285  //
286  // For the block diagonal variants we first block diagonalized and then apply "blocksize = 1" algorithms.
287 
288  // rowptr of filtered A
289  auto filtered_rowptr = rowptr_type("filtered_rowptr", lclA.numRows() + 1);
290  // Number of nonzeros of filtered A
291  LocalOrdinal nnz_filtered = 0;
292  // dropping decisions for each entry
293  auto results = Kokkos::View<DecisionType*, memory_space>("results", lclA.nnz()); // initialized to UNDECIDED
294  {
295  SubFactoryMonitor mDropping(*this, "Dropping decisions", currentLevel);
296 
297  std::string functorLabel = "MueLu::CoalesceDrop::CountEntries";
298 
299  // macro that applied dropping functors
300 #if !defined(HAVE_MUELU_DEBUG)
301 #define MueLu_runDroppingFunctors(...) \
302  { \
303  auto countingFunctor = MatrixConstruction::PointwiseCountingFunctor(lclA, results, filtered_rowptr, __VA_ARGS__); \
304  Kokkos::parallel_scan(functorLabel, range, countingFunctor, nnz_filtered); \
305  }
306 #else
307 #define MueLu_runDroppingFunctors(...) \
308  { \
309  auto debug = Misc::DebugFunctor(lclA, results); \
310  auto countingFunctor = MatrixConstruction::PointwiseCountingFunctor(lclA, results, filtered_rowptr, __VA_ARGS__, debug); \
311  Kokkos::parallel_scan(functorLabel, range, countingFunctor, nnz_filtered); \
312  }
313 #endif
314 
315  auto drop_boundaries = Misc::PointwiseDropBoundaryFunctor(lclA, boundaryNodes, results);
316 
317  if (threshold != zero) {
318  auto preserve_diagonals = Misc::KeepDiagonalFunctor(lclA, results);
319  auto mark_singletons_as_boundary = Misc::MarkSingletonFunctor(lclA, boundaryNodes, results);
320 
321  if (algo == "classical" || algo == "block diagonal classical") {
322  const auto SoC = Misc::SmoothedAggregationMeasure;
323 
324  if (algo == "block diagonal classical") {
325  auto BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
326  auto block_diagonalize = Misc::BlockDiagonalizeFunctor(*A, *BlockNumber, results);
327 
328  if (classicalAlgoStr == "default") {
329  auto classical_dropping = ClassicalDropping::make_drop_functor<SoC>(*A, threshold, results);
330 
331  if (aggregationMayCreateDirichlet) {
332  MueLu_runDroppingFunctors(block_diagonalize,
333  classical_dropping,
334  drop_boundaries,
335  preserve_diagonals,
336  mark_singletons_as_boundary);
337  } else {
338  MueLu_runDroppingFunctors(block_diagonalize,
339  classical_dropping,
340  drop_boundaries,
341  preserve_diagonals);
342  }
343  } else if (classicalAlgoStr == "unscaled cut") {
344  auto comparison = CutDrop::UnscaledComparison(*A, results);
345  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
346 
347  MueLu_runDroppingFunctors(block_diagonalize,
348  drop_boundaries,
349  preserve_diagonals,
350  cut_drop);
351  } else if (classicalAlgoStr == "scaled cut") {
352  auto comparison = CutDrop::make_scaled_comparison_functor<SoC>(*A, results);
353  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
354 
355  MueLu_runDroppingFunctors(block_diagonalize,
356  drop_boundaries,
357  preserve_diagonals,
358  cut_drop);
359  } else if (classicalAlgoStr == "scaled cut symmetric") {
360  auto comparison = CutDrop::make_scaled_comparison_functor<SoC>(*A, results);
361  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
362 
363  MueLu_runDroppingFunctors(block_diagonalize,
364  drop_boundaries,
365  preserve_diagonals,
366  cut_drop);
367 
368  auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
369 
370  MueLu_runDroppingFunctors(symmetrize);
371 
372  } else {
373  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: classical algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << classicalAlgoStr << "\"");
374  }
375  } else {
376  if (classicalAlgoStr == "default") {
377  auto classical_dropping = ClassicalDropping::make_drop_functor<SoC>(*A, threshold, results);
378 
379  if (aggregationMayCreateDirichlet) {
380  MueLu_runDroppingFunctors(classical_dropping,
381  drop_boundaries,
382  preserve_diagonals,
383  mark_singletons_as_boundary);
384  } else {
385  MueLu_runDroppingFunctors(classical_dropping,
386  drop_boundaries,
387  preserve_diagonals);
388  }
389  } else if (classicalAlgoStr == "unscaled cut") {
390  auto comparison = CutDrop::UnscaledComparison(*A, results);
391  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
392 
393  MueLu_runDroppingFunctors(drop_boundaries,
394  preserve_diagonals,
395  cut_drop);
396  } else if (classicalAlgoStr == "scaled cut") {
397  auto comparison = CutDrop::make_scaled_comparison_functor<SoC>(*A, results);
398  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
399 
400  MueLu_runDroppingFunctors(drop_boundaries,
401  preserve_diagonals,
402  cut_drop);
403  } else if (classicalAlgoStr == "scaled cut symmetric") {
404  auto comparison = CutDrop::make_scaled_comparison_functor<SoC>(*A, results);
405  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
406 
407  MueLu_runDroppingFunctors(drop_boundaries,
408  preserve_diagonals,
409  cut_drop);
410 
411  auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
412 
413  MueLu_runDroppingFunctors(symmetrize);
414 
415  } else {
416  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: classical algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << classicalAlgoStr << "\"");
417  }
418  }
419  } else if (algo == "signed classical" || algo == "block diagonal signed classical" || algo == "block diagonal colored signed classical") {
420  const auto SoC = Misc::SignedRugeStuebenMeasure;
421 
422  auto signed_classical_rs_dropping = ClassicalDropping::make_drop_functor<SoC>(*A, threshold, results);
423 
424  if (algo == "block diagonal signed classical" || algo == "block diagonal colored signed classical") {
425  auto BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
426  auto block_diagonalize = Misc::BlockDiagonalizeFunctor(*A, *BlockNumber, results);
427 
428  if (classicalAlgoStr == "default") {
429  if (aggregationMayCreateDirichlet) {
430  MueLu_runDroppingFunctors(block_diagonalize,
431  signed_classical_rs_dropping,
432  drop_boundaries,
433  preserve_diagonals,
434  mark_singletons_as_boundary);
435 
436  } else {
437  MueLu_runDroppingFunctors(block_diagonalize,
438  signed_classical_rs_dropping,
439  drop_boundaries,
440  preserve_diagonals);
441  }
442  } else {
443  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: classical algo\" must be default, not \"" << classicalAlgoStr << "\"");
444  }
445  } else {
446  if (classicalAlgoStr == "default") {
447  if (aggregationMayCreateDirichlet) {
448  MueLu_runDroppingFunctors(signed_classical_rs_dropping,
449  drop_boundaries,
450  preserve_diagonals,
451  mark_singletons_as_boundary);
452 
453  } else {
454  MueLu_runDroppingFunctors(signed_classical_rs_dropping,
455  drop_boundaries,
456  preserve_diagonals);
457  }
458  } else {
459  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: classical algo\" must be default, not \"" << classicalAlgoStr << "\"");
460  }
461  }
462  } else if (algo == "signed classical sa") {
463  const auto SoC = Misc::SignedSmoothedAggregationMeasure;
464  if (classicalAlgoStr == "default") {
465  auto signed_classical_sa_dropping = ClassicalDropping::make_drop_functor<SoC>(*A, threshold, results);
466 
467  if (aggregationMayCreateDirichlet) {
468  MueLu_runDroppingFunctors(signed_classical_sa_dropping,
469  drop_boundaries,
470  preserve_diagonals,
471  mark_singletons_as_boundary);
472 
473  } else {
474  MueLu_runDroppingFunctors(signed_classical_sa_dropping,
475  drop_boundaries,
476  preserve_diagonals);
477  }
478  } else {
479  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: classical algo\" must be default, not \"" << classicalAlgoStr << "\"");
480  }
481  } else if (algo == "distance laplacian" || algo == "block diagonal distance laplacian" || algo == "signed classical distance laplacian" || algo == "signed classical sa distance laplacian") {
482  using doubleMultiVector = Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO>;
483  auto coords = Get<RCP<doubleMultiVector>>(currentLevel, "Coordinates");
484 
485  if (algo == "block diagonal distance laplacian") {
486  const auto SoC = Misc::SmoothedAggregationMeasure;
487 
488  auto BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
489  auto block_diagonalize = Misc::BlockDiagonalizeFunctor(*A, *BlockNumber, results);
490 
491  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
492 
493  if (distanceLaplacianAlgoStr == "default") {
494  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
495 
496  if (aggregationMayCreateDirichlet) {
497  MueLu_runDroppingFunctors(block_diagonalize,
498  dist_laplacian_dropping,
499  drop_boundaries,
500  preserve_diagonals,
501  mark_singletons_as_boundary);
502  } else {
503  MueLu_runDroppingFunctors(block_diagonalize,
504  dist_laplacian_dropping,
505  drop_boundaries,
506  preserve_diagonals);
507  }
508  } else if (distanceLaplacianAlgoStr == "unscaled cut") {
509  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
510  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
511 
512  MueLu_runDroppingFunctors(block_diagonalize,
513  drop_boundaries,
514  preserve_diagonals,
515  cut_drop);
516  } else if (distanceLaplacianAlgoStr == "scaled cut") {
517  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
518  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
519 
520  MueLu_runDroppingFunctors(block_diagonalize,
521  drop_boundaries,
522  preserve_diagonals,
523  cut_drop);
524  } else if (distanceLaplacianAlgoStr == "scaled cut symmetric") {
525  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
526  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
527 
528  MueLu_runDroppingFunctors(block_diagonalize,
529  drop_boundaries,
530  cut_drop,
531  preserve_diagonals);
532 
533  auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
534 
535  MueLu_runDroppingFunctors(symmetrize);
536  } else {
537  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << distanceLaplacianAlgoStr << "\"");
538  }
539  } else if (algo == "distance laplacian") {
540  const auto SoC = Misc::SmoothedAggregationMeasure;
541 
542  if (distanceLaplacianAlgoStr == "default") {
543  if (distanceLaplacianMetric == "unweighted") {
544  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
545  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
546 
547  if (aggregationMayCreateDirichlet) {
548  MueLu_runDroppingFunctors(dist_laplacian_dropping,
549  drop_boundaries,
550  preserve_diagonals,
551  mark_singletons_as_boundary);
552  } else {
553  MueLu_runDroppingFunctors(dist_laplacian_dropping,
554  drop_boundaries,
555  preserve_diagonals);
556  }
557  } else if (distanceLaplacianMetric == "material") {
558  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
559  if (material->getNumVectors() == 1) {
560  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
561 
562  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
563  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
564 
565  if (aggregationMayCreateDirichlet) {
566  MueLu_runDroppingFunctors(dist_laplacian_dropping,
567  drop_boundaries,
568  preserve_diagonals,
569  mark_singletons_as_boundary);
570  } else {
571  MueLu_runDroppingFunctors(dist_laplacian_dropping,
572  drop_boundaries,
573  preserve_diagonals);
574  }
575  } else {
576  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
577 
578  {
579  std::stringstream ss;
580  ss << "material tensor mean =" << std::endl;
581  size_t k = 0;
582  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
583  ss << " ";
584  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
585  ss << material->getVector(k)->meanValue() << " ";
586  ++k;
587  }
588  ss << std::endl;
589  }
590  GetOStream(Runtime0) << ss.str();
591  }
592 
593  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
594  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
595 
596  if (aggregationMayCreateDirichlet) {
597  MueLu_runDroppingFunctors(dist_laplacian_dropping,
598  drop_boundaries,
599  preserve_diagonals,
600  mark_singletons_as_boundary);
601  } else {
602  MueLu_runDroppingFunctors(dist_laplacian_dropping,
603  drop_boundaries,
604  preserve_diagonals);
605  }
606  }
607  }
608 
609  } else if (distanceLaplacianAlgoStr == "unscaled cut") {
610  if (distanceLaplacianMetric == "unweighted") {
611  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
612  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
613  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
614 
615  MueLu_runDroppingFunctors(drop_boundaries,
616  preserve_diagonals,
617  cut_drop);
618  } else if (distanceLaplacianMetric == "material") {
619  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
620  if (material->getNumVectors() == 1) {
621  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
622 
623  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
624  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
625  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
626 
627  MueLu_runDroppingFunctors(drop_boundaries,
628  preserve_diagonals,
629  cut_drop);
630  } else {
631  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
632 
633  {
634  std::stringstream ss;
635  ss << "material tensor mean =" << std::endl;
636  size_t k = 0;
637  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
638  ss << " ";
639  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
640  ss << material->getVector(k)->meanValue() << " ";
641  ++k;
642  }
643  ss << std::endl;
644  }
645  GetOStream(Runtime0) << ss.str();
646  }
647 
648  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
649  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
650  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
651 
652  MueLu_runDroppingFunctors(drop_boundaries,
653  preserve_diagonals,
654  cut_drop);
655  }
656  }
657  } else if (distanceLaplacianAlgoStr == "scaled cut") {
658  if (distanceLaplacianMetric == "unweighted") {
659  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
660  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
661  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
662 
663  MueLu_runDroppingFunctors(drop_boundaries,
664  preserve_diagonals,
665  cut_drop);
666  } else if (distanceLaplacianMetric == "material") {
667  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
668  if (material->getNumVectors() == 1) {
669  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
670 
671  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
672  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
673  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
674 
675  MueLu_runDroppingFunctors(drop_boundaries,
676  preserve_diagonals,
677  cut_drop);
678  } else {
679  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
680 
681  {
682  std::stringstream ss;
683  ss << "material tensor mean =" << std::endl;
684  size_t k = 0;
685  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
686  ss << " ";
687  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
688  ss << material->getVector(k)->meanValue() << " ";
689  ++k;
690  }
691  ss << std::endl;
692  }
693  GetOStream(Runtime0) << ss.str();
694  }
695 
696  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
697  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
698  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
699 
700  MueLu_runDroppingFunctors(drop_boundaries,
701  preserve_diagonals,
702  cut_drop);
703  }
704  }
705  } else if (distanceLaplacianAlgoStr == "scaled cut symmetric") {
706  if (distanceLaplacianMetric == "unweighted") {
707  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
708  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
709  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
710 
711  MueLu_runDroppingFunctors(drop_boundaries,
712  preserve_diagonals,
713  cut_drop);
714  } else if (distanceLaplacianMetric == "material") {
715  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
716  if (material->getNumVectors() == 1) {
717  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
718 
719  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
720  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
721  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
722 
723  MueLu_runDroppingFunctors(drop_boundaries,
724  preserve_diagonals,
725  cut_drop);
726  } else {
727  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
728 
729  {
730  std::stringstream ss;
731  ss << "material tensor mean =" << std::endl;
732  size_t k = 0;
733  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
734  ss << " ";
735  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
736  ss << material->getVector(k)->meanValue() << " ";
737  ++k;
738  }
739  ss << std::endl;
740  }
741  GetOStream(Runtime0) << ss.str();
742  }
743 
744  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
745  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
746  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
747 
748  MueLu_runDroppingFunctors(drop_boundaries,
749  preserve_diagonals,
750  cut_drop);
751  }
752  }
753 
754  auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
755 
756  MueLu_runDroppingFunctors(symmetrize);
757  } else {
758  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << distanceLaplacianAlgoStr << "\"");
759  }
760  } else if (algo == "signed classical distance laplacian") {
761  const auto SoC = Misc::SignedRugeStuebenMeasure;
762  if (distanceLaplacianAlgoStr == "default") {
763  if (distanceLaplacianMetric == "unweighted") {
764  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
765  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
766 
767  if (aggregationMayCreateDirichlet) {
768  MueLu_runDroppingFunctors(dist_laplacian_dropping,
769  drop_boundaries,
770  preserve_diagonals,
771  mark_singletons_as_boundary);
772  } else {
773  MueLu_runDroppingFunctors(dist_laplacian_dropping,
774  drop_boundaries,
775  preserve_diagonals);
776  }
777  } else if (distanceLaplacianMetric == "material") {
778  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
779  if (material->getNumVectors() == 1) {
780  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
781 
782  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
783  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
784 
785  if (aggregationMayCreateDirichlet) {
786  MueLu_runDroppingFunctors(dist_laplacian_dropping,
787  drop_boundaries,
788  preserve_diagonals,
789  mark_singletons_as_boundary);
790  } else {
791  MueLu_runDroppingFunctors(dist_laplacian_dropping,
792  drop_boundaries,
793  preserve_diagonals);
794  }
795  } else {
796  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
797 
798  {
799  std::stringstream ss;
800  ss << "material tensor mean =" << std::endl;
801  size_t k = 0;
802  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
803  ss << " ";
804  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
805  ss << material->getVector(k)->meanValue() << " ";
806  ++k;
807  }
808  ss << std::endl;
809  }
810  GetOStream(Runtime0) << ss.str();
811  }
812 
813  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
814  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
815 
816  if (aggregationMayCreateDirichlet) {
817  MueLu_runDroppingFunctors(dist_laplacian_dropping,
818  drop_boundaries,
819  preserve_diagonals,
820  mark_singletons_as_boundary);
821  } else {
822  MueLu_runDroppingFunctors(dist_laplacian_dropping,
823  drop_boundaries,
824  preserve_diagonals);
825  }
826  }
827  }
828 
829  } else if (distanceLaplacianAlgoStr == "unscaled cut") {
830  if (distanceLaplacianMetric == "unweighted") {
831  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
832  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
833  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
834 
835  MueLu_runDroppingFunctors(drop_boundaries,
836  preserve_diagonals,
837  cut_drop);
838  } else if (distanceLaplacianMetric == "material") {
839  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
840  if (material->getNumVectors() == 1) {
841  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
842 
843  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
844  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
845  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
846 
847  MueLu_runDroppingFunctors(drop_boundaries,
848  preserve_diagonals,
849  cut_drop);
850  } else {
851  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
852 
853  {
854  std::stringstream ss;
855  ss << "material tensor mean =" << std::endl;
856  size_t k = 0;
857  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
858  ss << " ";
859  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
860  ss << material->getVector(k)->meanValue() << " ";
861  ++k;
862  }
863  ss << std::endl;
864  }
865  GetOStream(Runtime0) << ss.str();
866  }
867 
868  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
869  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
870  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
871 
872  MueLu_runDroppingFunctors(drop_boundaries,
873  preserve_diagonals,
874  cut_drop);
875  }
876  }
877  } else if (distanceLaplacianAlgoStr == "scaled cut") {
878  if (distanceLaplacianMetric == "unweighted") {
879  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
880  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
881  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
882 
883  MueLu_runDroppingFunctors(drop_boundaries,
884  preserve_diagonals,
885  cut_drop);
886  } else if (distanceLaplacianMetric == "material") {
887  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
888  if (material->getNumVectors() == 1) {
889  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
890 
891  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
892  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
893  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
894 
895  MueLu_runDroppingFunctors(drop_boundaries,
896  preserve_diagonals,
897  cut_drop);
898  } else {
899  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
900 
901  {
902  std::stringstream ss;
903  ss << "material tensor mean =" << std::endl;
904  size_t k = 0;
905  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
906  ss << " ";
907  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
908  ss << material->getVector(k)->meanValue() << " ";
909  ++k;
910  }
911  ss << std::endl;
912  }
913  GetOStream(Runtime0) << ss.str();
914  }
915 
916  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
917  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
918  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
919 
920  MueLu_runDroppingFunctors(drop_boundaries,
921  preserve_diagonals,
922  cut_drop);
923  }
924  }
925  } else if (distanceLaplacianAlgoStr == "scaled cut symmetric") {
926  if (distanceLaplacianMetric == "unweighted") {
927  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
928  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
929  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
930 
931  MueLu_runDroppingFunctors(drop_boundaries,
932  preserve_diagonals,
933  cut_drop);
934  } else if (distanceLaplacianMetric == "material") {
935  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
936  if (material->getNumVectors() == 1) {
937  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
938 
939  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
940  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
941  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
942 
943  MueLu_runDroppingFunctors(drop_boundaries,
944  preserve_diagonals,
945  cut_drop);
946  } else {
947  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
948 
949  {
950  std::stringstream ss;
951  ss << "material tensor mean =" << std::endl;
952  size_t k = 0;
953  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
954  ss << " ";
955  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
956  ss << material->getVector(k)->meanValue() << " ";
957  ++k;
958  }
959  ss << std::endl;
960  }
961  GetOStream(Runtime0) << ss.str();
962  }
963 
964  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
965  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
966  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
967 
968  MueLu_runDroppingFunctors(drop_boundaries,
969  preserve_diagonals,
970  cut_drop);
971  }
972  }
973 
974  auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
975 
976  MueLu_runDroppingFunctors(symmetrize);
977  } else {
978  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << distanceLaplacianAlgoStr << "\"");
979  }
980  } else if (algo == "signed classical sa distance laplacian") {
981  const auto SoC = Misc::SignedSmoothedAggregationMeasure;
982 
983  if (distanceLaplacianAlgoStr == "default") {
984  if (distanceLaplacianMetric == "unweighted") {
985  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
986  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
987 
988  if (aggregationMayCreateDirichlet) {
989  MueLu_runDroppingFunctors(dist_laplacian_dropping,
990  drop_boundaries,
991  preserve_diagonals,
992  mark_singletons_as_boundary);
993  } else {
994  MueLu_runDroppingFunctors(dist_laplacian_dropping,
995  drop_boundaries,
996  preserve_diagonals);
997  }
998  } else if (distanceLaplacianMetric == "material") {
999  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
1000  if (material->getNumVectors() == 1) {
1001  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
1002 
1003  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
1004  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
1005 
1006  if (aggregationMayCreateDirichlet) {
1007  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1008  drop_boundaries,
1009  preserve_diagonals,
1010  mark_singletons_as_boundary);
1011  } else {
1012  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1013  drop_boundaries,
1014  preserve_diagonals);
1015  }
1016  } else {
1017  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
1018 
1019  {
1020  std::stringstream ss;
1021  ss << "material tensor mean =" << std::endl;
1022  size_t k = 0;
1023  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
1024  ss << " ";
1025  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
1026  ss << material->getVector(k)->meanValue() << " ";
1027  ++k;
1028  }
1029  ss << std::endl;
1030  }
1031  GetOStream(Runtime0) << ss.str();
1032  }
1033 
1034  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
1035  auto dist_laplacian_dropping = DistanceLaplacian::make_drop_functor<SoC>(*A, threshold, dist2, results);
1036 
1037  if (aggregationMayCreateDirichlet) {
1038  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1039  drop_boundaries,
1040  preserve_diagonals,
1041  mark_singletons_as_boundary);
1042  } else {
1043  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1044  drop_boundaries,
1045  preserve_diagonals);
1046  }
1047  }
1048  }
1049 
1050  } else if (distanceLaplacianAlgoStr == "unscaled cut") {
1051  if (distanceLaplacianMetric == "unweighted") {
1052  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
1053  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
1054  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1055 
1056  MueLu_runDroppingFunctors(drop_boundaries,
1057  preserve_diagonals,
1058  cut_drop);
1059  } else if (distanceLaplacianMetric == "material") {
1060  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
1061  if (material->getNumVectors() == 1) {
1062  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
1063 
1064  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
1065  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
1066  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1067 
1068  MueLu_runDroppingFunctors(drop_boundaries,
1069  preserve_diagonals,
1070  cut_drop);
1071  } else {
1072  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
1073 
1074  {
1075  std::stringstream ss;
1076  ss << "material tensor mean =" << std::endl;
1077  size_t k = 0;
1078  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
1079  ss << " ";
1080  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
1081  ss << material->getVector(k)->meanValue() << " ";
1082  ++k;
1083  }
1084  ss << std::endl;
1085  }
1086  GetOStream(Runtime0) << ss.str();
1087  }
1088 
1089  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
1090  auto comparison = CutDrop::UnscaledDistanceLaplacianComparison(*A, dist2, results);
1091  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1092 
1093  MueLu_runDroppingFunctors(drop_boundaries,
1094  preserve_diagonals,
1095  cut_drop);
1096  }
1097  }
1098  } else if (distanceLaplacianAlgoStr == "scaled cut") {
1099  if (distanceLaplacianMetric == "unweighted") {
1100  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
1101  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
1102  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1103 
1104  MueLu_runDroppingFunctors(drop_boundaries,
1105  preserve_diagonals,
1106  cut_drop);
1107  } else if (distanceLaplacianMetric == "material") {
1108  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
1109  if (material->getNumVectors() == 1) {
1110  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
1111 
1112  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
1113  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
1114  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1115 
1116  MueLu_runDroppingFunctors(drop_boundaries,
1117  preserve_diagonals,
1118  cut_drop);
1119  } else {
1120  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
1121 
1122  {
1123  std::stringstream ss;
1124  ss << "material tensor mean =" << std::endl;
1125  size_t k = 0;
1126  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
1127  ss << " ";
1128  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
1129  ss << material->getVector(k)->meanValue() << " ";
1130  ++k;
1131  }
1132  ss << std::endl;
1133  }
1134  GetOStream(Runtime0) << ss.str();
1135  }
1136 
1137  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
1138  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
1139  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1140 
1141  MueLu_runDroppingFunctors(drop_boundaries,
1142  preserve_diagonals,
1143  cut_drop);
1144  }
1145  }
1146  } else if (distanceLaplacianAlgoStr == "scaled cut symmetric") {
1147  if (distanceLaplacianMetric == "unweighted") {
1148  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*A, coords);
1149  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
1150  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1151 
1152  MueLu_runDroppingFunctors(drop_boundaries,
1153  preserve_diagonals,
1154  cut_drop);
1155  } else if (distanceLaplacianMetric == "material") {
1156  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
1157  if (material->getNumVectors() == 1) {
1158  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
1159 
1160  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
1161  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
1162  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1163 
1164  MueLu_runDroppingFunctors(drop_boundaries,
1165  preserve_diagonals,
1166  cut_drop);
1167  } else {
1168  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
1169 
1170  {
1171  std::stringstream ss;
1172  ss << "material tensor mean =" << std::endl;
1173  size_t k = 0;
1174  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
1175  ss << " ";
1176  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
1177  ss << material->getVector(k)->meanValue() << " ";
1178  ++k;
1179  }
1180  ss << std::endl;
1181  }
1182  GetOStream(Runtime0) << ss.str();
1183  }
1184 
1185  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*A, coords, material);
1186  auto comparison = CutDrop::make_scaled_dlap_comparison_functor<SoC>(*A, dist2, results);
1187  auto cut_drop = CutDrop::CutDropFunctor(comparison, threshold);
1188 
1189  MueLu_runDroppingFunctors(drop_boundaries,
1190  preserve_diagonals,
1191  cut_drop);
1192  }
1193  }
1194 
1195  auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
1196 
1197  MueLu_runDroppingFunctors(symmetrize);
1198  } else {
1199  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << distanceLaplacianAlgoStr << "\"");
1200  }
1201  }
1202  } else if (algo == "block diagonal") {
1203  auto BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
1204  auto block_diagonalize = Misc::BlockDiagonalizeFunctor(*A, *BlockNumber, results);
1205 
1206  MueLu_runDroppingFunctors(block_diagonalize);
1207  } else {
1208  TEUCHOS_ASSERT(false);
1209  }
1210  } else {
1211  Kokkos::deep_copy(results, KEEP);
1212  // FIXME: This seems inconsistent
1213  // MueLu_runDroppingFunctors(drop_boundaries);
1214  auto no_op = Misc::NoOpFunctor<LocalOrdinal>();
1215  MueLu_runDroppingFunctors(no_op);
1216  }
1217 #undef MueLu_runDroppingFunctors
1218  }
1219  GO numDropped = lclA.nnz() - nnz_filtered;
1220  // We now know the number of entries of filtered A and have the final rowptr.
1221 
1223  // Pass 4: Create local matrix for filtered A
1224  //
1225  // Dropped entries are optionally lumped to the diagonal.
1226 
1227  RCP<Matrix> filteredA;
1228  RCP<LWGraph_kokkos> graph;
1229  {
1230  SubFactoryMonitor mFill(*this, "Filtered matrix fill", currentLevel);
1231 
1232  local_matrix_type lclFilteredA;
1233  local_graph_type lclGraph;
1234  if (reuseGraph) {
1235  filteredA = MatrixFactory::BuildCopy(A);
1236  lclFilteredA = filteredA->getLocalMatrixDevice();
1237 
1238  auto colidx = entries_type("entries", nnz_filtered);
1239  lclGraph = local_graph_type(colidx, filtered_rowptr);
1240  } else {
1241  auto colidx = entries_type("entries", nnz_filtered);
1242  auto values = values_type("values", nnz_filtered);
1243  lclFilteredA = local_matrix_type("filteredA",
1244  lclA.numRows(), lclA.numCols(),
1245  nnz_filtered,
1246  values, filtered_rowptr, colidx);
1247  }
1248 
1249  if (lumping) {
1250  if (reuseGraph) {
1251  auto fillFunctor = MatrixConstruction::PointwiseFillReuseFunctor<local_matrix_type, local_graph_type, true>(lclA, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1252  Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_reuse", range, fillFunctor);
1253  } else {
1254  auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor<local_matrix_type, true>(lclA, results, lclFilteredA, filteringDirichletThreshold);
1255  Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_noreuse", range, fillFunctor);
1256  }
1257  } else {
1258  if (reuseGraph) {
1259  auto fillFunctor = MatrixConstruction::PointwiseFillReuseFunctor<local_matrix_type, local_graph_type, false>(lclA, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1260  Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_reuse", range, fillFunctor);
1261  } else {
1262  auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor<local_matrix_type, false>(lclA, results, lclFilteredA, filteringDirichletThreshold);
1263  Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_noreuse", range, fillFunctor);
1264  }
1265  }
1266 
1267  if (!reuseGraph)
1268  filteredA = MatrixFactory::Build(lclFilteredA, A->getRowMap(), A->getColMap(), A->getDomainMap(), A->getRangeMap());
1269  filteredA->SetFixedBlockSize(A->GetFixedBlockSize());
1270 
1271  if (reuseEigenvalue) {
1272  // Reuse max eigenvalue from A
1273  // It is unclear what eigenvalue is the best for the smoothing, but we already may have
1274  // the D^{-1}A estimate in A, may as well use it.
1275  // NOTE: ML does that too
1276  filteredA->SetMaxEigenvalueEstimate(A->GetMaxEigenvalueEstimate());
1277  } else {
1278  filteredA->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
1279  }
1280 
1281  if (!reuseGraph) {
1282  // Use graph of filteredA as graph.
1283  lclGraph = filteredA->getCrsGraph()->getLocalGraphDevice();
1284  }
1285  graph = rcp(new LWGraph_kokkos(lclGraph, filteredA->getRowMap(), filteredA->getColMap(), "amalgamated graph of A"));
1286  graph->SetBoundaryNodeMap(boundaryNodes);
1287  }
1288 
1289  LO dofsPerNode = 1;
1290  Set(currentLevel, "DofsPerNode", dofsPerNode);
1291  Set(currentLevel, "Graph", graph);
1292  Set(currentLevel, "A", filteredA);
1293 
1294  return std::make_tuple(numDropped, boundaryNodes);
1295 }
1296 
1297 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1298 std::tuple<GlobalOrdinal, typename MueLu::LWGraph_kokkos<LocalOrdinal, GlobalOrdinal, Node>::boundary_nodes_type> CoalesceDropFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1299  BuildVector(Level& currentLevel) const {
1300  FactoryMonitor m(*this, "Build", currentLevel);
1301 
1302  using MatrixType = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
1303  using GraphType = Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>;
1304  using local_matrix_type = typename MatrixType::local_matrix_type;
1305  using local_graph_type = typename GraphType::local_graph_type;
1306  using rowptr_type = typename local_graph_type::row_map_type::non_const_type;
1307  using entries_type = typename local_graph_type::entries_type::non_const_type;
1308  using values_type = typename local_matrix_type::values_type::non_const_type;
1309  using device_type = typename Node::device_type;
1310  using memory_space = typename device_type::memory_space;
1311 
1312  typedef Teuchos::ScalarTraits<SC> STS;
1313  typedef typename STS::magnitudeType MT;
1314  const MT zero = Teuchos::ScalarTraits<MT>::zero();
1315 
1316  auto A = Get<RCP<Matrix>>(currentLevel, "A");
1317 
1318  /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
1319  blkSize is the number of storage blocks that must kept together during the amalgamation process.
1320 
1321  Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
1322 
1323  numPDEs = blkSize * storageblocksize.
1324 
1325  If numPDEs==1
1326  Matrix is point storage (classical CRS storage). storageblocksize=1 and blkSize=1
1327  No other values makes sense.
1328 
1329  If numPDEs>1
1330  If matrix uses point storage, then storageblocksize=1 and blkSize=numPDEs.
1331  If matrix uses block storage, with block size of n, then storageblocksize=n, and blkSize=numPDEs/n.
1332  Thus far, only storageblocksize=numPDEs and blkSize=1 has been tested.
1333  */
1334 
1335  TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0, Exceptions::RuntimeError, "A->GetFixedBlockSize() needs to be a multiple of A->GetStorageBlockSize()");
1336  LO blkSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
1337 
1338  auto amalInfo = Get<RCP<AmalgamationInfo>>(currentLevel, "UnAmalgamationInfo");
1339 
1340  const RCP<const Map> rowMap = A->getRowMap();
1341  const RCP<const Map> colMap = A->getColMap();
1342 
1343  // build a node row map (uniqueMap = non-overlapping) and a node column map
1344  // (nonUniqueMap = overlapping). The arrays rowTranslation and colTranslation
1345  // stored in the AmalgamationInfo class container contain the local node id
1346  // given a local dof id. The data is calculated in the AmalgamationFactory and
1347  // stored in the variable "UnAmalgamationInfo" (which is of type AmalagamationInfo)
1348  const RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
1349  const RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
1350  Array<LO> rowTranslationArray = *(amalInfo->getRowTranslation()); // TAW should be transform that into a View?
1351  Array<LO> colTranslationArray = *(amalInfo->getColTranslation());
1352 
1353  Kokkos::View<LO*, Kokkos::MemoryUnmanaged>
1354  rowTranslationView(rowTranslationArray.getRawPtr(), rowTranslationArray.size());
1355  Kokkos::View<LO*, Kokkos::MemoryUnmanaged>
1356  colTranslationView(colTranslationArray.getRawPtr(), colTranslationArray.size());
1357 
1358  // get number of local nodes
1359  LO numNodes = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
1360  typedef typename Kokkos::View<LocalOrdinal*, typename Node::device_type> id_translation_type;
1361  id_translation_type rowTranslation("dofId2nodeId", rowTranslationArray.size());
1362  id_translation_type colTranslation("ov_dofId2nodeId", colTranslationArray.size());
1363  Kokkos::deep_copy(rowTranslation, rowTranslationView);
1364  Kokkos::deep_copy(colTranslation, colTranslationView);
1365 
1366  // extract striding information
1367  blkSize = A->GetFixedBlockSize(); //< the full block size (number of dofs per node in strided map)
1368  LocalOrdinal blkId = -1; //< the block id within a strided map or -1 if it is a full block map
1369  LocalOrdinal blkPartSize = A->GetFixedBlockSize(); //< stores block size of part blkId (or the full block size)
1370  if (A->IsView("stridedMaps") == true) {
1371  const RCP<const Map> myMap = A->getRowMap("stridedMaps");
1372  const RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
1373  TEUCHOS_TEST_FOR_EXCEPTION(strMap.is_null() == true, Exceptions::RuntimeError, "Map is not of type stridedMap");
1374  blkSize = Teuchos::as<const LocalOrdinal>(strMap->getFixedBlockSize());
1375  blkId = strMap->getStridedBlockId();
1376  if (blkId > -1)
1377  blkPartSize = Teuchos::as<LocalOrdinal>(strMap->getStridingData()[blkId]);
1378  }
1379 
1380  TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getLocalNumElements() % blkPartSize != 0, MueLu::Exceptions::RuntimeError, "MueLu::CoalesceDropFactory: Number of local elements is " << A->getRowMap()->getLocalNumElements() << " but should be a multiple of " << blkPartSize);
1381 
1383  // Process parameterlist
1384  const ParameterList& pL = GetParameterList();
1385 
1386  // Boundary detection
1387  const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<double>("aggregation: Dirichlet threshold")));
1388  const typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<double>("aggregation: row sum drop tol"));
1389  const LocalOrdinal dirichletNonzeroThreshold = 1;
1390  const bool useGreedyDirichlet = pL.get<bool>("aggregation: greedy Dirichlet");
1391  TEUCHOS_TEST_FOR_EXCEPTION(rowSumTol > zero, MueLu::Exceptions::RuntimeError, "MueLu::CoalesceDropFactory: RowSum is not implemented for vectorial problems.");
1392 
1393  // Dropping
1394  const std::string algo = pL.get<std::string>("aggregation: drop scheme");
1395  std::string classicalAlgoStr = pL.get<std::string>("aggregation: classical algo");
1396  std::string distanceLaplacianAlgoStr = pL.get<std::string>("aggregation: distance laplacian algo");
1397  std::string distanceLaplacianMetric = pL.get<std::string>("aggregation: distance laplacian metric");
1398  MT threshold;
1399  // If we're doing the ML-style halving of the drop tol at each level, we do that here.
1400  if (pL.get<bool>("aggregation: use ml scaling of drop tol"))
1401  threshold = pL.get<double>("aggregation: drop tol") / pow(2.0, currentLevel.GetLevelID());
1402  else
1403  threshold = as<MT>(pL.get<double>("aggregation: drop tol"));
1404  bool aggregationMayCreateDirichlet = pL.get<bool>("aggregation: dropping may create Dirichlet");
1405 
1406  // Fill
1407  const bool lumping = pL.get<bool>("filtered matrix: use lumping");
1408  const bool reuseGraph = pL.get<bool>("filtered matrix: reuse graph");
1409  const bool reuseEigenvalue = pL.get<bool>("filtered matrix: reuse eigenvalue");
1410 
1411  const bool useRootStencil = pL.get<bool>("filtered matrix: use root stencil");
1412  const bool useSpreadLumping = pL.get<bool>("filtered matrix: use spread lumping");
1413 
1414  const MT filteringDirichletThreshold = as<MT>(pL.get<double>("filtered matrix: Dirichlet threshold"));
1415 
1416  TEUCHOS_ASSERT(!useRootStencil);
1417  TEUCHOS_ASSERT(!useSpreadLumping);
1418 
1419  if (algo == "classical") {
1420  GetOStream(Runtime0) << "algorithm = \"" << algo << "\" classical algorithm = \"" << classicalAlgoStr << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
1421  } else if (algo == "distance laplacian") {
1422  GetOStream(Runtime0) << "algorithm = \"" << algo << "\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr << "\" distance laplacian metric = \"" << distanceLaplacianMetric << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
1423  } else
1424  GetOStream(Runtime0) << "algorithm = \"" << algo << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
1425 
1427  // We perform four sweeps over the rows of A:
1428  // Pass 1: detection of boundary nodes
1429  // Pass 2: diagonal extraction
1430  // Pass 3: drop decision for each entry and construction of the rowptr of the filtered matrix
1431  // Pass 4: fill of the filtered matrix
1432  //
1433  // Pass 1 and 3 apply a sequence of criteria to each row of the matrix.
1434 
1435  // TODO: We could merge pass 1 and 2.
1436 
1437  auto crsA = toCrsMatrix(A);
1438  auto lclA = crsA->getLocalMatrixDevice();
1439  auto range = range_type(0, numNodes);
1440 
1442  // Pass 1: Detect boundary nodes
1443  //
1444  // The following criteria are available:
1445  // - BoundaryDetection::VectorDirichletFunctor
1446  // Marks rows as Dirichlet based on value threshold and number of off-diagonal entries
1447 
1448  // Dirichlet nodes
1449  auto boundaryNodes = boundary_nodes_type("boundaryNodes", numNodes); // initialized to false
1450  {
1451  SubFactoryMonitor mBoundary(*this, "Boundary detection", currentLevel);
1452 
1453 #define MueLu_runBoundaryFunctors(...) \
1454  { \
1455  auto boundaries = BoundaryDetection::BoundaryFunctor(lclA, __VA_ARGS__); \
1456  Kokkos::parallel_for("CoalesceDrop::BoundaryDetection", range, boundaries); \
1457  }
1458 
1459  if (useGreedyDirichlet) {
1460  auto dirichlet_detection = BoundaryDetection::VectorDirichletFunctor<local_matrix_type, true>(lclA, blkPartSize, boundaryNodes, dirichletThreshold, dirichletNonzeroThreshold);
1461  MueLu_runBoundaryFunctors(dirichlet_detection);
1462  } else {
1463  auto dirichlet_detection = BoundaryDetection::VectorDirichletFunctor<local_matrix_type, false>(lclA, blkPartSize, boundaryNodes, dirichletThreshold, dirichletNonzeroThreshold);
1464  MueLu_runBoundaryFunctors(dirichlet_detection);
1465  }
1466 #undef MueLu_runBoundaryFunctors
1467  }
1468  // In what follows, boundaryNodes can still still get modified if aggregationMayCreateDirichlet == true.
1469  // Otherwise we're now done with it now.
1470 
1472  // Pass 2 & 3: Diagonal extraction and determine dropping and construct
1473  // rowptr of filtered matrix
1474  //
1475  // The following criteria are available:
1476  // - Misc::VectorDropBoundaryFunctor
1477  // Drop all rows that have been marked as Dirichlet
1478  // - Misc::DropOffRankFunctor
1479  // Drop all entries that are off-rank
1480  // - ClassicalDropping::DropFunctor
1481  // Classical dropping
1482  // - DistanceLaplacian::VectorDropFunctor
1483  // Distance Laplacian dropping
1484  // - Misc::KeepDiagonalFunctor
1485  // Mark diagonal as KEEP
1486  // - Misc::MarkSingletonFunctor
1487  // Mark singletons after dropping as Dirichlet
1488  // - Misc::BlockDiagonalizeVectorFunctor
1489  // Drop coupling between blocks
1490 
1491  // rowptr of filtered A
1492  auto filtered_rowptr = rowptr_type("rowptr", lclA.numRows() + 1);
1493  auto graph_rowptr = rowptr_type("rowptr", numNodes + 1);
1494  // Number of nonzeros of filtered A and graph
1495  Kokkos::pair<LocalOrdinal, LocalOrdinal> nnz = {0, 0};
1496 
1497  // dropping decisions for each entry
1498  auto results = Kokkos::View<DecisionType*, memory_space>("results", lclA.nnz()); // initialized to UNDECIDED
1499  {
1500  SubFactoryMonitor mDropping(*this, "Dropping decisions", currentLevel);
1501 
1502  std::string functorLabel = "MueLu::CoalesceDrop::CountEntries";
1503 
1504 #if !defined(HAVE_MUELU_DEBUG)
1505 #define MueLu_runDroppingFunctors(...) \
1506  { \
1507  auto countingFunctor = MatrixConstruction::VectorCountingFunctor(lclA, blkPartSize, colTranslation, results, filtered_rowptr, graph_rowptr, __VA_ARGS__); \
1508  Kokkos::parallel_scan(functorLabel, range, countingFunctor, nnz); \
1509  }
1510 #else
1511 #define MueLu_runDroppingFunctors(...) \
1512  { \
1513  auto debug = Misc::DebugFunctor(lclA, results); \
1514  auto countingFunctor = MatrixConstruction::VectorCountingFunctor(lclA, blkPartSize, colTranslation, results, filtered_rowptr, graph_rowptr, __VA_ARGS__, debug); \
1515  Kokkos::parallel_scan(functorLabel, range, countingFunctor, nnz); \
1516  }
1517 #endif
1518 
1519  auto drop_boundaries = Misc::VectorDropBoundaryFunctor(lclA, rowTranslation, boundaryNodes, results);
1520 
1521  if (threshold != zero) {
1522  auto preserve_diagonals = Misc::KeepDiagonalFunctor(lclA, results);
1523  auto mark_singletons_as_boundary = Misc::MarkSingletonVectorFunctor(lclA, rowTranslation, boundaryNodes, results);
1524 
1525  if (algo == "classical" || algo == "block diagonal classical") {
1526  const auto SoC = Misc::SmoothedAggregationMeasure;
1527 
1528  if (classicalAlgoStr == "default") {
1529  auto classical_dropping = ClassicalDropping::make_drop_functor<SoC>(*A, threshold, results);
1530 
1531  if (algo == "block diagonal classical") {
1532  auto BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
1533  auto block_diagonalize = Misc::BlockDiagonalizeVectorFunctor(*A, *BlockNumber, nonUniqueMap, results, rowTranslation, colTranslation);
1534 
1535  if (aggregationMayCreateDirichlet) {
1536  MueLu_runDroppingFunctors(block_diagonalize,
1537  classical_dropping,
1538  // drop_boundaries,
1539  preserve_diagonals,
1540  mark_singletons_as_boundary);
1541  } else {
1542  MueLu_runDroppingFunctors(block_diagonalize,
1543  classical_dropping,
1544  // drop_boundaries,
1545  preserve_diagonals);
1546  }
1547  } else {
1548  if (aggregationMayCreateDirichlet) {
1549  MueLu_runDroppingFunctors(classical_dropping,
1550  // drop_boundaries,
1551  preserve_diagonals,
1552  mark_singletons_as_boundary);
1553  } else {
1554  MueLu_runDroppingFunctors(classical_dropping,
1555  // drop_boundaries,
1556  preserve_diagonals);
1557  }
1558  }
1559  } else if (classicalAlgoStr == "unscaled cut") {
1560  TEUCHOS_ASSERT(false);
1561  } else if (classicalAlgoStr == "scaled cut") {
1562  TEUCHOS_ASSERT(false);
1563  } else if (classicalAlgoStr == "scaled cut symmetric") {
1564  TEUCHOS_ASSERT(false);
1565  } else {
1566  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: classical algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << classicalAlgoStr << "\"");
1567  }
1568  } else if (algo == "signed classical" || algo == "block diagonal colored signed classical" || algo == "block diagonal signed classical") {
1569  const auto SoC = Misc::SignedRugeStuebenMeasure;
1570 
1571  auto signed_classical_rs_dropping = ClassicalDropping::make_drop_functor<SoC>(*A, threshold, results);
1572 
1573  if (algo == "block diagonal signed classical" || algo == "block diagonal colored signed classical") {
1574  auto BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
1575  auto block_diagonalize = Misc::BlockDiagonalizeVectorFunctor(*A, *BlockNumber, nonUniqueMap, results, rowTranslation, colTranslation);
1576 
1577  if (aggregationMayCreateDirichlet) {
1578  MueLu_runDroppingFunctors(block_diagonalize,
1579  signed_classical_rs_dropping,
1580  // drop_boundaries,
1581  preserve_diagonals,
1582  mark_singletons_as_boundary);
1583 
1584  } else {
1585  MueLu_runDroppingFunctors(block_diagonalize,
1586  signed_classical_rs_dropping,
1587  // drop_boundaries,
1588  preserve_diagonals);
1589  }
1590  } else {
1591  if (aggregationMayCreateDirichlet) {
1592  MueLu_runDroppingFunctors(signed_classical_rs_dropping,
1593  // drop_boundaries,
1594  preserve_diagonals,
1595  mark_singletons_as_boundary);
1596 
1597  } else {
1598  MueLu_runDroppingFunctors(signed_classical_rs_dropping,
1599  // drop_boundaries,
1600  preserve_diagonals);
1601  }
1602  }
1603  } else if (algo == "signed classical sa") {
1604  const auto SoC = Misc::SignedSmoothedAggregationMeasure;
1605 
1606  auto signed_classical_sa_dropping = ClassicalDropping::make_drop_functor<SoC>(*A, threshold, results);
1607 
1608  if (aggregationMayCreateDirichlet) {
1609  MueLu_runDroppingFunctors(signed_classical_sa_dropping,
1610  // drop_boundaries,
1611  preserve_diagonals,
1612  mark_singletons_as_boundary);
1613 
1614  } else {
1615  MueLu_runDroppingFunctors(signed_classical_sa_dropping,
1616  // drop_boundaries,
1617  preserve_diagonals);
1618  }
1619  } else if (algo == "distance laplacian" || algo == "block diagonal distance laplacian") {
1620  using doubleMultiVector = Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO>;
1621  auto coords = Get<RCP<doubleMultiVector>>(currentLevel, "Coordinates");
1622 
1623  bool use_block_algorithm = (algo == "block diagonal classical" || algo == "block diagonal distance laplacian");
1624  Array<double> dlap_weights = pL.get<Array<double>>("aggregation: distance laplacian directional weights");
1625  enum { NO_WEIGHTS = 0,
1626  SINGLE_WEIGHTS,
1627  BLOCK_WEIGHTS };
1628  int use_dlap_weights = NO_WEIGHTS;
1629  if (algo == "distance laplacian") {
1630  LO dim = (LO)coords->getNumVectors();
1631  // If anything isn't 1.0 we need to turn on the weighting
1632  bool non_unity = false;
1633  for (LO i = 0; !non_unity && i < (LO)dlap_weights.size(); i++) {
1634  if (dlap_weights[i] != 1.0) {
1635  non_unity = true;
1636  }
1637  }
1638  if (non_unity) {
1639  LO blocksize = use_block_algorithm ? as<LO>(pL.get<int>("aggregation: block diagonal: interleaved blocksize")) : 1;
1640  if ((LO)dlap_weights.size() == dim)
1641  use_dlap_weights = SINGLE_WEIGHTS;
1642  else if ((LO)dlap_weights.size() == blocksize * dim)
1643  use_dlap_weights = BLOCK_WEIGHTS;
1644  else {
1645  TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError,
1646  "length of 'aggregation: distance laplacian directional weights' must equal the coordinate dimension OR the coordinate dimension times the blocksize");
1647  }
1648  if (GetVerbLevel() & Statistics1)
1649  GetOStream(Statistics1) << "Using distance laplacian weights: " << dlap_weights << std::endl;
1650  }
1651  }
1652  TEUCHOS_TEST_FOR_EXCEPTION(use_dlap_weights != NO_WEIGHTS, Exceptions::RuntimeError, "Only the NO_WEIGHTS option is implemented for distance laplacian ");
1653 
1654  RCP<Matrix> mergedA;
1655  {
1656  // Construct merged A.
1657 
1658  auto merged_rowptr = rowptr_type("rowptr", numNodes + 1);
1659  LocalOrdinal nnz_merged = 0;
1660 
1661  auto functor = MatrixConstruction::MergeCountFunctor(lclA, blkPartSize, colTranslation, merged_rowptr);
1662  Kokkos::parallel_scan("MergeCount", range, functor, nnz_merged);
1663 
1664  local_graph_type lclMergedGraph;
1665  auto colidx_merged = entries_type("entries", nnz_merged);
1666  auto values_merged = values_type("values", nnz_merged);
1667 
1668  local_matrix_type lclMergedA = local_matrix_type("mergedA",
1669  numNodes, nonUniqueMap->getLocalNumElements(),
1670  nnz_merged,
1671  values_merged, merged_rowptr, colidx_merged);
1672 
1673  auto fillFunctor = MatrixConstruction::MergeFillFunctor<local_matrix_type>(lclA, blkSize, colTranslation, lclMergedA);
1674  Kokkos::parallel_for("MueLu::CoalesceDrop::MergeFill", range, fillFunctor);
1675 
1676  mergedA = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(lclMergedA, uniqueMap, nonUniqueMap, uniqueMap, uniqueMap);
1677  }
1678 
1679  if (algo == "block diagonal distance laplacian") {
1680  auto BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
1681  auto block_diagonalize = Misc::BlockDiagonalizeVectorFunctor(*A, *BlockNumber, nonUniqueMap, results, rowTranslation, colTranslation);
1682 
1683  if (distanceLaplacianAlgoStr == "default") {
1684  const auto SoC = Misc::SmoothedAggregationMeasure;
1685 
1686  if (distanceLaplacianMetric == "unweighted") {
1687  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*mergedA, coords);
1688  auto dist_laplacian_dropping = DistanceLaplacian::make_vector_drop_functor<SoC>(*A, *mergedA, threshold, dist2, results, rowTranslation, colTranslation);
1689 
1690  if (aggregationMayCreateDirichlet) {
1691  MueLu_runDroppingFunctors(block_diagonalize,
1692  dist_laplacian_dropping,
1693  // drop_boundaries,
1694  preserve_diagonals,
1695  mark_singletons_as_boundary);
1696  } else {
1697  MueLu_runDroppingFunctors(block_diagonalize,
1698  dist_laplacian_dropping,
1699  // drop_boundaries,
1700  preserve_diagonals);
1701  }
1702  } else if (distanceLaplacianMetric == "material") {
1703  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
1704  if (material->getNumVectors() == 1) {
1705  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
1706 
1707  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
1708  auto dist_laplacian_dropping = DistanceLaplacian::make_vector_drop_functor<SoC>(*A, *mergedA, threshold, dist2, results, rowTranslation, colTranslation);
1709 
1710  if (aggregationMayCreateDirichlet) {
1711  MueLu_runDroppingFunctors(block_diagonalize,
1712  dist_laplacian_dropping,
1713  drop_boundaries,
1714  preserve_diagonals,
1715  mark_singletons_as_boundary);
1716  } else {
1717  MueLu_runDroppingFunctors(block_diagonalize,
1718  dist_laplacian_dropping,
1719  drop_boundaries,
1720  preserve_diagonals);
1721  }
1722  } else {
1723  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
1724 
1725  {
1726  std::stringstream ss;
1727  ss << "material tensor mean =" << std::endl;
1728  size_t k = 0;
1729  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
1730  ss << " ";
1731  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
1732  ss << material->getVector(k)->meanValue() << " ";
1733  ++k;
1734  }
1735  ss << std::endl;
1736  }
1737  GetOStream(Runtime0) << ss.str();
1738  }
1739 
1740  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*mergedA, coords, material);
1741  auto dist_laplacian_dropping = DistanceLaplacian::make_vector_drop_functor<SoC>(*A, *mergedA, threshold, dist2, results, rowTranslation, colTranslation);
1742 
1743  if (aggregationMayCreateDirichlet) {
1744  MueLu_runDroppingFunctors(block_diagonalize,
1745  dist_laplacian_dropping,
1746  drop_boundaries,
1747  preserve_diagonals,
1748  mark_singletons_as_boundary);
1749  } else {
1750  MueLu_runDroppingFunctors(block_diagonalize,
1751  dist_laplacian_dropping,
1752  drop_boundaries,
1753  preserve_diagonals);
1754  }
1755  }
1756  }
1757  } else if (distanceLaplacianAlgoStr == "unscaled cut") {
1758  TEUCHOS_ASSERT(false);
1759  } else if (distanceLaplacianAlgoStr == "scaled cut") {
1760  TEUCHOS_ASSERT(false);
1761  } else if (distanceLaplacianAlgoStr == "scaled cut symmetric") {
1762  TEUCHOS_ASSERT(false);
1763  } else {
1764  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << distanceLaplacianAlgoStr << "\"");
1765  }
1766  } else {
1767  if (distanceLaplacianAlgoStr == "default") {
1768  const auto SoC = Misc::SmoothedAggregationMeasure;
1769 
1770  if (distanceLaplacianMetric == "unweighted") {
1771  auto dist2 = DistanceLaplacian::UnweightedDistanceFunctor(*mergedA, coords);
1772  auto dist_laplacian_dropping = DistanceLaplacian::make_vector_drop_functor<SoC>(*A, *mergedA, threshold, dist2, results, rowTranslation, colTranslation);
1773 
1774  if (aggregationMayCreateDirichlet) {
1775  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1776  // drop_boundaries,
1777  preserve_diagonals,
1778  mark_singletons_as_boundary);
1779  } else {
1780  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1781  // drop_boundaries,
1782  preserve_diagonals);
1783  }
1784  } else if (distanceLaplacianMetric == "material") {
1785  auto material = Get<RCP<MultiVector>>(currentLevel, "Material");
1786  if (material->getNumVectors() == 1) {
1787  GetOStream(Runtime0) << "material scalar mean = " << material->getVector(0)->meanValue() << std::endl;
1788 
1789  auto dist2 = DistanceLaplacian::ScalarMaterialDistanceFunctor(*A, coords, material);
1790  auto dist_laplacian_dropping = DistanceLaplacian::make_vector_drop_functor<SoC>(*A, *mergedA, threshold, dist2, results, rowTranslation, colTranslation);
1791 
1792  if (aggregationMayCreateDirichlet) {
1793  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1794  drop_boundaries,
1795  preserve_diagonals,
1796  mark_singletons_as_boundary);
1797  } else {
1798  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1799  drop_boundaries,
1800  preserve_diagonals);
1801  }
1802  } else {
1803  TEUCHOS_TEST_FOR_EXCEPTION(coords->getNumVectors() * coords->getNumVectors() != material->getNumVectors(), Exceptions::RuntimeError, "Need \"Material\" to have spatialDim^2 vectors.");
1804 
1805  {
1806  std::stringstream ss;
1807  ss << "material tensor mean =" << std::endl;
1808  size_t k = 0;
1809  for (size_t i = 0; i < coords->getNumVectors(); ++i) {
1810  ss << " ";
1811  for (size_t j = 0; j < coords->getNumVectors(); ++j) {
1812  ss << material->getVector(k)->meanValue() << " ";
1813  ++k;
1814  }
1815  ss << std::endl;
1816  }
1817  GetOStream(Runtime0) << ss.str();
1818  }
1819 
1820  auto dist2 = DistanceLaplacian::TensorMaterialDistanceFunctor(*mergedA, coords, material);
1821  auto dist_laplacian_dropping = DistanceLaplacian::make_vector_drop_functor<SoC>(*A, *mergedA, threshold, dist2, results, rowTranslation, colTranslation);
1822 
1823  if (aggregationMayCreateDirichlet) {
1824  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1825  drop_boundaries,
1826  preserve_diagonals,
1827  mark_singletons_as_boundary);
1828  } else {
1829  MueLu_runDroppingFunctors(dist_laplacian_dropping,
1830  drop_boundaries,
1831  preserve_diagonals);
1832  }
1833  }
1834  }
1835  } else if (distanceLaplacianAlgoStr == "unscaled cut") {
1836  TEUCHOS_ASSERT(false);
1837  } else if (distanceLaplacianAlgoStr == "scaled cut") {
1838  TEUCHOS_ASSERT(false);
1839  } else if (distanceLaplacianAlgoStr == "scaled cut symmetric") {
1840  TEUCHOS_ASSERT(false);
1841  } else {
1842  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut|scaled cut symmetric), not \"" << distanceLaplacianAlgoStr << "\"");
1843  }
1844  }
1845 
1846  } else if (algo == "block diagonal") {
1847  auto BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
1848  auto block_diagonalize = Misc::BlockDiagonalizeVectorFunctor(*A, *BlockNumber, nonUniqueMap, results, rowTranslation, colTranslation);
1849 
1850  MueLu_runDroppingFunctors(block_diagonalize);
1851  } else {
1852  TEUCHOS_ASSERT(false);
1853  }
1854  } else {
1855  Kokkos::deep_copy(results, KEEP);
1856  // MueLu_runDroppingFunctors(drop_boundaries);
1857  auto no_op = Misc::NoOpFunctor<LocalOrdinal>();
1858  MueLu_runDroppingFunctors(no_op);
1859  }
1860 #undef MueLu_runDroppingFunctors
1861  }
1862  LocalOrdinal nnz_filtered = nnz.first;
1863  LocalOrdinal nnz_graph = nnz.second;
1864  GO numTotal = lclA.nnz();
1865  GO numDropped = numTotal - nnz_filtered;
1866  // We now know the number of entries of filtered A and have the final rowptr.
1867 
1869  // Pass 4: Create local matrix for filtered A
1870  //
1871  // Dropped entries are optionally lumped to the diagonal.
1872 
1873  RCP<Matrix> filteredA;
1874  RCP<LWGraph_kokkos> graph;
1875  {
1876  SubFactoryMonitor mFill(*this, "Filtered matrix fill", currentLevel);
1877 
1878  local_matrix_type lclFilteredA;
1879  if (reuseGraph) {
1880  lclFilteredA = local_matrix_type("filteredA", lclA.graph, lclA.numCols());
1881  } else {
1882  auto colidx = entries_type("entries", nnz_filtered);
1883  auto values = values_type("values", nnz_filtered);
1884  lclFilteredA = local_matrix_type("filteredA",
1885  lclA.numRows(), lclA.numCols(),
1886  nnz_filtered,
1887  values, filtered_rowptr, colidx);
1888  }
1889 
1890  local_graph_type lclGraph;
1891  {
1892  auto colidx = entries_type("entries", nnz_graph);
1893  lclGraph = local_graph_type(colidx, graph_rowptr);
1894  }
1895 
1896  if (lumping) {
1897  if (reuseGraph) {
1898  auto fillFunctor = MatrixConstruction::VectorFillFunctor<local_matrix_type, true, true>(lclA, blkPartSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1899  Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_reuse", range, fillFunctor);
1900  } else {
1901  auto fillFunctor = MatrixConstruction::VectorFillFunctor<local_matrix_type, true, false>(lclA, blkPartSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1902  Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_noreuse", range, fillFunctor);
1903  }
1904  } else {
1905  if (reuseGraph) {
1906  auto fillFunctor = MatrixConstruction::VectorFillFunctor<local_matrix_type, false, true>(lclA, blkSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1907  Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_reuse", range, fillFunctor);
1908  } else {
1909  auto fillFunctor = MatrixConstruction::VectorFillFunctor<local_matrix_type, false, false>(lclA, blkSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1910  Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_noreuse", range, fillFunctor);
1911  }
1912  }
1913 
1914  filteredA = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(lclFilteredA, A->getRowMap(), A->getColMap(), A->getDomainMap(), A->getRangeMap());
1915  filteredA->SetFixedBlockSize(blkSize);
1916 
1917  if (reuseEigenvalue) {
1918  // Reuse max eigenvalue from A
1919  // It is unclear what eigenvalue is the best for the smoothing, but we already may have
1920  // the D^{-1}A estimate in A, may as well use it.
1921  // NOTE: ML does that too
1922  filteredA->SetMaxEigenvalueEstimate(A->GetMaxEigenvalueEstimate());
1923  } else {
1924  filteredA->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
1925  }
1926 
1927  graph = rcp(new LWGraph_kokkos(lclGraph, uniqueMap, nonUniqueMap, "amalgamated graph of A"));
1928  graph->SetBoundaryNodeMap(boundaryNodes);
1929  }
1930 
1931  LO dofsPerNode = blkSize;
1932 
1933  Set(currentLevel, "DofsPerNode", dofsPerNode);
1934  Set(currentLevel, "Graph", graph);
1935  Set(currentLevel, "A", filteredA);
1936 
1937  return std::make_tuple(numDropped, boundaryNodes);
1938 }
1939 
1940 } // namespace MueLu
1941 #endif // MUELU_COALESCEDROPFACTORY_KOKKOS_DEF_HPP
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:25
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:13
MueLu::LWGraph_kokkos
Lightweight MueLu representation of a compressed row storage graph.
Definition: MueLu_LWGraph_kokkos_decl.hpp:30
MueLu::LWGraphBase::SetBoundaryNodeMap
KOKKOS_INLINE_FUNCTION void SetBoundaryNodeMap(const boundary_nodes_type bndry)
Set boolean array indicating which rows correspond to Dirichlet boundaries.
Definition: MueLu_LWGraphBase.hpp:234
MueLu_MatrixConstruction.hpp
Teuchos::ParameterEntry::setValidator
void setValidator(RCP< const ParameterEntryValidator > const &validator)
MueLu::Misc::VectorDropBoundaryFunctor
Functor that drops boundary nodes for a blockSize &gt; 1 problem.
Definition: MueLu_DroppingCommon.hpp:86
GO
GlobalOrdinal GO
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
MueLu_Exceptions.hpp
MueLu_DroppingCommon.hpp
MueLu::Misc::MarkSingletonVectorFunctor
Functor that marks singletons (all off-diagonal entries in a row are dropped) as boundary.
Definition: MueLu_DroppingCommon.hpp:238
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:165
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:34
MueLu_runBoundaryFunctors
#define MueLu_runBoundaryFunctors(...)
LO
LocalOrdinal LO
MueLu::Runtime0
One-liner description of what is happening.
Definition: MueLu_VerbosityLevel.hpp:25
MueLu_Level.hpp
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Xpetra_Matrix.hpp
MueLu_MasterList.hpp
MueLu::Misc::MarkSingletonFunctor
Functor that marks singletons (all off-diagonal entries in a row are dropped) as boundary.
Definition: MueLu_DroppingCommon.hpp:194
MueLu_CutDrop.hpp
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::Misc::SymmetrizeFunctor
Functor that symmetrizes the dropping decisions.
Definition: MueLu_DroppingCommon.hpp:436
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:14
MueLu::Misc::PointwiseDropBoundaryFunctor
Functor that drops boundary nodes for a blockSize == 1 problem.
Definition: MueLu_DroppingCommon.hpp:48
MueLu_BoundaryDetection.hpp
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:218
MueLu::CoalesceDropFactory_kokkos::boundary_nodes_type
typename MueLu::LWGraph_kokkos< LocalOrdinal, GlobalOrdinal, Node >::boundary_nodes_type boundary_nodes_type
Definition: MueLu_CoalesceDropFactory_kokkos_decl.hpp:106
MueLu::CoalesceDropFactory_kokkos::DeclareInput
void DeclareInput(Level &currentLevel) const
Input.
Definition: MueLu_CoalesceDropFactory_kokkos_def.hpp:88
MueLu::MatrixConstruction::PointwiseFillReuseFunctor
Functor that fills the filtered matrix while reusing the graph of the matrix before dropping...
Definition: MueLu_MatrixConstruction.hpp:248
MueLu_CoalesceDropFactory_kokkos_decl.hpp
MueLu_ClassicalDropping.hpp
MueLu::CutDrop::CutDropFunctor
Order each row by a criterion, compare the ratio of values and drop all entries once the ratio is bel...
Definition: MueLu_CutDrop.hpp:497
MueLu::BoundaryDetection::PointDirichletFunctor
Functor for marking nodes as Dirichlet.
Definition: MueLu_BoundaryDetection.hpp:34
MueLu::Misc::BlockDiagonalizeFunctor
Functor that drops all entries that are not on the block diagonal.
Definition: MueLu_DroppingCommon.hpp:286
MueLu::Misc::KeepDiagonalFunctor
Functor that marks diagonal as kept, unless the are already marked as boundary.
Definition: MueLu_DroppingCommon.hpp:127
MueLu::CoalesceDropFactory_kokkos::range_type
Kokkos::RangePolicy< local_ordinal_type, execution_space > range_type
Definition: MueLu_CoalesceDropFactory_kokkos_decl.hpp:102
MueLu::CoalesceDropFactory_kokkos::BuildVector
std::tuple< GlobalOrdinal, boundary_nodes_type > BuildVector(Level &currentLevel) const
Definition: MueLu_CoalesceDropFactory_kokkos_def.hpp:1299
MueLu::CoalesceDropFactory_kokkos::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_CoalesceDropFactory_kokkos_def.hpp:41
MueLu::CoalesceDropFactory_kokkos::Build
void Build(Level &currentLevel) const
Build an object with this factory.
Definition: MueLu_CoalesceDropFactory_kokkos_def.hpp:109
MueLu::BoundaryDetection::RowSumFunctor
Functor for marking nodes as Dirichlet based on rowsum.
Definition: MueLu_BoundaryDetection.hpp:149
Teuchos::Array::size
size_type size() const
MueLu_Monitor.hpp
MueLu_DistanceLaplacianDropping.hpp
SET_VALID_ENTRY
#define SET_VALID_ENTRY(name)
MueLu::BoundaryDetection::VectorDirichletFunctor
Functor for marking nodes as Dirichlet in a block operator.
Definition: MueLu_BoundaryDetection.hpp:84
MueLu::CoalesceDropFactory_kokkos::BuildScalar
std::tuple< GlobalOrdinal, boundary_nodes_type > BuildScalar(Level &currentLevel) const
Definition: MueLu_CoalesceDropFactory_kokkos_def.hpp:153
MueLu::MatrixConstruction::PointwiseFillNoReuseFunctor
Functor does not reuse the graph of the matrix for a problem with blockSize == 1. ...
Definition: MueLu_MatrixConstruction.hpp:323
NO
Node NO
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:36
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
MueLu::Misc::BlockDiagonalizeVectorFunctor
Functor that drops all entries that are not on the block diagonal.
Definition: MueLu_DroppingCommon.hpp:339
Xpetra::MultiVector::getNumVectors
virtual size_t getNumVectors() const =0
Teuchos::ParameterList::getEntry
ParameterEntry & getEntry(const std::string &name)
MueLu::CutDrop::UnscaledComparison
Orders entries of row by .
Definition: MueLu_CutDrop.hpp:34
MueLu_runDroppingFunctors
#define MueLu_runDroppingFunctors(...)
Teuchos::RCP::is_null
bool is_null() const