MueLu Version of the Day
Loading...
Searching...
No Matches
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 <string>
17#include <tuple>
18
19#include "Xpetra_Matrix.hpp"
20
22
23#include "MueLu_AmalgamationInfo.hpp"
24#include "MueLu_Exceptions.hpp"
25#include "MueLu_Level.hpp"
26#include "MueLu_LWGraph_kokkos.hpp"
27#include "MueLu_MasterList.hpp"
28#include "MueLu_Monitor.hpp"
29#include "MueLu_Utilities.hpp"
30
34
35// The different dropping algorithms are split up over several translation units. This speeds up compilation and also avoids launch latencies on GPU.
37#include "MueLu_ScalarDroppingClassical.hpp"
38#include "MueLu_ScalarDroppingDistanceLaplacian.hpp"
40#include "MueLu_VectorDroppingClassical.hpp"
41#include "MueLu_VectorDroppingDistanceLaplacian.hpp"
42
43namespace MueLu {
44
45template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
47 RCP<ParameterList> validParamList = rcp(new ParameterList());
48
49#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
50 SET_VALID_ENTRY("aggregation: use blocking");
51 SET_VALID_ENTRY("aggregation: drop tol");
52 SET_VALID_ENTRY("aggregation: use ml scaling of drop tol");
53 SET_VALID_ENTRY("aggregation: Dirichlet threshold");
54 SET_VALID_ENTRY("aggregation: greedy Dirichlet");
55 SET_VALID_ENTRY("aggregation: row sum drop tol");
56 SET_VALID_ENTRY("aggregation: strength-of-connection: matrix");
57 SET_VALID_ENTRY("aggregation: strength-of-connection: measure");
58 SET_VALID_ENTRY("aggregation: drop scheme");
59 SET_VALID_ENTRY("aggregation: block diagonal: interleaved blocksize");
60 SET_VALID_ENTRY("aggregation: distance laplacian metric");
61 SET_VALID_ENTRY("aggregation: distance laplacian directional weights");
62 SET_VALID_ENTRY("aggregation: dropping may create Dirichlet");
63#ifdef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
64 SET_VALID_ENTRY("aggregation: distance laplacian algo");
65 SET_VALID_ENTRY("aggregation: classical algo");
66#endif
67 SET_VALID_ENTRY("aggregation: symmetrize graph after dropping");
68 SET_VALID_ENTRY("aggregation: coloring: use color graph");
69 SET_VALID_ENTRY("aggregation: coloring: localize color graph");
70
71 SET_VALID_ENTRY("filtered matrix: use lumping");
72 SET_VALID_ENTRY("filtered matrix: reuse graph");
73 SET_VALID_ENTRY("filtered matrix: reuse eigenvalue");
74
75 SET_VALID_ENTRY("filtered matrix: use root stencil");
76 SET_VALID_ENTRY("filtered matrix: lumping choice");
77 SET_VALID_ENTRY("filtered matrix: use spread lumping");
78 SET_VALID_ENTRY("filtered matrix: spread lumping diag dom growth factor");
79 SET_VALID_ENTRY("filtered matrix: spread lumping diag dom cap");
80 SET_VALID_ENTRY("filtered matrix: Dirichlet threshold");
81 SET_VALID_ENTRY("filtered matrix: count negative diagonals");
82
83#undef SET_VALID_ENTRY
84 validParamList->set<bool>("lightweight wrap", true, "Experimental option for lightweight graph access");
85#ifndef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
86 validParamList->getEntry("aggregation: drop scheme").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("point-wise", "cut-drop"))));
87#else
88 validParamList->getEntry("aggregation: drop scheme").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("point-wise", "cut-drop", "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"))));
89 validParamList->getEntry("aggregation: classical algo").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("default", "unscaled cut", "scaled cut", "scaled cut symmetric"))));
90 validParamList->getEntry("aggregation: distance laplacian algo").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("default", "unscaled cut", "scaled cut", "scaled cut symmetric"))));
91#endif
92 validParamList->getEntry("aggregation: strength-of-connection: matrix").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("A", "distance laplacian", "MinvA"))));
93 validParamList->getEntry("aggregation: strength-of-connection: measure").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("smoothed aggregation", "signed smoothed aggregation", "signed ruge-stueben", "unscaled"))));
94 validParamList->getEntry("aggregation: distance laplacian metric").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("unweighted", "material"))));
95
96 validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory of the matrix A");
97 validParamList->set<RCP<const FactoryBase>>("UnAmalgamationInfo", Teuchos::null, "Generating factory for UnAmalgamationInfo");
98 validParamList->set<RCP<const FactoryBase>>("Coordinates", Teuchos::null, "Generating factory for Coordinates");
99 validParamList->set<RCP<const FactoryBase>>("BlockNumber", Teuchos::null, "Generating factory for BlockNumber");
100 validParamList->set<RCP<const FactoryBase>>("Material", Teuchos::null, "Generating factory for Material");
101 validParamList->set<RCP<const FactoryBase>>("M", Teuchos::null, "Generating factory for M");
102 validParamList->set<RCP<const FactoryBase>>("Minv", Teuchos::null, "Generating factory for Minv");
103 validParamList->set<RCP<const FactoryBase>>("MinvA", Teuchos::null, "Generating factory for MinvA");
104
105 // Make sure we don't recursively validate options for project auxiliary matrices
106 ParameterList norecurse;
107 norecurse.disableRecursiveValidation();
108 validParamList->set<ParameterList>("project auxiliary matrices", norecurse, "matrices that will be projected on coarse levels");
109
110 return validParamList;
111}
112
113template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
115 Input(currentLevel, "A");
116 Input(currentLevel, "UnAmalgamationInfo");
117
118 const ParameterList& pL = GetParameterList();
119
120 std::string socUsesMatrix = pL.get<std::string>("aggregation: strength-of-connection: matrix");
121 bool needCoords = (socUsesMatrix == "distance laplacian");
122 const bool needM = (socUsesMatrix == "MinvA");
123#ifdef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
124 std::string droppingMethod = pL.get<std::string>("aggregation: drop scheme");
125 needCoords |= (droppingMethod.find("distance laplacian") != std::string::npos);
126#endif
127 if (needCoords) {
128 Input(currentLevel, "Coordinates");
129 std::string distLaplMetric = pL.get<std::string>("aggregation: distance laplacian metric");
130 if (distLaplMetric == "material")
131 Input(currentLevel, "Material");
132 }
133 if (needM && (currentLevel.GetLevelID() != 0)) {
134 if (pL.isSublist("project auxiliary matrices")) {
135 auto projectList = pL.sublist("project auxiliary matrices");
136 if (projectList.isParameter("M")) Input(currentLevel, "M");
137 if (projectList.isParameter("Minv")) Input(currentLevel, "Minv");
138 if (projectList.isParameter("MinvA")) Input(currentLevel, "MinvA");
139 }
140 }
141
142 bool useBlocking = pL.get<bool>("aggregation: use blocking");
143#ifdef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
144 useBlocking |= (droppingMethod.find("block diagonal") != std::string::npos);
145#endif
146 if (useBlocking) {
147 Input(currentLevel, "BlockNumber");
148 }
149}
150
151template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
153 Build(Level& currentLevel) const {
154 auto A = Get<RCP<Matrix>>(currentLevel, "A");
155 TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0, Exceptions::RuntimeError, "A->GetFixedBlockSize() needs to be a multiple of A->GetStorageBlockSize()");
156 LO blkSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
157
158 std::tuple<GlobalOrdinal, boundary_nodes_type> results;
159 if (blkSize == 1)
160 results = BuildScalar(currentLevel);
161 else
162 results = BuildVector(currentLevel);
163
164 if (GetVerbLevel() & Statistics1) {
165 GlobalOrdinal numDropped = std::get<0>(results);
166 auto boundaryNodes = std::get<1>(results);
167
168 GO numLocalBoundaryNodes = 0;
169
170 Kokkos::parallel_reduce(
171 "MueLu:CoalesceDropF:Build:bnd", range_type(0, boundaryNodes.extent(0)),
172 KOKKOS_LAMBDA(const LO i, GO& n) {
173 if (boundaryNodes(i))
174 n++;
175 },
176 numLocalBoundaryNodes);
177
178 if (IsPrint(Statistics1)) {
179 auto comm = A->getRowMap()->getComm();
180
181 std::vector<GlobalOrdinal> localStats = {numLocalBoundaryNodes, numDropped};
182 std::vector<GlobalOrdinal> globalStats(2);
183 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, 2, localStats.data(), globalStats.data());
184
185 GO numGlobalTotal = A->getGlobalNumEntries();
186 GO numGlobalBoundaryNodes = globalStats[0];
187 GO numGlobalDropped = globalStats[1];
188
189 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
190 if (numGlobalTotal != 0) {
191 GetOStream(Statistics1) << "Number of dropped entries: "
192 << numGlobalDropped << "/" << numGlobalTotal
193 << " (" << 100 * Teuchos::as<double>(numGlobalDropped) / Teuchos::as<double>(numGlobalTotal) << "%)" << std::endl;
194 }
195 }
196 }
197}
198
199template <class local_matrix_type, class boundary_nodes_view, class... Functors>
200void runBoundaryFunctors(local_matrix_type& lclA, boundary_nodes_view& boundaryNodes, Functors&... functors) {
201 using local_ordinal_type = typename local_matrix_type::ordinal_type;
202 using execution_space = typename local_matrix_type::execution_space;
203 using range_type = Kokkos::RangePolicy<local_ordinal_type, execution_space>;
204 auto range = range_type(0, boundaryNodes.extent(0));
205 auto boundaries = BoundaryDetection::BoundaryFunctor(lclA, functors...);
206 Kokkos::parallel_for("CoalesceDrop::BoundaryDetection", range, boundaries);
207}
208
209template <class magnitudeType>
210void translateOldAlgoParam(const Teuchos::ParameterList& pL, std::string& droppingMethod, bool& useBlocking, std::string& socUsesMatrix, std::string& socUsesMeasure, bool& symmetrizeDroppedGraph, bool& generateColoringGraph, magnitudeType& threshold, MueLu::MatrixConstruction::lumpingType& lumpingChoice) {
211 std::set<std::string> validDroppingMethods = {"piece-wise", "cut-drop"};
212
213 if (!pL.get<bool>("filtered matrix: use lumping")) lumpingChoice = MueLu::MatrixConstruction::no_lumping;
214
215 if (validDroppingMethods.find(droppingMethod) == validDroppingMethods.end()) {
216 std::string algo = droppingMethod;
217 std::string classicalAlgoStr = pL.get<std::string>("aggregation: classical algo");
218 std::string distanceLaplacianAlgoStr = pL.get<std::string>("aggregation: distance laplacian algo");
219
220 // Remove prefix "block diagonal" from algo
221 if (algo.find("block diagonal") == 0) {
222 useBlocking = true;
223 algo = algo.substr(14);
224 if (algo != "") {
225 algo = algo.substr(1);
226 }
227 }
228
229 if ((algo == "classical") || (algo == "signed classical sa") || (algo == "signed classical") || (algo == "colored signed classical")) {
230 socUsesMatrix = "A";
231
232 if (algo == "classical") {
233 socUsesMeasure = "smoothed aggregation";
234 } else if (algo == "signed classical sa") {
235 socUsesMeasure = "signed smoothed aggregation";
236 } else if (algo == "signed classical") {
237 socUsesMeasure = "signed ruge-stueben";
238 } else if (algo == "colored signed classical") {
239 socUsesMeasure = "signed ruge-stueben";
240 generateColoringGraph = true;
241 }
242
243 if (classicalAlgoStr == "default")
244 droppingMethod = "point-wise";
245 else if (classicalAlgoStr == "unscaled cut") {
246 socUsesMeasure = "unscaled";
247 droppingMethod = "cut-drop";
248 } else if (classicalAlgoStr == "scaled cut") {
249 droppingMethod = "cut-drop";
250 } else if (classicalAlgoStr == "scaled cut symmetric") {
251 droppingMethod = "cut-drop";
252 symmetrizeDroppedGraph = true;
253 }
254 } else if ((algo == "distance laplacian") || (algo == "signed classical sa distance laplacian") || (algo == "signed classical distance laplacian")) {
255 socUsesMatrix = "distance laplacian";
256
257 if (algo == "distance laplacian") {
258 socUsesMeasure = "smoothed aggregation";
259 } else if (algo == "signed classical sa distance laplacian") {
260 socUsesMeasure = "signed smoothed aggregation";
261 } else if (algo == "signed classical distance laplacian") {
262 socUsesMeasure = "signed ruge-stueben";
263 }
264
265 if (distanceLaplacianAlgoStr == "default")
266 droppingMethod = "point-wise";
267 else if (distanceLaplacianAlgoStr == "unscaled cut") {
268 socUsesMeasure = "unscaled";
269 droppingMethod = "cut-drop";
270 } else if (distanceLaplacianAlgoStr == "scaled cut") {
271 droppingMethod = "cut-drop";
272 } else if (distanceLaplacianAlgoStr == "scaled cut symmetric") {
273 droppingMethod = "cut-drop";
274 symmetrizeDroppedGraph = true;
275 }
276 } else if (algo == "") {
277 // algo was "block diagonal", but we process and remove the "block diagonal" part
278 socUsesMatrix = "A";
279 threshold = Teuchos::ScalarTraits<magnitudeType>::zero();
280 }
281 }
282}
283
284template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
285std::tuple<GlobalOrdinal, typename MueLu::LWGraph_kokkos<LocalOrdinal, GlobalOrdinal, Node>::boundary_nodes_type> CoalesceDropFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
286 BuildScalar(Level& currentLevel) const {
287 FactoryMonitor m(*this, "BuildScalar", currentLevel);
288
289 using MatrixType = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
290 using GraphType = Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>;
291 using local_matrix_type = typename MatrixType::local_matrix_device_type;
292 using local_graph_type = typename GraphType::local_graph_device_type;
293 using rowptr_type = typename local_graph_type::row_map_type::non_const_type;
294 using entries_type = typename local_graph_type::entries_type::non_const_type;
295 using values_type = typename local_matrix_type::values_type::non_const_type;
296 using device_type = typename Node::device_type;
297 using memory_space = typename device_type::memory_space;
298 using results_view_type = Kokkos::View<DecisionType*, memory_space>;
299 using magnitudeType = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
300 using doubleMultiVector = Xpetra::MultiVector<magnitudeType, LO, GO, NO>;
301
302 typedef Teuchos::ScalarTraits<Scalar> STS;
303 const magnitudeType zero = Teuchos::ScalarTraits<magnitudeType>::zero();
304
305 auto A = Get<RCP<Matrix>>(currentLevel, "A");
306
308 // Process parameterlist
309 const ParameterList& pL = GetParameterList();
310
311 // Boundary detection
312 const magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<double>("aggregation: Dirichlet threshold")));
313 const magnitudeType rowSumTol = as<magnitudeType>(pL.get<double>("aggregation: row sum drop tol"));
314 const LocalOrdinal dirichletNonzeroThreshold = 1;
315
316 // Dropping
317 bool useBlocking = pL.get<bool>("aggregation: use blocking");
318 std::string droppingMethod = pL.get<std::string>("aggregation: drop scheme");
319 std::string socUsesMatrix = pL.get<std::string>("aggregation: strength-of-connection: matrix");
320 std::string socUsesMeasure = pL.get<std::string>("aggregation: strength-of-connection: measure");
321 std::string distanceLaplacianMetric = pL.get<std::string>("aggregation: distance laplacian metric");
322 bool symmetrizeDroppedGraph = pL.get<bool>("aggregation: symmetrize graph after dropping");
323 magnitudeType threshold;
324 // If we're doing the ML-style halving of the drop tol at each level, we do that here.
325 if (pL.get<bool>("aggregation: use ml scaling of drop tol"))
326 threshold = pL.get<double>("aggregation: drop tol") / pow(2.0, currentLevel.GetLevelID());
327 else
328 threshold = as<magnitudeType>(pL.get<double>("aggregation: drop tol"));
329 bool aggregationMayCreateDirichlet = pL.get<bool>("aggregation: dropping may create Dirichlet");
330
331 // Fill
332 const bool reuseGraph = pL.get<bool>("filtered matrix: reuse graph");
333 const bool reuseEigenvalue = pL.get<bool>("filtered matrix: reuse eigenvalue");
334
335 const bool useRootStencil = pL.get<bool>("filtered matrix: use root stencil");
336 const bool useSpreadLumping = pL.get<bool>("filtered matrix: use spread lumping");
337 const std::string lumpingChoiceString = pL.get<std::string>("filtered matrix: lumping choice");
339 if (lumpingChoiceString == "diag lumping")
341 else if (lumpingChoiceString == "distributed lumping")
343
344 const magnitudeType filteringDirichletThreshold = as<magnitudeType>(pL.get<double>("filtered matrix: Dirichlet threshold"));
345
346 // coloring graph
347 bool generateColoringGraph = pL.get<bool>("aggregation: coloring: use color graph");
348 const bool localizeColoringGraph = pL.get<bool>("aggregation: coloring: localize color graph");
349 const bool symmetrizeColoringGraph = true;
350
351#ifdef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
352 translateOldAlgoParam(pL, droppingMethod, useBlocking, socUsesMatrix, socUsesMeasure, symmetrizeDroppedGraph, generateColoringGraph, threshold, lumpingChoice);
353#endif
354
355 {
356 std::stringstream ss;
357 ss << "dropping scheme = \"" << droppingMethod << "\", strength-of-connection measure = \"" << socUsesMeasure << "\", strength-of-connection matrix = \"" << socUsesMatrix << "\", ";
358 if (socUsesMatrix == "distance laplacian")
359 ss << "distance laplacian metric = \"" << distanceLaplacianMetric << "\", ";
360 ss << "threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << ", useBlocking = " << useBlocking;
361 ss << ", symmetrizeDroppedGraph = " << symmetrizeDroppedGraph << std::endl;
362
363 GetOStream(Runtime0) << ss.str();
364 }
365
366 TEUCHOS_ASSERT(!useRootStencil);
367 TEUCHOS_ASSERT(!useSpreadLumping);
368 TEUCHOS_ASSERT((lumpingChoice != MueLu::MatrixConstruction::distributed_lumping) || !reuseGraph);
369 if (droppingMethod == "cut-drop")
370 TEUCHOS_TEST_FOR_EXCEPTION(threshold > 1.0, Exceptions::RuntimeError, "For cut-drop algorithms, \"aggregation: drop tol\" = " << threshold << ", needs to be <= 1.0");
371
373 // We perform four sweeps over the rows of A:
374 // Pass 1: detection of boundary nodes
375 // Pass 2: diagonal extraction
376 // Pass 3: drop decision for each entry and construction of the rowptr of the filtered matrix
377 // Pass 4: fill of the filtered matrix
378 //
379 // Pass 1 and 3 apply a sequence of criteria to each row of the matrix.
380
381 // TODO: We could merge pass 1 and 2.
382
383 // Compute matrix that is used for dropping
384 RCP<Matrix> A_drop;
385 if (threshold != zero) {
386 if ((socUsesMatrix == "A") || (socUsesMatrix == "MinvA")) {
387 // Get matrix used for dropping
388 if (socUsesMatrix == "A")
389 A_drop = A;
390 else if (socUsesMatrix == "MinvA") {
391 bool storeMinvOnLevel = false, storeMinvAOnLevel = false;
392 if (pL.isSublist("project auxiliary matrices")) {
393 auto projectList = pL.sublist("project auxiliary matrices");
394 if (projectList.isParameter("Minv")) storeMinvOnLevel = true;
395 if (projectList.isParameter("MinvA")) storeMinvAOnLevel = true;
396 // project list appears to be empty, so default behavior is to store MinvA
397 if (!projectList.isParameter("M") && !storeMinvOnLevel && !storeMinvAOnLevel) storeMinvAOnLevel = true;
398 } else
399 storeMinvAOnLevel = true; // default behavior is to project MinvA
400
401 if (IsAvailable(currentLevel, "MinvA")) {
402 A_drop = Get<RCP<Matrix>>(currentLevel, "MinvA");
403 } else {
404 RCP<Matrix> Minv;
405 if (IsAvailable(currentLevel, "Minv")) {
406 Minv = Get<RCP<Matrix>>(currentLevel, "Minv");
407 } else { // get M and create Minv
408 if (currentLevel.GetLevelID() == 0) {
409 auto M = currentLevel.Get<RCP<Matrix>>("M", NoFactory::get());
410 // Create Minv via sparse approximate inverse
411 Minv = Utilities::SPAI(M);
412 } else {
413 auto M = Get<RCP<Matrix>>(currentLevel, "M");
414 // Create Minv via sparse approximate inverse
415 Minv = Utilities::SPAI(M);
416 }
417 if (storeMinvOnLevel) currentLevel.Set("Minv", Minv);
418 } // finished if/else (currentLevel.IsAvailable("Minv", *mtf)
419
420 // build MinvA matrix with same sparsity pattern as A.
421 A_drop = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(A);
422 auto params = Teuchos::rcp(new Teuchos::ParameterList());
423 params->set("MM Throw For Non-Existent Entries", false);
424 A_drop = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*Minv, false, *A, false, A_drop, GetOStream(Statistics2), true, true, std::string("MinvA"), params);
425 if (storeMinvAOnLevel) currentLevel.Set("MinvA", A_drop);
426 } // finished if/else (currentLevel.IsAvailable("MinvA", NoFactory::get()))
427 } // else if (socUsesMatrix == "MinvA") {
428 }
429 }
430
431 auto crsA = toCrsMatrix(A);
432 auto lclA = crsA->getLocalMatrixDevice();
433 auto range = range_type(0, lclA.numRows());
434
436 // Pass 1: Detect boundary nodes
437 //
438 // The following criteria are available:
439 // - BoundaryDetection::PointDirichletFunctor
440 // Marks rows as Dirichlet based on value threshold and number of off-diagonal entries
441 // - BoundaryDetection::RowSumFunctor
442 // Marks rows as Dirichlet bases on row-sum criterion
443
444 // Dirichlet nodes
445 auto boundaryNodes = boundary_nodes_type("boundaryNodes", lclA.numRows()); // initialized to false
446 {
447 SubFactoryMonitor mBoundary(*this, "Boundary detection", currentLevel);
448
449 // macro that applies boundary detection functors
450 auto dirichlet_detection = BoundaryDetection::PointDirichletFunctor(lclA, boundaryNodes, dirichletThreshold, dirichletNonzeroThreshold);
451
452 if (rowSumTol <= 0.) {
453 runBoundaryFunctors(lclA, boundaryNodes, dirichlet_detection);
454 } else {
455 auto apply_rowsum = BoundaryDetection::RowSumFunctor(lclA, boundaryNodes, rowSumTol);
456 runBoundaryFunctors(lclA, boundaryNodes, dirichlet_detection, apply_rowsum);
457 }
458 }
459 // In what follows, boundaryNodes can still still get modified if aggregationMayCreateDirichlet == true.
460 // Otherwise we're now done with it now.
461
463 // Pass 2 & 3: Diagonal extraction and determine dropping and construct
464 // rowptr of filtered matrix
465 //
466 // The following criteria are available:
467 // - Misc::PointwiseDropBoundaryFunctor
468 // Drop all rows that have been marked as Dirichlet
469 // - Misc::DropOffRankFunctor
470 // Drop all entries that are off-rank
471 // - ClassicalDropping::DropFunctor
472 // Classical dropping
473 // - DistanceLaplacian::DropFunctor
474 // Distance Laplacian dropping
475 // - Misc::KeepDiagonalFunctor
476 // Mark diagonal as KEEP
477 // - Misc::MarkSingletonFunctor
478 // Mark singletons after dropping as Dirichlet
479 // - Misc::BlockDiagonalizeFunctor
480 // Drop coupling between blocks
481 //
482 // For the block diagonal variants we first block diagonalized and then apply "blocksize = 1" algorithms.
483
484 // rowptr of filtered A
485 auto filtered_rowptr = rowptr_type("filtered_rowptr", lclA.numRows() + 1);
486 // Number of nonzeros of filtered A
487 LocalOrdinal nnz_filtered = 0;
488 // dropping decisions for each entry
489 auto results = results_view_type("results", lclA.nnz()); // initialized to UNDECIDED
490 {
491 SubFactoryMonitor mDropping(*this, "Dropping decisions", currentLevel);
492
493 if (threshold != zero) {
494 if ((socUsesMatrix == "A") || (socUsesMatrix == "MinvA")) {
495 if (socUsesMeasure == "unscaled") {
496 ScalarDroppingClassical<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::UnscaledMeasure>::runDroppingFunctors_on_A(*A_drop, results, filtered_rowptr, nnz_filtered, boundaryNodes, droppingMethod, threshold,
497 aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, currentLevel, *this);
498 } else if (socUsesMeasure == "smoothed aggregation") {
499 ScalarDroppingClassical<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SmoothedAggregationMeasure>::runDroppingFunctors_on_A(*A_drop, results, filtered_rowptr, nnz_filtered, boundaryNodes, droppingMethod, threshold,
500 aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, currentLevel, *this);
501 } else if (socUsesMeasure == "signed ruge-stueben") {
502 ScalarDroppingClassical<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SignedRugeStuebenMeasure>::runDroppingFunctors_on_A(*A_drop, results, filtered_rowptr, nnz_filtered, boundaryNodes, droppingMethod, threshold,
503 aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, currentLevel, *this);
504 } else if (socUsesMeasure == "signed smoothed aggregation") {
505 ScalarDroppingClassical<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SignedSmoothedAggregationMeasure>::runDroppingFunctors_on_A(*A_drop, results, filtered_rowptr, nnz_filtered, boundaryNodes, droppingMethod, threshold,
506 aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, currentLevel, *this);
507 }
508 } else if (socUsesMatrix == "distance laplacian") {
509 auto coords = Get<RCP<doubleMultiVector>>(currentLevel, "Coordinates");
510 if (socUsesMeasure == "unscaled") {
511 ScalarDroppingDistanceLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::UnscaledMeasure>::runDroppingFunctors_on_dlap(*A, results, filtered_rowptr, nnz_filtered, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, distanceLaplacianMetric, currentLevel, *this);
512 } else if (socUsesMeasure == "smoothed aggregation") {
513 ScalarDroppingDistanceLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SmoothedAggregationMeasure>::runDroppingFunctors_on_dlap(*A, results, filtered_rowptr, nnz_filtered, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, distanceLaplacianMetric, currentLevel, *this);
514 } else if (socUsesMeasure == "signed ruge-stueben") {
515 ScalarDroppingDistanceLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SignedRugeStuebenMeasure>::runDroppingFunctors_on_dlap(*A, results, filtered_rowptr, nnz_filtered, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, distanceLaplacianMetric, currentLevel, *this);
516 } else if (socUsesMeasure == "signed smoothed aggregation") {
517 ScalarDroppingDistanceLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SignedSmoothedAggregationMeasure>::runDroppingFunctors_on_dlap(*A, results, filtered_rowptr, nnz_filtered, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, distanceLaplacianMetric, currentLevel, *this);
518 }
519 }
520 } else {
521 Kokkos::deep_copy(results, KEEP);
522
523 if (symmetrizeDroppedGraph) {
524 auto drop_boundaries = Misc::PointwiseSymmetricDropBoundaryFunctor(*A, boundaryNodes, results);
525 ScalarDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, results, filtered_rowptr, nnz_filtered, useBlocking, currentLevel, *this, drop_boundaries);
526 } else {
527 auto no_op = Misc::NoOpFunctor<LocalOrdinal>();
528 ScalarDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, results, filtered_rowptr, nnz_filtered, useBlocking, currentLevel, *this, no_op);
529 }
530 }
531
532 if (symmetrizeDroppedGraph) {
533 auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
534 ScalarDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, results, filtered_rowptr, nnz_filtered, useBlocking, currentLevel, *this, symmetrize);
535 }
536 }
537 GO numDropped = lclA.nnz() - nnz_filtered;
538 GO numGlobalDropped;
539 Teuchos::reduceAll(*A->getRowMap()->getComm(), Teuchos::REDUCE_SUM, 1, &numDropped, &numGlobalDropped);
540 // We now know the number of entries of filtered A and have the final rowptr.
541
543 // Pass 4: Create local matrix for filtered A
544 //
545 // Dropped entries are optionally lumped to the diagonal.
546
547 RCP<Matrix> filteredA;
548 RCP<LWGraph_kokkos> graph;
549 if (numGlobalDropped > 0) {
550 SubFactoryMonitor mFill(*this, "Filtered matrix fill", currentLevel);
551
552 local_matrix_type lclFilteredA;
553 local_graph_type lclGraph;
554 if (reuseGraph) {
555 filteredA = MatrixFactory::BuildCopy(A);
556 lclFilteredA = filteredA->getLocalMatrixDevice();
557
558 auto colidx = entries_type("entries", nnz_filtered);
559 lclGraph = local_graph_type(colidx, filtered_rowptr);
560 } else {
561 auto colidx = entries_type("entries", nnz_filtered);
562 auto values = values_type("values", nnz_filtered);
563 lclFilteredA = local_matrix_type("filteredA",
564 lclA.numRows(), lclA.numCols(),
565 nnz_filtered,
566 values, filtered_rowptr, colidx);
567 }
568
569 if (lumpingChoice != MueLu::MatrixConstruction::no_lumping) {
570 if (reuseGraph) {
571 auto fillFunctor = MatrixConstruction::PointwiseFillReuseFunctor<local_matrix_type, local_graph_type, true>(lclA, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
572 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_reuse", range, fillFunctor);
573 } else {
574 if (lumpingChoice == MueLu::MatrixConstruction::diag_lumping) {
575 auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor<local_matrix_type, MueLu::MatrixConstruction::diag_lumping>(lclA, results, lclFilteredA, filteringDirichletThreshold);
576 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_noreuse", range, fillFunctor);
577 } else if (lumpingChoice == MueLu::MatrixConstruction::distributed_lumping) {
578 auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor<local_matrix_type, MueLu::MatrixConstruction::distributed_lumping>(lclA, results, lclFilteredA, filteringDirichletThreshold);
579 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_noreuse", range, fillFunctor);
580 }
581 }
582 } else {
583 if (reuseGraph) {
584 auto fillFunctor = MatrixConstruction::PointwiseFillReuseFunctor<local_matrix_type, local_graph_type, false>(lclA, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
585 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_reuse", range, fillFunctor);
586 } else {
587 auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor<local_matrix_type, MueLu::MatrixConstruction::no_lumping>(lclA, results, lclFilteredA, filteringDirichletThreshold);
588 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_noreuse", range, fillFunctor);
589 }
590 }
591
592 if (!reuseGraph)
593 filteredA = MatrixFactory::Build(lclFilteredA, A->getRowMap(), A->getColMap(), A->getDomainMap(), A->getRangeMap());
594 filteredA->SetFixedBlockSize(A->GetFixedBlockSize());
595
596 if (reuseEigenvalue) {
597 // Reuse max eigenvalue from A
598 // It is unclear what eigenvalue is the best for the smoothing, but we already may have
599 // the D^{-1}A estimate in A, may as well use it.
600 // NOTE: ML does that too
601 filteredA->SetMaxEigenvalueEstimate(A->GetMaxEigenvalueEstimate());
602 } else {
603 filteredA->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
604 }
605
606 if (!reuseGraph) {
607 // Use graph of filteredA as graph.
608 lclGraph = filteredA->getCrsGraph()->getLocalGraphDevice();
609 }
610 graph = rcp(new LWGraph_kokkos(lclGraph, filteredA->getRowMap(), filteredA->getColMap(), "amalgamated graph of A"));
611 graph->SetBoundaryNodeMap(boundaryNodes);
612 } else {
613 filteredA = A;
614 graph = rcp(new LWGraph_kokkos(filteredA->getCrsGraph()->getLocalGraphDevice(), filteredA->getRowMap(), filteredA->getColMap(), "amalgamated graph of A"));
615 graph->SetBoundaryNodeMap(boundaryNodes);
616 }
617
618 // Construct a second graph for coloring
619 if (generateColoringGraph) {
620 SubFactoryMonitor mColoringGraph(*this, "Construct coloring graph", currentLevel);
621
622 filtered_rowptr = rowptr_type("rowptr_coloring_graph", lclA.numRows() + 1);
623 if (localizeColoringGraph) {
624 auto drop_offrank = Misc::DropOffRankFunctor(lclA, results);
625 ScalarDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, results, filtered_rowptr, nnz_filtered, useBlocking, currentLevel, *this, drop_offrank);
626 }
627 if (symmetrizeColoringGraph) {
628 auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
629 ScalarDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, results, filtered_rowptr, nnz_filtered, useBlocking, currentLevel, *this, symmetrize);
630 }
631 auto colidx = entries_type("entries_coloring_graph", nnz_filtered);
632 auto lclGraph = local_graph_type(colidx, filtered_rowptr);
633 auto graphConstruction = MatrixConstruction::GraphConstruction<local_matrix_type, local_graph_type>(lclA, results, lclGraph);
634 Kokkos::parallel_for("MueLu::CoalesceDrop::Construct_coloring_graph", range, graphConstruction);
635
636 auto colorGraph = rcp(new LWGraph_kokkos(lclGraph, filteredA->getRowMap(), filteredA->getColMap(), "coloring graph"));
637 Set(currentLevel, "Coloring Graph", colorGraph);
638 }
639
640 if (pL.get<bool>("filtered matrix: count negative diagonals")) {
641 // Count the negative diagonals (and display that information)
643 GetOStream(Runtime0) << "CoalesceDrop: Negative diagonals: " << neg_count << std::endl;
644 }
645
646 LO dofsPerNode = 1;
647 Set(currentLevel, "DofsPerNode", dofsPerNode);
648 Set(currentLevel, "Graph", graph);
649 Set(currentLevel, "A", filteredA);
650
651 return std::make_tuple(numDropped, boundaryNodes);
652}
653
654template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
655std::tuple<GlobalOrdinal, typename MueLu::LWGraph_kokkos<LocalOrdinal, GlobalOrdinal, Node>::boundary_nodes_type> CoalesceDropFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
656 BuildVector(Level& currentLevel) const {
657 FactoryMonitor m(*this, "BuildVector", currentLevel);
658
659 using MatrixType = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
660 using GraphType = Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>;
661 using local_matrix_type = typename MatrixType::local_matrix_device_type;
662 using local_graph_type = typename GraphType::local_graph_device_type;
663 using rowptr_type = typename local_graph_type::row_map_type::non_const_type;
664 using entries_type = typename local_graph_type::entries_type::non_const_type;
665 using values_type = typename local_matrix_type::values_type::non_const_type;
666 using device_type = typename Node::device_type;
667 using memory_space = typename device_type::memory_space;
668 using results_view_type = Kokkos::View<DecisionType*, memory_space>;
669 using magnitudeType = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
670 using doubleMultiVector = Xpetra::MultiVector<magnitudeType, LO, GO, NO>;
671
672 typedef Teuchos::ScalarTraits<Scalar> STS;
673 const magnitudeType zero = Teuchos::ScalarTraits<magnitudeType>::zero();
674
675 auto A = Get<RCP<Matrix>>(currentLevel, "A");
676
677 /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
678 blkSize is the number of storage blocks that must kept together during the amalgamation process.
679
680 Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
681
682 numPDEs = blkSize * storageblocksize.
683
684 If numPDEs==1
685 Matrix is point storage (classical CRS storage). storageblocksize=1 and blkSize=1
686 No other values makes sense.
687
688 If numPDEs>1
689 If matrix uses point storage, then storageblocksize=1 and blkSize=numPDEs.
690 If matrix uses block storage, with block size of n, then storageblocksize=n, and blkSize=numPDEs/n.
691 Thus far, only storageblocksize=numPDEs and blkSize=1 has been tested.
692 */
693
694 TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0, Exceptions::RuntimeError, "A->GetFixedBlockSize() needs to be a multiple of A->GetStorageBlockSize()");
695 LO blkSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
696
697 auto amalInfo = Get<RCP<AmalgamationInfo>>(currentLevel, "UnAmalgamationInfo");
698
699 const RCP<const Map> rowMap = A->getRowMap();
700 const RCP<const Map> colMap = A->getColMap();
701
702 // build a node row map (uniqueMap = non-overlapping) and a node column map
703 // (nonUniqueMap = overlapping). The arrays rowTranslation and colTranslation
704 // stored in the AmalgamationInfo class container contain the local node id
705 // given a local dof id. The data is calculated in the AmalgamationFactory and
706 // stored in the variable "UnAmalgamationInfo" (which is of type AmalagamationInfo)
707 const RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
708 const RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
709 Array<LO> rowTranslationArray = *(amalInfo->getRowTranslation()); // TAW should be transform that into a View?
710 Array<LO> colTranslationArray = *(amalInfo->getColTranslation());
711
712 Kokkos::View<LO*, Kokkos::MemoryUnmanaged>
713 rowTranslationView(rowTranslationArray.getRawPtr(), rowTranslationArray.size());
714 Kokkos::View<LO*, Kokkos::MemoryUnmanaged>
715 colTranslationView(colTranslationArray.getRawPtr(), colTranslationArray.size());
716
717 // get number of local nodes
718 LO numNodes = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
719 typedef typename Kokkos::View<LocalOrdinal*, typename Node::device_type> id_translation_type;
720 id_translation_type rowTranslation("dofId2nodeId", rowTranslationArray.size());
721 id_translation_type colTranslation("ov_dofId2nodeId", colTranslationArray.size());
722 Kokkos::deep_copy(rowTranslation, rowTranslationView);
723 Kokkos::deep_copy(colTranslation, colTranslationView);
724
725 // extract striding information
726 blkSize = A->GetFixedBlockSize(); //< the full block size (number of dofs per node in strided map)
727 LocalOrdinal blkId = -1; //< the block id within a strided map or -1 if it is a full block map
728 LocalOrdinal blkPartSize = A->GetFixedBlockSize(); //< stores block size of part blkId (or the full block size)
729 if (A->IsView("stridedMaps") == true) {
730 const RCP<const Map> myMap = A->getRowMap("stridedMaps");
731 const RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
732 TEUCHOS_TEST_FOR_EXCEPTION(strMap.is_null() == true, Exceptions::RuntimeError, "Map is not of type stridedMap");
733 blkSize = Teuchos::as<const LocalOrdinal>(strMap->getFixedBlockSize());
734 blkId = strMap->getStridedBlockId();
735 if (blkId > -1)
736 blkPartSize = Teuchos::as<LocalOrdinal>(strMap->getStridingData()[blkId]);
737 }
738
739 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);
740
742 // Process parameterlist
743 const ParameterList& pL = GetParameterList();
744
745 // Boundary detection
746 const magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<double>("aggregation: Dirichlet threshold")));
747 const magnitudeType rowSumTol = as<magnitudeType>(pL.get<double>("aggregation: row sum drop tol"));
748 const LocalOrdinal dirichletNonzeroThreshold = 1;
749 const bool useGreedyDirichlet = pL.get<bool>("aggregation: greedy Dirichlet");
750 TEUCHOS_TEST_FOR_EXCEPTION(rowSumTol > zero, MueLu::Exceptions::RuntimeError, "MueLu::CoalesceDropFactory: RowSum is not implemented for vectorial problems.");
751
752 // Dropping
753 bool useBlocking = pL.get<bool>("aggregation: use blocking");
754 std::string droppingMethod = pL.get<std::string>("aggregation: drop scheme");
755 std::string socUsesMatrix = pL.get<std::string>("aggregation: strength-of-connection: matrix");
756 std::string socUsesMeasure = pL.get<std::string>("aggregation: strength-of-connection: measure");
757 std::string distanceLaplacianMetric = pL.get<std::string>("aggregation: distance laplacian metric");
758 bool symmetrizeDroppedGraph = pL.get<bool>("aggregation: symmetrize graph after dropping");
759 magnitudeType threshold;
760 // If we're doing the ML-style halving of the drop tol at each level, we do that here.
761 if (pL.get<bool>("aggregation: use ml scaling of drop tol"))
762 threshold = pL.get<double>("aggregation: drop tol") / pow(2.0, currentLevel.GetLevelID());
763 else
764 threshold = as<magnitudeType>(pL.get<double>("aggregation: drop tol"));
765 bool aggregationMayCreateDirichlet = pL.get<bool>("aggregation: dropping may create Dirichlet");
766
767 // Fill
768 const bool reuseGraph = pL.get<bool>("filtered matrix: reuse graph");
769 const bool reuseEigenvalue = pL.get<bool>("filtered matrix: reuse eigenvalue");
770
771 const bool useRootStencil = pL.get<bool>("filtered matrix: use root stencil");
772 const bool useSpreadLumping = pL.get<bool>("filtered matrix: use spread lumping");
773 const std::string lumpingChoiceString = pL.get<std::string>("filtered matrix: lumping choice");
775 if (lumpingChoiceString == "diag lumping")
777 else if (lumpingChoiceString == "distributed lumping")
779
780 const magnitudeType filteringDirichletThreshold = as<magnitudeType>(pL.get<double>("filtered matrix: Dirichlet threshold"));
781
782 // coloring graph
783 bool generateColoringGraph = pL.get<bool>("aggregation: coloring: use color graph");
784 const bool localizeColoringGraph = pL.get<bool>("aggregation: coloring: localize color graph");
785 const bool symmetrizeColoringGraph = true;
786
787#ifdef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
788 translateOldAlgoParam(pL, droppingMethod, useBlocking, socUsesMatrix, socUsesMeasure, symmetrizeDroppedGraph, generateColoringGraph, threshold, lumpingChoice);
789#endif
790 {
791 std::stringstream ss;
792 ss << "dropping scheme = \"" << droppingMethod << "\", strength-of-connection measure = \"" << socUsesMeasure << "\", strength-of-connection matrix = \"" << socUsesMatrix << "\", ";
793 if (socUsesMatrix == "distance laplacian")
794 ss << "distance laplacian metric = \"" << distanceLaplacianMetric << "\", ";
795 ss << "threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << ", useBlocking = " << useBlocking;
796 ss << ", symmetrizeDroppedGraph = " << symmetrizeDroppedGraph << std::endl;
797
798 GetOStream(Runtime0) << ss.str();
799 }
800
801 TEUCHOS_ASSERT(!useRootStencil);
802 TEUCHOS_ASSERT(!useSpreadLumping);
803 TEUCHOS_ASSERT((lumpingChoice != MueLu::MatrixConstruction::distributed_lumping) || !reuseGraph);
804 if (droppingMethod == "cut-drop")
805 TEUCHOS_TEST_FOR_EXCEPTION(threshold > 1.0, Exceptions::RuntimeError, "For cut-drop algorithms, \"aggregation: drop tol\" = " << threshold << ", needs to be <= 1.0");
806
808 // We perform four sweeps over the rows of A:
809 // Pass 1: detection of boundary nodes
810 // Pass 2: diagonal extraction
811 // Pass 3: drop decision for each entry and construction of the rowptr of the filtered matrix
812 // Pass 4: fill of the filtered matrix
813 //
814 // Pass 1 and 3 apply a sequence of criteria to each row of the matrix.
815
816 // TODO: We could merge pass 1 and 2.
817
818 auto crsA = toCrsMatrix(A);
819 auto lclA = crsA->getLocalMatrixDevice();
820 auto range = range_type(0, numNodes);
821
823 // Pass 1: Detect boundary nodes
824 //
825 // The following criteria are available:
826 // - BoundaryDetection::VectorDirichletFunctor
827 // Marks rows as Dirichlet based on value threshold and number of off-diagonal entries
828
829 // Dirichlet nodes
830 auto boundaryNodes = boundary_nodes_type("boundaryNodes", numNodes); // initialized to false
831 {
832 SubFactoryMonitor mBoundary(*this, "Boundary detection", currentLevel);
833
834 if (useGreedyDirichlet) {
835 auto dirichlet_detection = BoundaryDetection::VectorDirichletFunctor<local_matrix_type, true>(lclA, blkPartSize, boundaryNodes, dirichletThreshold, dirichletNonzeroThreshold);
836 runBoundaryFunctors(lclA, boundaryNodes, dirichlet_detection);
837 } else {
838 auto dirichlet_detection = BoundaryDetection::VectorDirichletFunctor<local_matrix_type, false>(lclA, blkPartSize, boundaryNodes, dirichletThreshold, dirichletNonzeroThreshold);
839 runBoundaryFunctors(lclA, boundaryNodes, dirichlet_detection);
840 }
841 }
842 // In what follows, boundaryNodes can still still get modified if aggregationMayCreateDirichlet == true.
843 // Otherwise we're now done with it now.
844
846 // Pass 2 & 3: Diagonal extraction and determine dropping and construct
847 // rowptr of filtered matrix
848 //
849 // The following criteria are available:
850 // - Misc::VectorDropBoundaryFunctor
851 // Drop all rows that have been marked as Dirichlet
852 // - Misc::DropOffRankFunctor
853 // Drop all entries that are off-rank
854 // - ClassicalDropping::DropFunctor
855 // Classical dropping
856 // - DistanceLaplacian::VectorDropFunctor
857 // Distance Laplacian dropping
858 // - Misc::KeepDiagonalFunctor
859 // Mark diagonal as KEEP
860 // - Misc::MarkSingletonFunctor
861 // Mark singletons after dropping as Dirichlet
862
863 // rowptr of filtered A
864 auto filtered_rowptr = rowptr_type("rowptr", lclA.numRows() + 1);
865 auto graph_rowptr = rowptr_type("rowptr", numNodes + 1);
866 // Number of nonzeros of filtered A and graph
867 Kokkos::pair<LocalOrdinal, LocalOrdinal> nnz = {0, 0};
868
869 // dropping decisions for each entry
870 auto results = results_view_type("results", lclA.nnz()); // initialized to UNDECIDED
871
872 RCP<Matrix> mergedA;
873 {
874 SubFactoryMonitor mDropping(*this, "Dropping decisions", currentLevel);
875
876 {
877 // Construct merged A.
878
879 auto merged_rowptr = rowptr_type("rowptr", numNodes + 1);
880 LocalOrdinal nnz_merged = 0;
881
882 auto functor = MatrixConstruction::MergeCountFunctor(lclA, blkPartSize, colTranslation, merged_rowptr);
883 Kokkos::parallel_scan("MergeCount", range, functor, nnz_merged);
884
885 local_graph_type lclMergedGraph;
886 auto colidx_merged = entries_type("entries", nnz_merged);
887 auto values_merged = values_type("values", nnz_merged);
888
889 local_matrix_type lclMergedA = local_matrix_type("mergedA",
890 numNodes, nonUniqueMap->getLocalNumElements(),
891 nnz_merged,
892 values_merged, merged_rowptr, colidx_merged);
893
894 auto fillFunctor = MatrixConstruction::MergeFillFunctor<local_matrix_type>(lclA, blkSize, colTranslation, lclMergedA);
895 Kokkos::parallel_for("MueLu::CoalesceDrop::MergeFill", range, fillFunctor);
896
897 mergedA = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(lclMergedA, uniqueMap, nonUniqueMap, uniqueMap, uniqueMap);
898 }
899
900 if (threshold != zero) {
901 if (socUsesMatrix == "A") {
902 if (socUsesMeasure == "unscaled") {
903 VectorDroppingClassical<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::UnscaledMeasure>::runDroppingFunctors_on_A(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, currentLevel, *this);
904 } else if (socUsesMeasure == "smoothed aggregation") {
905 VectorDroppingClassical<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SmoothedAggregationMeasure>::runDroppingFunctors_on_A(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, currentLevel, *this);
906 } else if (socUsesMeasure == "signed ruge-stueben") {
907 VectorDroppingClassical<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SignedRugeStuebenMeasure>::runDroppingFunctors_on_A(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, currentLevel, *this);
908 } else if (socUsesMeasure == "signed smoothed aggregation") {
909 VectorDroppingClassical<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SignedSmoothedAggregationMeasure>::runDroppingFunctors_on_A(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, currentLevel, *this);
910 }
911 } else if (socUsesMatrix == "distance laplacian") {
912 auto coords = Get<RCP<doubleMultiVector>>(currentLevel, "Coordinates");
913
914 Array<double> dlap_weights = pL.get<Array<double>>("aggregation: distance laplacian directional weights");
915 LocalOrdinal interleaved_blocksize = as<LocalOrdinal>(pL.get<int>("aggregation: block diagonal: interleaved blocksize"));
916 if (socUsesMeasure == "distance laplacian") {
917 LO dim = (LO)coords->getNumVectors();
918 // If anything isn't 1.0 we need to turn on the weighting
919 bool non_unity = false;
920 for (LO i = 0; !non_unity && i < (LO)dlap_weights.size(); i++) {
921 if (dlap_weights[i] != 1.0) {
922 non_unity = true;
923 }
924 }
925 if (non_unity) {
926 if ((LO)dlap_weights.size() == dim) {
927 distanceLaplacianMetric = "weighted";
928 } else if ((LO)dlap_weights.size() == interleaved_blocksize * dim)
929 distanceLaplacianMetric = "block weighted";
930 else {
931 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError,
932 "length of 'aggregation: distance laplacian directional weights' must equal the coordinate dimension OR the coordinate dimension times the blocksize");
933 }
934 if (GetVerbLevel() & Statistics1)
935 GetOStream(Statistics1) << "Using distance laplacian weights: " << dlap_weights << std::endl;
936 }
937 }
938
939 if (socUsesMeasure == "unscaled") {
940 VectorDroppingDistanceLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::UnscaledMeasure>::runDroppingFunctors_on_dlap(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, distanceLaplacianMetric, dlap_weights, interleaved_blocksize, currentLevel, *this);
941 } else if (socUsesMeasure == "smoothed aggregation") {
942 VectorDroppingDistanceLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SmoothedAggregationMeasure>::runDroppingFunctors_on_dlap(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, distanceLaplacianMetric, dlap_weights, interleaved_blocksize, currentLevel, *this);
943 } else if (socUsesMeasure == "signed ruge-stueben") {
944 VectorDroppingDistanceLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SignedRugeStuebenMeasure>::runDroppingFunctors_on_dlap(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, distanceLaplacianMetric, dlap_weights, interleaved_blocksize, currentLevel, *this);
945 } else if (socUsesMeasure == "signed smoothed aggregation") {
946 VectorDroppingDistanceLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node, Misc::SignedSmoothedAggregationMeasure>::runDroppingFunctors_on_dlap(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, boundaryNodes, droppingMethod, threshold, aggregationMayCreateDirichlet, symmetrizeDroppedGraph, useBlocking, distanceLaplacianMetric, dlap_weights, interleaved_blocksize, currentLevel, *this);
947 }
948 }
949 } else {
950 Kokkos::deep_copy(results, KEEP);
951
952 auto no_op = Misc::NoOpFunctor<LocalOrdinal>();
953 VectorDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, useBlocking, currentLevel, *this, no_op);
954 }
955
956 if (symmetrizeDroppedGraph) {
957 auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
958 VectorDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, useBlocking, currentLevel, *this, symmetrize);
959 }
960 }
961 LocalOrdinal nnz_filtered = nnz.first;
962 LocalOrdinal nnz_graph = nnz.second;
963 GO numTotal = lclA.nnz();
964 GO numDropped = numTotal - nnz_filtered;
965 GO numGlobalDropped;
966 Teuchos::reduceAll(*A->getRowMap()->getComm(), Teuchos::REDUCE_SUM, 1, &numDropped, &numGlobalDropped);
967 // We now know the number of entries of filtered A and have the final rowptr.
968
970 // Pass 4: Create local matrix for filtered A
971 //
972 // Dropped entries are optionally lumped to the diagonal.
973
974 RCP<Matrix> filteredA;
975 RCP<LWGraph_kokkos> graph;
976 if (numGlobalDropped > 0) {
977 SubFactoryMonitor mFill(*this, "Filtered matrix fill", currentLevel);
978
979 local_matrix_type lclFilteredA;
980 if (reuseGraph) {
981 lclFilteredA = local_matrix_type("filteredA", lclA.graph, lclA.numCols());
982 } else {
983 auto colidx = entries_type("entries", nnz_filtered);
984 auto values = values_type("values", nnz_filtered);
985 lclFilteredA = local_matrix_type("filteredA",
986 lclA.numRows(), lclA.numCols(),
987 nnz_filtered,
988 values, filtered_rowptr, colidx);
989 }
990
991 local_graph_type lclGraph;
992 {
993 auto colidx = entries_type("entries", nnz_graph);
994 lclGraph = local_graph_type(colidx, graph_rowptr);
995 }
996
997 if (lumpingChoice != MueLu::MatrixConstruction::no_lumping) {
998 if (reuseGraph) {
999 auto fillFunctor = MatrixConstruction::VectorFillFunctor<local_matrix_type, true, true>(lclA, blkPartSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1000 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_reuse", range, fillFunctor);
1001 } else {
1002 auto fillFunctor = MatrixConstruction::VectorFillFunctor<local_matrix_type, true, false>(lclA, blkPartSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1003 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_noreuse", range, fillFunctor);
1004 }
1005 } else {
1006 if (reuseGraph) {
1007 auto fillFunctor = MatrixConstruction::VectorFillFunctor<local_matrix_type, false, true>(lclA, blkSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1008 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_reuse", range, fillFunctor);
1009 } else {
1010 auto fillFunctor = MatrixConstruction::VectorFillFunctor<local_matrix_type, false, false>(lclA, blkSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold);
1011 Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_noreuse", range, fillFunctor);
1012 }
1013 }
1014
1015 filteredA = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(lclFilteredA, A->getRowMap(), A->getColMap(), A->getDomainMap(), A->getRangeMap());
1016 filteredA->SetFixedBlockSize(blkSize);
1017
1018 if (reuseEigenvalue) {
1019 // Reuse max eigenvalue from A
1020 // It is unclear what eigenvalue is the best for the smoothing, but we already may have
1021 // the D^{-1}A estimate in A, may as well use it.
1022 // NOTE: ML does that too
1023 filteredA->SetMaxEigenvalueEstimate(A->GetMaxEigenvalueEstimate());
1024 } else {
1025 filteredA->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
1026 }
1027
1028 graph = rcp(new LWGraph_kokkos(lclGraph, uniqueMap, nonUniqueMap, "amalgamated graph of A"));
1029 graph->SetBoundaryNodeMap(boundaryNodes);
1030 } else {
1031 filteredA = A;
1032 graph = rcp(new LWGraph_kokkos(mergedA->getCrsGraph()->getLocalGraphDevice(), uniqueMap, nonUniqueMap, "amalgamated graph of A"));
1033 graph->SetBoundaryNodeMap(boundaryNodes);
1034 }
1035
1036 // Construct a second graph for coloring
1037 if (generateColoringGraph) {
1038 SubFactoryMonitor mColoringGraph(*this, "Construct coloring graph", currentLevel);
1039
1040 filtered_rowptr = rowptr_type("rowptr_coloring_graph", lclA.numRows() + 1);
1041 graph_rowptr = rowptr_type("rowptr", numNodes + 1);
1042 if (localizeColoringGraph) {
1043 auto drop_offrank = Misc::DropOffRankFunctor(lclA, results);
1044 VectorDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, useBlocking, currentLevel, *this, drop_offrank);
1045 }
1046 if (symmetrizeColoringGraph) {
1047 auto symmetrize = Misc::SymmetrizeFunctor(lclA, results);
1048 VectorDroppingBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::template runDroppingFunctors<>(*A, *mergedA, blkPartSize, rowTranslation, colTranslation, results, filtered_rowptr, graph_rowptr, nnz, useBlocking, currentLevel, *this, symmetrize);
1049 }
1050 auto colidx = entries_type("entries_coloring_graph", nnz_filtered);
1051 auto lclGraph = local_graph_type(colidx, filtered_rowptr);
1052 auto graphConstruction = MatrixConstruction::GraphConstruction<local_matrix_type, local_graph_type>(lclA, results, lclGraph);
1053 Kokkos::parallel_for("MueLu::CoalesceDrop::Construct_coloring_graph", range, graphConstruction);
1054
1055 auto colorGraph = rcp(new LWGraph_kokkos(lclGraph, filteredA->getRowMap(), filteredA->getColMap(), "coloring graph"));
1056 Set(currentLevel, "Coloring Graph", colorGraph);
1057 }
1058
1059 LO dofsPerNode = blkSize;
1060
1061 Set(currentLevel, "DofsPerNode", dofsPerNode);
1062 Set(currentLevel, "Graph", graph);
1063 Set(currentLevel, "A", filteredA);
1064
1065 return std::make_tuple(numDropped, boundaryNodes);
1066}
1067
1068} // namespace MueLu
1069#endif // MUELU_COALESCEDROPFACTORY_KOKKOS_DEF_HPP
#define SET_VALID_ENTRY(name)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Functor that serially applies sub-functors to rows.
Functor for marking nodes as Dirichlet based on rowsum.
Functor for marking nodes as Dirichlet in a block operator.
void DeclareInput(Level &currentLevel) const
Input.
Kokkos::RangePolicy< local_ordinal_type, execution_space > range_type
void Build(Level &currentLevel) const
Build an object with this factory.
typename MueLu::LWGraph_kokkos< LocalOrdinal, GlobalOrdinal, Node >::boundary_nodes_type boundary_nodes_type
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
std::tuple< GlobalOrdinal, boundary_nodes_type > BuildVector(Level &currentLevel) const
std::tuple< GlobalOrdinal, boundary_nodes_type > BuildScalar(Level &currentLevel) const
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Lightweight MueLu representation of a compressed row storage graph.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Functor does not reuse the graph of the matrix for a problem with blockSize == 1.
Functor that fills the filtered matrix while reusing the graph of the matrix before dropping,...
Functor that drops off-rank entries.
Functor that drops boundary nodes for a blockSize == 1 problem.
Functor that symmetrizes the dropping decisions.
static const NoFactory * get()
static void runDroppingFunctors_on_A(matrix_type &A, results_view &results, rowptr_type &filtered_rowptr, LocalOrdinal &nnz_filtered, boundary_nodes_type &boundaryNodes, const std::string &droppingMethod, const magnitudeType threshold, const bool aggregationMayCreateDirichlet, const bool symmetrizeDroppedGraph, const bool useBlocking, Level &level, const Factory &factory)
static void runDroppingFunctors_on_dlap(matrix_type &A, results_view &results, rowptr_type &filtered_rowptr, LocalOrdinal &nnz_filtered, boundary_nodes_type &boundaryNodes, const std::string &droppingMethod, const magnitudeType threshold, const bool aggregationMayCreateDirichlet, const bool symmetrizeDroppedGraph, const bool useBlocking, const std::string &distanceLaplacianMetric, Level &level, const Factory &factory)
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > SPAI(const RCP< Matrix > &original)
Creates a sparse approximate inverse of a matrix with the same nonzero pattern as the input matrix.
static GlobalOrdinal CountNegativeDiagonalEntries(const Matrix &A)
Counts the number of negative diagonal entries.
static void runDroppingFunctors_on_A(matrix_type &A, matrix_type &mergedA, LocalOrdinal blkPartSize, block_indices_view_type &rowTranslation, block_indices_view_type &colTranslation, results_view &results, rowptr_type &filtered_rowptr, rowptr_type &graph_rowptr, nnz_count_type &nnz, boundary_nodes_type &boundaryNodes, const std::string &droppingMethod, const magnitudeType threshold, const bool aggregationMayCreateDirichlet, const bool symmetrizeDroppedGraph, const bool useBlocking, Level &level, const Factory &factory)
static void runDroppingFunctors_on_dlap(matrix_type &A, matrix_type &mergedA, LocalOrdinal blkPartSize, block_indices_view_type &rowTranslation, block_indices_view_type &colTranslation, results_view &results, rowptr_type &filtered_rowptr, rowptr_type &graph_rowptr, nnz_count_type &nnz, boundary_nodes_type &boundaryNodes, const std::string &droppingMethod, const magnitudeType threshold, const bool aggregationMayCreateDirichlet, const bool symmetrizeDroppedGraph, const bool useBlocking, const std::string &distanceLaplacianMetric, Teuchos::Array< double > &dlap_weights, LocalOrdinal interleaved_blocksize, Level &level, const Factory &factory)
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.
@ Runtime0
One-liner description of what is happening.
void runBoundaryFunctors(local_matrix_type &lclA, boundary_nodes_view &boundaryNodes, Functors &... functors)
void translateOldAlgoParam(const Teuchos::ParameterList &pL, std::string &droppingMethod, bool &useBlocking, std::string &socUsesMatrix, std::string &socUsesMeasure, bool &symmetrizeDroppedGraph, bool &generateColoringGraph, magnitudeType &threshold, MueLu::MatrixConstruction::lumpingType &lumpingChoice)