MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_CoalesceDropFactory_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_DEF_HPP
11#define MUELU_COALESCEDROPFACTORY_DEF_HPP
12
13#include <Xpetra_CrsGraphFactory.hpp>
14#include <Xpetra_CrsGraph.hpp>
15#include <Xpetra_ImportFactory.hpp>
16#include <Xpetra_ExportFactory.hpp>
17#include <Xpetra_MapFactory.hpp>
18#include <Xpetra_Map.hpp>
19#include <Xpetra_Matrix.hpp>
20#include <Xpetra_MultiVectorFactory.hpp>
21#include <Xpetra_MultiVector.hpp>
22#include <Xpetra_StridedMap.hpp>
23#include <Xpetra_VectorFactory.hpp>
24#include <Xpetra_Vector.hpp>
25
26#include <Xpetra_IO.hpp>
27
28#include <Kokkos_NestedSort.hpp>
29#include <Kokkos_StdAlgorithms.hpp>
30#include "MueLu_CoalesceDropFactory_decl.hpp"
31
32#include "MueLu_AmalgamationFactory.hpp"
33#include "MueLu_AmalgamationInfo.hpp"
34#include "MueLu_Exceptions.hpp"
35#include "MueLu_LWGraph.hpp"
36
37#include "MueLu_Level.hpp"
38#include "MueLu_MasterList.hpp"
39#include "MueLu_Monitor.hpp"
40#include "MueLu_PreDropFunctionConstVal.hpp"
41#include "MueLu_Utilities.hpp"
42
43#include "Tpetra_CrsGraphTransposer.hpp"
44
45#include <algorithm>
46#include <cstdlib>
47#include <string>
48
49// If defined, read environment variables.
50// Should be removed once we are confident that this works.
51//#define DJS_READ_ENV_VARIABLES
52
53namespace MueLu {
54
55namespace Details {
56template <class real_type, class LO>
57struct DropTol {
58 DropTol() = default;
59 DropTol(DropTol const&) = default;
60 DropTol(DropTol&&) = default;
61
62 DropTol& operator=(DropTol const&) = default;
63 DropTol& operator=(DropTol&&) = default;
64
65 DropTol(real_type val_, real_type diag_, LO col_, bool drop_)
66 : val{val_}
67 , diag{diag_}
68 , col{col_}
69 , drop{drop_} {}
70
71 real_type val{Teuchos::ScalarTraits<real_type>::zero()};
72 real_type diag{Teuchos::ScalarTraits<real_type>::zero()};
73 LO col{Teuchos::OrdinalTraits<LO>::invalid()};
74 bool drop{true};
75
76 // CMS: Auxillary information for debugging info
77 // real_type aux_val {Teuchos::ScalarTraits<real_type>::nan()};
78};
79} // namespace Details
80
85
86template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
87RCP<const ParameterList> CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
88 RCP<ParameterList> validParamList = rcp(new ParameterList());
89
90#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
91 SET_VALID_ENTRY("aggregation: drop tol");
92 SET_VALID_ENTRY("aggregation: use ml scaling of drop tol");
93 SET_VALID_ENTRY("aggregation: Dirichlet threshold");
94 SET_VALID_ENTRY("aggregation: greedy Dirichlet");
95 SET_VALID_ENTRY("aggregation: row sum drop tol");
96 SET_VALID_ENTRY("aggregation: drop scheme");
97 SET_VALID_ENTRY("aggregation: block diagonal: interleaved blocksize");
98 SET_VALID_ENTRY("aggregation: distance laplacian directional weights");
99 SET_VALID_ENTRY("aggregation: dropping may create Dirichlet");
100
101 {
102 // "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)
103 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"))));
104 }
105 SET_VALID_ENTRY("aggregation: distance laplacian algo");
106 SET_VALID_ENTRY("aggregation: classical algo");
107 SET_VALID_ENTRY("aggregation: coloring: localize color graph");
108#undef SET_VALID_ENTRY
109 validParamList->set<bool>("lightweight wrap", true, "Experimental option for lightweight graph access");
110
111 validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory of the matrix A");
112 validParamList->set<RCP<const FactoryBase>>("UnAmalgamationInfo", Teuchos::null, "Generating factory for UnAmalgamationInfo");
113 validParamList->set<RCP<const FactoryBase>>("Coordinates", Teuchos::null, "Generating factory for Coordinates");
114 validParamList->set<RCP<const FactoryBase>>("BlockNumber", Teuchos::null, "Generating factory for BlockNUmber");
115
116 return validParamList;
117}
118
119template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
122
123template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
124void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
125 Input(currentLevel, "A");
126 Input(currentLevel, "UnAmalgamationInfo");
127
128 const ParameterList& pL = GetParameterList();
129 if (pL.get<bool>("lightweight wrap") == true) {
130 std::string algo = pL.get<std::string>("aggregation: drop scheme");
131 if (algo == "distance laplacian" || algo == "block diagonal distance laplacian") {
132 Input(currentLevel, "Coordinates");
133 }
134 if (algo == "signed classical sa")
135 ;
136 else if (algo.find("block diagonal") != std::string::npos || algo.find("signed classical") != std::string::npos) {
137 Input(currentLevel, "BlockNumber");
138 }
139 }
140}
141
142template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
143void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& currentLevel) const {
144 FactoryMonitor m(*this, "Build", currentLevel);
145
146 typedef Teuchos::ScalarTraits<SC> STS;
147 typedef typename STS::magnitudeType real_type;
148 typedef Xpetra::MultiVector<real_type, LO, GO, NO> RealValuedMultiVector;
149 typedef Xpetra::MultiVectorFactory<real_type, LO, GO, NO> RealValuedMultiVectorFactory;
150
151 if (predrop_ != Teuchos::null)
152 GetOStream(Parameters0) << predrop_->description();
153
154 RCP<Matrix> realA = Get<RCP<Matrix>>(currentLevel, "A");
155 RCP<AmalgamationInfo> amalInfo = Get<RCP<AmalgamationInfo>>(currentLevel, "UnAmalgamationInfo");
156 const ParameterList& pL = GetParameterList();
157 bool doExperimentalWrap = pL.get<bool>("lightweight wrap");
158
159 GetOStream(Parameters0) << "lightweight wrap = " << doExperimentalWrap << std::endl;
160 std::string algo = pL.get<std::string>("aggregation: drop scheme");
161 const bool aggregationMayCreateDirichlet = pL.get<bool>("aggregation: dropping may create Dirichlet");
162
163 RCP<RealValuedMultiVector> Coords;
164 RCP<Matrix> A;
165
166 bool use_block_algorithm = false;
167 LO interleaved_blocksize = as<LO>(pL.get<int>("aggregation: block diagonal: interleaved blocksize"));
168 bool useSignedClassicalRS = false;
169 bool useSignedClassicalSA = false;
170 bool generateColoringGraph = false;
171
172 // NOTE: If we're doing blockDiagonal, we'll not want to do rowSum twice (we'll do it
173 // in the block diagonalization). So we'll clobber the rowSumTol with -1.0 in this case
174 typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<double>("aggregation: row sum drop tol"));
175
176 RCP<LocalOrdinalVector> ghostedBlockNumber;
177
178 if (algo == "distance laplacian") {
179 // Grab the coordinates for distance laplacian
180 Coords = Get<RCP<RealValuedMultiVector>>(currentLevel, "Coordinates");
181 A = realA;
182 } else if (algo == "signed classical sa") {
183 useSignedClassicalSA = true;
184 algo = "classical";
185 A = realA;
186 } else if (algo == "signed classical" || algo == "block diagonal colored signed classical" || algo == "block diagonal signed classical") {
187 useSignedClassicalRS = true;
188 // if(realA->GetFixedBlockSize() > 1) {
189 RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
190 // Ghost the column block numbers if we need to
191 RCP<const Import> importer = realA->getCrsGraph()->getImporter();
192 if (!importer.is_null()) {
193 SubFactoryMonitor m1(*this, "Block Number import", currentLevel);
194 ghostedBlockNumber = Xpetra::VectorFactory<LO, LO, GO, NO>::Build(importer->getTargetMap());
195 ghostedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
196 } else {
197 ghostedBlockNumber = BlockNumber;
198 }
199 // }
200 if (algo == "block diagonal colored signed classical")
201 generateColoringGraph = true;
202 algo = "classical";
203 A = realA;
204
205 } else if (algo == "block diagonal") {
206 // Handle the "block diagonal" filtering and then leave
207 BlockDiagonalize(currentLevel, realA, false);
208 return;
209 } else if (algo == "block diagonal classical" || algo == "block diagonal distance laplacian") {
210 // Handle the "block diagonal" filtering, and then continue onward
211 use_block_algorithm = true;
212 RCP<Matrix> filteredMatrix = BlockDiagonalize(currentLevel, realA, true);
213 if (algo == "block diagonal distance laplacian") {
214 // We now need to expand the coordinates by the interleaved blocksize
215 RCP<RealValuedMultiVector> OldCoords = Get<RCP<RealValuedMultiVector>>(currentLevel, "Coordinates");
216 if (OldCoords->getLocalLength() != realA->getLocalNumRows()) {
217 LO dim = (LO)OldCoords->getNumVectors();
218 Coords = RealValuedMultiVectorFactory::Build(realA->getRowMap(), dim);
219 for (LO k = 0; k < dim; k++) {
220 ArrayRCP<const real_type> old_vec = OldCoords->getData(k);
221 ArrayRCP<real_type> new_vec = Coords->getDataNonConst(k);
222 for (LO i = 0; i < (LO)OldCoords->getLocalLength(); i++) {
223 LO new_base = i * dim;
224 for (LO j = 0; j < interleaved_blocksize; j++)
225 new_vec[new_base + j] = old_vec[i];
226 }
227 }
228 } else {
229 Coords = OldCoords;
230 }
231 algo = "distance laplacian";
232 } else if (algo == "block diagonal classical") {
233 algo = "classical";
234 }
235 // All cases
236 A = filteredMatrix;
237 rowSumTol = -1.0;
238 } else {
239 A = realA;
240 }
241
242 // Distance Laplacian weights
243 Array<double> dlap_weights = pL.get<Array<double>>("aggregation: distance laplacian directional weights");
244 enum { NO_WEIGHTS = 0,
245 SINGLE_WEIGHTS,
246 BLOCK_WEIGHTS };
247 int use_dlap_weights = NO_WEIGHTS;
248 if (algo == "distance laplacian") {
249 LO dim = (LO)Coords->getNumVectors();
250 // If anything isn't 1.0 we need to turn on the weighting
251 bool non_unity = false;
252 for (LO i = 0; !non_unity && i < (LO)dlap_weights.size(); i++) {
253 if (dlap_weights[i] != 1.0) {
254 non_unity = true;
255 }
256 }
257 if (non_unity) {
258 LO blocksize = use_block_algorithm ? as<LO>(pL.get<int>("aggregation: block diagonal: interleaved blocksize")) : 1;
259 if ((LO)dlap_weights.size() == dim)
260 use_dlap_weights = SINGLE_WEIGHTS;
261 else if ((LO)dlap_weights.size() == blocksize * dim)
262 use_dlap_weights = BLOCK_WEIGHTS;
263 else {
264 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError,
265 "length of 'aggregation: distance laplacian directional weights' must equal the coordinate dimension OR the coordinate dimension times the blocksize");
266 }
267 if (GetVerbLevel() & Statistics1)
268 GetOStream(Statistics1) << "Using distance laplacian weights: " << dlap_weights << std::endl;
269 }
270 }
271
272 // decide wether to use the fast-track code path for standard maps or the somewhat slower
273 // code path for non-standard maps
274 /*bool bNonStandardMaps = false;
275 if (A->IsView("stridedMaps") == true) {
276 Teuchos::RCP<const Map> myMap = A->getRowMap("stridedMaps");
277 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
278 TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
279 if (strMap->getStridedBlockId() != -1 || strMap->getOffset() > 0)
280 bNonStandardMaps = true;
281 }*/
282
283 if (doExperimentalWrap) {
284 TEUCHOS_TEST_FOR_EXCEPTION(predrop_ != null && algo != "classical", Exceptions::RuntimeError, "Dropping function must not be provided for \"" << algo << "\" algorithm");
285 TEUCHOS_TEST_FOR_EXCEPTION(algo != "classical" && algo != "distance laplacian" && algo != "signed classical", Exceptions::RuntimeError, "\"algorithm\" must be one of (classical|distance laplacian|signed classical)");
286
287 SC threshold;
288 // If we're doing the ML-style halving of the drop tol at each level, we do that here.
289 if (pL.get<bool>("aggregation: use ml scaling of drop tol"))
290 threshold = pL.get<double>("aggregation: drop tol") / pow(2.0, currentLevel.GetLevelID());
291 else
292 threshold = as<SC>(pL.get<double>("aggregation: drop tol"));
293
294 std::string distanceLaplacianAlgoStr = pL.get<std::string>("aggregation: distance laplacian algo");
295 std::string classicalAlgoStr = pL.get<std::string>("aggregation: classical algo");
296 real_type realThreshold = STS::magnitude(threshold); // CMS: Rename this to "magnitude threshold" sometime
297
299 // Remove this bit once we are confident that cut-based dropping works.
300#ifdef HAVE_MUELU_DEBUG
301 int distanceLaplacianCutVerbose = 0;
302#endif
303#ifdef DJS_READ_ENV_VARIABLES
304 if (getenv("MUELU_DROP_TOLERANCE_MODE")) {
305 distanceLaplacianAlgoStr = std::string(getenv("MUELU_DROP_TOLERANCE_MODE"));
306 }
307
308 if (getenv("MUELU_DROP_TOLERANCE_THRESHOLD")) {
309 auto tmp = atoi(getenv("MUELU_DROP_TOLERANCE_THRESHOLD"));
310 realThreshold = 1e-4 * tmp;
311 }
312
313#ifdef HAVE_MUELU_DEBUG
314 if (getenv("MUELU_DROP_TOLERANCE_VERBOSE")) {
315 distanceLaplacianCutVerbose = atoi(getenv("MUELU_DROP_TOLERANCE_VERBOSE"));
316 }
317#endif
318#endif
320
321 decisionAlgoType distanceLaplacianAlgo = defaultAlgo;
322 decisionAlgoType classicalAlgo = defaultAlgo;
323 if (algo == "distance laplacian") {
324 if (distanceLaplacianAlgoStr == "default")
325 distanceLaplacianAlgo = defaultAlgo;
326 else if (distanceLaplacianAlgoStr == "unscaled cut")
327 distanceLaplacianAlgo = unscaled_cut;
328 else if (distanceLaplacianAlgoStr == "scaled cut")
329 distanceLaplacianAlgo = scaled_cut;
330 else if (distanceLaplacianAlgoStr == "scaled cut symmetric")
331 distanceLaplacianAlgo = scaled_cut_symmetric;
332 else
333 TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut), not \"" << distanceLaplacianAlgoStr << "\"");
334 GetOStream(Runtime0) << "algorithm = \"" << algo << "\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
335 } else if (algo == "classical") {
336 if (classicalAlgoStr == "default")
337 classicalAlgo = defaultAlgo;
338 else if (classicalAlgoStr == "unscaled cut")
339 classicalAlgo = unscaled_cut;
340 else if (classicalAlgoStr == "scaled cut")
341 classicalAlgo = scaled_cut;
342 else
343 TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: classical algo\" must be one of (default|unscaled cut|scaled cut), not \"" << classicalAlgoStr << "\"");
344 GetOStream(Runtime0) << "algorithm = \"" << algo << "\" classical algorithm = \"" << classicalAlgoStr << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
345
346 } else
347 GetOStream(Runtime0) << "algorithm = \"" << algo << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
348
349 if (((algo == "classical") && (classicalAlgoStr.find("scaled") != std::string::npos)) || ((algo == "distance laplacian") && (distanceLaplacianAlgoStr.find("scaled") != std::string::npos)))
350 TEUCHOS_TEST_FOR_EXCEPTION(realThreshold > 1.0, Exceptions::RuntimeError, "For cut-drop algorithms, \"aggregation: drop tol\" = " << threshold << ", needs to be <= 1.0");
351
352 Set<bool>(currentLevel, "Filtering", (threshold != STS::zero()));
353
354 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<double>("aggregation: Dirichlet threshold")));
355
356 // NOTE: We don't support signed classical RS or SA with cut drop at present
357 TEUCHOS_TEST_FOR_EXCEPTION(useSignedClassicalRS && classicalAlgo != defaultAlgo, Exceptions::RuntimeError, "\"aggregation: classical algo\" != default is not supported for scalled classical aggregation");
358 TEUCHOS_TEST_FOR_EXCEPTION(useSignedClassicalSA && classicalAlgo != defaultAlgo, Exceptions::RuntimeError, "\"aggregation: classical algo\" != default is not supported for scalled classical sa aggregation");
359
360 GO numDropped = 0, numTotal = 0;
361 std::string graphType = "unamalgamated"; // for description purposes only
362
363 /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
364 BlockSize is the number of storage blocks that must kept together during the amalgamation process.
365
366 Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
367
368 numPDEs = BlockSize * storageblocksize.
369
370 If numPDEs==1
371 Matrix is point storage (classical CRS storage). storageblocksize=1 and BlockSize=1
372 No other values makes sense.
373
374 If numPDEs>1
375 If matrix uses point storage, then storageblocksize=1 and BlockSize=numPDEs.
376 If matrix uses block storage, with block size of n, then storageblocksize=n, and BlockSize=numPDEs/n.
377 Thus far, only storageblocksize=numPDEs and BlockSize=1 has been tested.
378 */
379 TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0, Exceptions::RuntimeError, "A->GetFixedBlockSize() needs to be a multiple of A->GetStorageBlockSize()");
380 const LO BlockSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
381
382 /************************** RS or SA-style Classical Dropping (and variants) **************************/
383 if (algo == "classical") {
384 if (predrop_ == null) {
385 // ap: this is a hack: had to declare predrop_ as mutable
386 predrop_ = rcp(new PreDropFunctionConstVal(threshold));
387 }
388
389 if (predrop_ != null) {
390 RCP<PreDropFunctionConstVal> predropConstVal = rcp_dynamic_cast<PreDropFunctionConstVal>(predrop_);
391 TEUCHOS_TEST_FOR_EXCEPTION(predropConstVal == Teuchos::null, Exceptions::BadCast,
392 "MueLu::CoalesceFactory::Build: cast to PreDropFunctionConstVal failed.");
393 // If a user provided a predrop function, it overwrites the XML threshold parameter
394 SC newt = predropConstVal->GetThreshold();
395 if (newt != threshold) {
396 GetOStream(Warnings0) << "switching threshold parameter from " << threshold << " (list) to " << newt << " (user function" << std::endl;
397 threshold = newt;
398 }
399 }
400 // At this points we either have
401 // (predrop_ != null)
402 // Therefore, it is sufficient to check only threshold
403 if (BlockSize == 1 && threshold == STS::zero() && !useSignedClassicalRS && !useSignedClassicalSA && A->hasCrsGraph()) {
404 // Case 1: scalar problem, no dropping => just use matrix graph
405 RCP<LWGraph> graph = rcp(new LWGraph(A->getCrsGraph(), "graph of A"));
406 // Detect and record rows that correspond to Dirichlet boundary conditions
407 auto boundaryNodes = MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows_kokkos_host(*A, dirichletThreshold);
408 if (rowSumTol > 0.)
409 Utilities::ApplyRowSumCriterionHost(*A, rowSumTol, boundaryNodes);
410
411 graph->SetBoundaryNodeMap(boundaryNodes);
412 numTotal = A->getLocalNumEntries();
413
414 if (GetVerbLevel() & Statistics1) {
415 GO numLocalBoundaryNodes = 0;
416 GO numGlobalBoundaryNodes = 0;
417 for (size_t i = 0; i < boundaryNodes.size(); ++i)
418 if (boundaryNodes[i])
419 numLocalBoundaryNodes++;
420 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
421 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
422 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
423 }
424
425 Set(currentLevel, "DofsPerNode", 1);
426 Set(currentLevel, "Graph", graph);
427
428 } else if ((BlockSize == 1 && threshold != STS::zero()) ||
429 (BlockSize == 1 && threshold == STS::zero() && !A->hasCrsGraph()) ||
430 (BlockSize == 1 && useSignedClassicalRS) ||
431 (BlockSize == 1 && useSignedClassicalSA)) {
432 // Case 2: scalar problem with dropping => record the column indices of undropped entries, but still use original
433 // graph's map information, e.g., whether index is local
434 // OR a matrix without a CrsGraph
435
436 // allocate space for the local graph
437 typename LWGraph::row_type::non_const_type rows("rows", A->getLocalNumRows() + 1);
438 typename LWGraph::entries_type::non_const_type columns("columns", A->getLocalNumEntries());
439
440 using MT = typename STS::magnitudeType;
441 RCP<Vector> ghostedDiag;
442 ArrayRCP<const SC> ghostedDiagVals;
443 ArrayRCP<const SC> negMaxOffDiagonal;
444 // RS style needs the max negative off-diagonal, SA style needs the diagonal
445 if (useSignedClassicalRS) {
446 if (ghostedBlockNumber.is_null()) {
447 auto negMaxOffDiagonalVec = MueLu::Utilities<SC, LO, GO, NO>::GetMatrixMaxMinusOffDiagonal(*A);
448 negMaxOffDiagonal = negMaxOffDiagonalVec->getData(0);
449 if (GetVerbLevel() & Statistics1)
450 GetOStream(Statistics1) << "Calculated max point off-diagonal" << std::endl;
451 } else {
452 auto negMaxOffDiagonalVec = MueLu::Utilities<SC, LO, GO, NO>::GetMatrixMaxMinusOffDiagonal(*A, *ghostedBlockNumber);
453 negMaxOffDiagonal = negMaxOffDiagonalVec->getData(0);
454 if (GetVerbLevel() & Statistics1)
455 GetOStream(Statistics1) << "Calculating max block off-diagonal" << std::endl;
456 }
457 } else {
458 ghostedDiag = MueLu::Utilities<SC, LO, GO, NO>::GetMatrixOverlappedDiagonal(*A);
459 if (classicalAlgo == defaultAlgo) {
460 ghostedDiagVals = ghostedDiag->getData(0);
461 }
462 }
463 auto boundaryNodes = MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows_kokkos_host(*A, dirichletThreshold);
464 if (rowSumTol > 0.) {
465 if (ghostedBlockNumber.is_null()) {
466 if (GetVerbLevel() & Statistics1)
467 GetOStream(Statistics1) << "Applying point row sum criterion." << std::endl;
468 Utilities::ApplyRowSumCriterionHost(*A, rowSumTol, boundaryNodes);
469 } else {
470 if (GetVerbLevel() & Statistics1)
471 GetOStream(Statistics1) << "Applying block row sum criterion." << std::endl;
472 Utilities::ApplyRowSumCriterionHost(*A, *ghostedBlockNumber, rowSumTol, boundaryNodes);
473 }
474 }
475
476 ArrayRCP<const LO> g_block_id;
477 if (!ghostedBlockNumber.is_null())
478 g_block_id = ghostedBlockNumber->getData(0);
479
480 LO realnnz = 0;
481 rows(0) = 0;
482 if (classicalAlgo == defaultAlgo) {
483 SubFactoryMonitor m1(*this, "Classical RS/SA", currentLevel);
484 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
485 size_t nnz = A->getNumEntriesInLocalRow(row);
486 bool rowIsDirichlet = boundaryNodes[row];
487 ArrayView<const LO> indices;
488 ArrayView<const SC> vals;
489 A->getLocalRowView(row, indices, vals);
490
491 // FIXME the current predrop function uses the following
492 // FIXME if(std::abs(vals[k]) > std::abs(threshold_) || grow == gcid )
493 // FIXME but the threshold doesn't take into account the rows' diagonal entries
494 // FIXME For now, hardwiring the dropping in here
495
496 LO rownnz = 0;
497 if (useSignedClassicalRS) {
498 // Signed classical RS style
499 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
500 LO col = indices[colID];
501 MT max_neg_aik = realThreshold * STS::real(negMaxOffDiagonal[row]);
502 MT neg_aij = -STS::real(vals[colID]);
503 /* if(row==1326) printf("A(%d,%d) = %6.4e, block = (%d,%d) neg_aij = %6.4e max_neg_aik = %6.4e\n",row,col,vals[colID],
504 g_block_id.is_null() ? -1 : g_block_id[row],
505 g_block_id.is_null() ? -1 : g_block_id[col],
506 neg_aij, max_neg_aik);*/
507 if ((!rowIsDirichlet && (g_block_id.is_null() || g_block_id[row] == g_block_id[col]) && neg_aij > max_neg_aik) || row == col) {
508 columns[realnnz++] = col;
509 rownnz++;
510 } else
511 numDropped++;
512 }
513 rows(row + 1) = realnnz;
514 } else if (useSignedClassicalSA) {
515 // Signed classical SA style
516 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
517 LO col = indices[colID];
518
519 bool is_nonpositive = STS::real(vals[colID]) <= 0;
520 MT aiiajj = STS::magnitude(threshold * threshold * ghostedDiagVals[col] * ghostedDiagVals[row]); // eps^2*|a_ii|*|a_jj|
521 MT aij = is_nonpositive ? STS::magnitude(vals[colID] * vals[colID]) : (-STS::magnitude(vals[colID] * vals[colID])); // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
522 /*
523 if(row==1326) printf("A(%d,%d) = %6.4e, raw_aij = %6.4e aij = %6.4e aiiajj = %6.4e\n",row,col,vals[colID],
524 vals[colID],aij, aiiajj);
525 */
526
527 if ((!rowIsDirichlet && aij > aiiajj) || row == col) {
528 columns(realnnz++) = col;
529 rownnz++;
530 } else
531 numDropped++;
532 }
533 rows[row + 1] = realnnz;
534 } else {
535 // Standard abs classical
536 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
537 LO col = indices[colID];
538 MT aiiajj = STS::magnitude(threshold * threshold * ghostedDiagVals[col] * ghostedDiagVals[row]); // eps^2*|a_ii|*|a_jj|
539 MT aij = STS::magnitude(vals[colID] * vals[colID]); // |a_ij|^2
540
541 if ((!rowIsDirichlet && aij > aiiajj) || row == col) {
542 columns(realnnz++) = col;
543 rownnz++;
544 } else
545 numDropped++;
546 }
547 rows(row + 1) = realnnz;
548 }
549 } // end for row
550 } else {
551 /* Cut Algorithm */
552 SubFactoryMonitor m1(*this, "Cut Drop", currentLevel);
553 using ExecSpace = typename Node::execution_space;
554 using TeamPol = Kokkos::TeamPolicy<ExecSpace>;
555 using TeamMem = typename TeamPol::member_type;
556 using ATS = KokkosKernels::ArithTraits<Scalar>;
557 using impl_scalar_type = typename ATS::val_type;
558 using implATS = KokkosKernels::ArithTraits<impl_scalar_type>;
559
560 // move from host to device
561 auto ghostedDiagValsView = Kokkos::subview(ghostedDiag->getLocalViewDevice(Tpetra::Access::ReadOnly), Kokkos::ALL(), 0);
562 auto thresholdKokkos = static_cast<impl_scalar_type>(threshold);
563 auto realThresholdKokkos = implATS::magnitude(thresholdKokkos);
564 auto columnsDevice = Kokkos::create_mirror_view(ExecSpace(), columns);
565
566 auto A_device = A->getLocalMatrixDevice();
567 RCP<LWGraph> graph = rcp(new LWGraph(A->getCrsGraph(), "graph of A"));
568 RCP<const Import> importer = A->getCrsGraph()->getImporter();
569 RCP<LocalOrdinalVector> boundaryNodesVector = Xpetra::VectorFactory<LO, LO, GO, NO>::Build(graph->GetDomainMap());
570 RCP<LocalOrdinalVector> boundaryColumnVector;
571 for (size_t i = 0; i < graph->GetNodeNumVertices(); i++) {
572 boundaryNodesVector->getDataNonConst(0)[i] = boundaryNodes[i];
573 }
574 if (!importer.is_null()) {
575 boundaryColumnVector = Xpetra::VectorFactory<LO, LO, GO, NO>::Build(graph->GetImportMap());
576 boundaryColumnVector->doImport(*boundaryNodesVector, *importer, Xpetra::INSERT);
577 } else {
578 boundaryColumnVector = boundaryNodesVector;
579 }
580 auto boundaryColumn = boundaryColumnVector->getLocalViewDevice(Tpetra::Access::ReadOnly);
581 auto boundary = Kokkos::subview(boundaryColumn, Kokkos::ALL(), 0);
582
583 Kokkos::View<LO*, ExecSpace> rownnzView("rownnzView", A_device.numRows());
584 auto drop_views = Kokkos::View<bool*, ExecSpace>("drop_views", A_device.nnz());
585 auto index_views = Kokkos::View<size_t*, ExecSpace>("index_views", A_device.nnz());
586
587 Kokkos::parallel_reduce(
588 "classical_cut", TeamPol(A_device.numRows(), Kokkos::AUTO), KOKKOS_LAMBDA(const TeamMem& teamMember, LO& globalnnz, GO& totalDropped) {
589 LO row = teamMember.league_rank();
590 auto rowView = A_device.rowConst(row);
591 size_t nnz = rowView.length;
592
593 auto drop_view = Kokkos::subview(drop_views, Kokkos::make_pair(A_device.graph.row_map(row), A_device.graph.row_map(row + 1)));
594 auto index_view = Kokkos::subview(index_views, Kokkos::make_pair(A_device.graph.row_map(row), A_device.graph.row_map(row + 1)));
595
596 // find magnitudes
597 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, (LO)nnz), [&](const LO colID) {
598 index_view(colID) = colID;
599 LO col = rowView.colidx(colID);
600 // ignore diagonals for now, they are checked again later
601 // Don't aggregate boundaries
602 if (row == col || boundary(col)) {
603 drop_view(colID) = true;
604 } else {
605 drop_view(colID) = false;
606 }
607 });
608
609 size_t dropStart = nnz;
610 if (classicalAlgo == unscaled_cut) {
611 // push diagonals and boundaries to the right, sort everything else by aij on the left
612 Kokkos::Experimental::sort_team(teamMember, index_view, [=](size_t& x, size_t& y) -> bool {
613 if (drop_view(x) || drop_view(y)) {
614 return drop_view(x) < drop_view(y);
615 } else {
616 auto x_aij = implATS::magnitude(rowView.value(x) * rowView.value(x));
617 auto y_aij = implATS::magnitude(rowView.value(y) * rowView.value(y));
618 return x_aij > y_aij;
619 }
620 });
621
622 // find index where dropping starts
623 Kokkos::parallel_reduce(
624 Kokkos::TeamThreadRange(teamMember, 1, nnz), [=](size_t i, size_t& min) {
625 auto const& x = index_view(i - 1);
626 auto const& y = index_view(i);
627 typename implATS::magnitudeType x_aij = 0;
628 typename implATS::magnitudeType y_aij = 0;
629 if (!drop_view(x)) {
630 x_aij = implATS::magnitude(rowView.value(x) * rowView.value(x));
631 }
632 if (!drop_view(y)) {
633 y_aij = implATS::magnitude(rowView.value(y) * rowView.value(y));
634 }
635
636 if (realThresholdKokkos * realThresholdKokkos * x_aij > y_aij) {
637 if (i < min) {
638 min = i;
639 }
640 }
641 },
642 Kokkos::Min<size_t>(dropStart));
643 } else if (classicalAlgo == scaled_cut) {
644 // push diagonals and boundaries to the right, sort everything else by aij/aiiajj on the left
645 Kokkos::Experimental::sort_team(teamMember, index_view, [=](size_t& x, size_t& y) -> bool {
646 if (drop_view(x) || drop_view(y)) {
647 return drop_view(x) < drop_view(y);
648 } else {
649 auto x_aij = implATS::magnitude(rowView.value(x) * rowView.value(x));
650 auto y_aij = implATS::magnitude(rowView.value(y) * rowView.value(y));
651 auto x_aiiajj = implATS::magnitude(ghostedDiagValsView(rowView.colidx(x)) * ghostedDiagValsView(row));
652 auto y_aiiajj = implATS::magnitude(ghostedDiagValsView(rowView.colidx(y)) * ghostedDiagValsView(row));
653 return (x_aij / x_aiiajj) > (y_aij / y_aiiajj);
654 }
655 });
656
657 // find index where dropping starts
658 Kokkos::parallel_reduce(
659 Kokkos::TeamThreadRange(teamMember, 1, nnz), [=](size_t i, size_t& min) {
660 auto const& x = index_view(i - 1);
661 auto const& y = index_view(i);
662 typename implATS::magnitudeType x_val = 0;
663 typename implATS::magnitudeType y_val = 0;
664 if (!drop_view(x)) {
665 typename implATS::magnitudeType x_aij = implATS::magnitude(rowView.value(x) * rowView.value(x));
666 typename implATS::magnitudeType x_aiiajj = implATS::magnitude(ghostedDiagValsView(rowView.colidx(x)) * ghostedDiagValsView(row));
667 x_val = x_aij / x_aiiajj;
668 }
669 if (!drop_view(y)) {
670 typename implATS::magnitudeType y_aij = implATS::magnitude(rowView.value(y) * rowView.value(y));
671 typename implATS::magnitudeType y_aiiajj = implATS::magnitude(ghostedDiagValsView(rowView.colidx(y)) * ghostedDiagValsView(row));
672 y_val = y_aij / y_aiiajj;
673 }
674
675 if (realThresholdKokkos * realThresholdKokkos * x_val > y_val) {
676 if (i < min) {
677 min = i;
678 }
679 }
680 },
681 Kokkos::Min<size_t>(dropStart));
682 }
683
684 // drop everything to the right of where values stop passing threshold
685 if (dropStart < nnz) {
686 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, dropStart, nnz), [=](size_t i) {
687 drop_view(index_view(i)) = true;
688 });
689 }
690
691 LO rownnz = 0;
692 GO rowDropped = 0;
693 Kokkos::parallel_reduce(
694 Kokkos::TeamThreadRange(teamMember, nnz), [=](const size_t idxID, LO& keep, GO& drop) {
695 LO col = rowView.colidx(idxID);
696 // don't drop diagonal
697 if (row == col || !drop_view(idxID)) {
698 columnsDevice(A_device.graph.row_map(row) + idxID) = col;
699 keep++;
700 } else {
701 columnsDevice(A_device.graph.row_map(row) + idxID) = -1;
702 drop++;
703 }
704 },
705 rownnz, rowDropped);
706
707 Kokkos::single(Kokkos::PerTeam(teamMember), [&]() {
708 globalnnz += rownnz;
709 totalDropped += rowDropped;
710 rownnzView(row) = rownnz;
711 });
712 },
713 realnnz, numDropped);
714
715 // update column indices so that kept indices are aligned to the left for subview that happens later on
716 Kokkos::Experimental::remove(ExecSpace(), columnsDevice, -1);
717 Kokkos::deep_copy(columns, columnsDevice);
718
719 // update row indices by adding up new # of nnz in each row
720 auto rowsDevice = Kokkos::create_mirror_view(ExecSpace(), rows);
721 Kokkos::parallel_scan(
722 Kokkos::RangePolicy<ExecSpace>(0, A_device.numRows()), KOKKOS_LAMBDA(const int i, LO& partial_sum, bool is_final) {
723 partial_sum += rownnzView(i);
724 if (is_final) rowsDevice(i + 1) = partial_sum;
725 });
726 Kokkos::deep_copy(rows, rowsDevice);
727 }
728
729 numTotal = A->getLocalNumEntries();
730
731 if (aggregationMayCreateDirichlet) {
732 // If the only element remaining after filtering is diagonal, mark node as boundary
733 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
734 if (rows[row + 1] - rows[row] <= 1)
735 boundaryNodes[row] = true;
736 }
737 }
738
739 RCP<LWGraph> graph = rcp(new LWGraph(rows, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), A->getRowMap(), A->getColMap(), "thresholded graph of A"));
740 graph->SetBoundaryNodeMap(boundaryNodes);
741 if (GetVerbLevel() & Statistics1) {
742 GO numLocalBoundaryNodes = 0;
743 GO numGlobalBoundaryNodes = 0;
744 for (size_t i = 0; i < boundaryNodes.size(); ++i)
745 if (boundaryNodes(i))
746 numLocalBoundaryNodes++;
747 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
748 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
749 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
750 }
751 Set(currentLevel, "Graph", graph);
752 Set(currentLevel, "DofsPerNode", 1);
753
754 // If we're doing signed classical, we might want to block-diagonalize *after* the dropping
755 if (generateColoringGraph) {
756 RCP<LWGraph> colorGraph;
757 RCP<const Import> importer = A->getCrsGraph()->getImporter();
758 BlockDiagonalizeGraph(graph, ghostedBlockNumber, colorGraph, importer);
759 Set(currentLevel, "Coloring Graph", colorGraph);
760 // #define CMS_DUMP
761#ifdef CMS_DUMP
762 {
763 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("m_regular_graph." + std::to_string(currentLevel.GetLevelID()), *rcp_dynamic_cast<LWGraph>(graph)->GetCrsGraph());
764 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("m_color_graph." + std::to_string(currentLevel.GetLevelID()), *rcp_dynamic_cast<LWGraph>(colorGraph)->GetCrsGraph());
765 // int rank = graph->GetDomainMap()->getComm()->getRank();
766 // {
767 // std::ofstream ofs(std::string("m_color_graph_") + std::to_string(currentLevel.GetLevelID())+std::string("_") + std::to_string(rank) + std::string(".dat"),std::ofstream::out);
768 // RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(ofs));
769 // colorGraph->print(*fancy,Debug);
770 // }
771 // {
772 // std::ofstream ofs(std::string("m_regular_graph_") + std::to_string(currentLevel.GetLevelID())+std::string("_") + std::to_string(rank) + std::string(".dat"),std::ofstream::out);
773 // RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(ofs));
774 // graph->print(*fancy,Debug);
775 // }
776 }
777#endif
778 } // end generateColoringGraph
779 } else if (BlockSize > 1 && threshold == STS::zero()) {
780 // Case 3: Multiple DOF/node problem without dropping
781 const RCP<const Map> rowMap = A->getRowMap();
782 const RCP<const Map> colMap = A->getColMap();
783
784 graphType = "amalgamated";
785
786 // build node row map (uniqueMap) and node column map (nonUniqueMap)
787 // the arrays rowTranslation and colTranslation contain the local node id
788 // given a local dof id. The data is calculated by the AmalgamationFactory and
789 // stored in the variable container "UnAmalgamationInfo"
790 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
791 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
792 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
793 Array<LO> colTranslation = *(amalInfo->getColTranslation());
794
795 // get number of local nodes
796 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
797
798 // Allocate space for the local graph
799 typename LWGraph::row_type::non_const_type rows("rows", numRows + 1);
800 typename LWGraph::entries_type::non_const_type columns("columns", A->getLocalNumEntries());
801
802 typename LWGraph::boundary_nodes_type amalgBoundaryNodes("amalgBoundaryNodes", numRows);
803 Kokkos::deep_copy(amalgBoundaryNodes, false);
804
805 // Detect and record rows that correspond to Dirichlet boundary conditions
806 // TODO If we use ArrayRCP<LO>, then we can record boundary nodes as usual. Size
807 // TODO the array one bigger than the number of local rows, and the last entry can
808 // TODO hold the actual number of boundary nodes. Clever, huh?
809 ArrayRCP<bool> pointBoundaryNodes;
810 pointBoundaryNodes = Teuchos::arcp_const_cast<bool>(MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows(*A, dirichletThreshold));
811 if (rowSumTol > 0.)
812 Utilities::ApplyRowSumCriterion(*A, rowSumTol, pointBoundaryNodes);
813
814 // extract striding information
815 LO blkSize = A->GetFixedBlockSize(); //< the full block size (number of dofs per node in strided map)
816 LO blkId = -1; //< the block id within the strided map (or -1 if it is a full block map)
817 LO blkPartSize = A->GetFixedBlockSize(); //< stores the size of the block within the strided map
818 if (A->IsView("stridedMaps") == true) {
819 Teuchos::RCP<const Map> myMap = A->getRowMap("stridedMaps");
820 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
821 TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
822 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
823 blkId = strMap->getStridedBlockId();
824 if (blkId > -1)
825 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
826 }
827
828 // loop over all local nodes
829 LO realnnz = 0;
830 rows(0) = 0;
831 Array<LO> indicesExtra;
832 for (LO row = 0; row < numRows; row++) {
833 ArrayView<const LO> indices;
834 indicesExtra.resize(0);
835
836 // The amalgamated row is marked as Dirichlet iff all point rows are Dirichlet
837 // Note, that pointBoundaryNodes lives on the dofmap (and not the node map).
838 // Therefore, looping over all dofs is fine here. We use blkPartSize as we work
839 // with local ids.
840 // TODO: Here we have different options of how to define a node to be a boundary (or Dirichlet)
841 // node.
842 bool isBoundary = false;
843 if (pL.get<bool>("aggregation: greedy Dirichlet") == true) {
844 for (LO j = 0; j < blkPartSize; j++) {
845 if (pointBoundaryNodes[row * blkPartSize + j]) {
846 isBoundary = true;
847 break;
848 }
849 }
850 } else {
851 isBoundary = true;
852 for (LO j = 0; j < blkPartSize; j++) {
853 if (!pointBoundaryNodes[row * blkPartSize + j]) {
854 isBoundary = false;
855 break;
856 }
857 }
858 }
859
860 // Merge rows of A
861 // The array indicesExtra contains local column node ids for the current local node "row"
862 if (!isBoundary)
863 MergeRows(*A, row, indicesExtra, colTranslation);
864 else
865 indicesExtra.push_back(row);
866 indices = indicesExtra;
867 numTotal += indices.size();
868
869 // add the local column node ids to the full columns array which
870 // contains the local column node ids for all local node rows
871 LO nnz = indices.size(), rownnz = 0;
872 for (LO colID = 0; colID < nnz; colID++) {
873 LO col = indices[colID];
874 columns(realnnz++) = col;
875 rownnz++;
876 }
877
878 if (rownnz == 1) {
879 // If the only element remaining after filtering is diagonal, mark node as boundary
880 // FIXME: this should really be replaced by the following
881 // if (indices.size() == 1 && indices[0] == row)
882 // boundaryNodes[row] = true;
883 // We do not do it this way now because there is no framework for distinguishing isolated
884 // and boundary nodes in the aggregation algorithms
885 amalgBoundaryNodes[row] = true;
886 }
887 rows(row + 1) = realnnz;
888 } // for (LO row = 0; row < numRows; row++)
889
890 RCP<LWGraph> graph = rcp(new LWGraph(rows, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), uniqueMap, nonUniqueMap, "amalgamated graph of A"));
891 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
892
893 if (GetVerbLevel() & Statistics1) {
894 GO numLocalBoundaryNodes = 0;
895 GO numGlobalBoundaryNodes = 0;
896
897 for (size_t i = 0; i < amalgBoundaryNodes.size(); ++i)
898 if (amalgBoundaryNodes(i))
899 numLocalBoundaryNodes++;
900
901 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
902 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
903 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes
904 << " agglomerated Dirichlet nodes" << std::endl;
905 }
906
907 Set(currentLevel, "Graph", graph);
908 Set(currentLevel, "DofsPerNode", blkSize); // full block size
909
910 } else if (BlockSize > 1 && threshold != STS::zero()) {
911 // Case 4: Multiple DOF/node problem with dropping
912 const RCP<const Map> rowMap = A->getRowMap();
913 const RCP<const Map> colMap = A->getColMap();
914 graphType = "amalgamated";
915
916 // build node row map (uniqueMap) and node column map (nonUniqueMap)
917 // the arrays rowTranslation and colTranslation contain the local node id
918 // given a local dof id. The data is calculated by the AmalgamationFactory and
919 // stored in the variable container "UnAmalgamationInfo"
920 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
921 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
922 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
923 Array<LO> colTranslation = *(amalInfo->getColTranslation());
924
925 // get number of local nodes
926 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
927
928 // Allocate space for the local graph
929 typename LWGraph::row_type::non_const_type rows("rows", numRows + 1);
930 typename LWGraph::entries_type::non_const_type columns("columns", A->getLocalNumEntries());
931
932 typename LWGraph::boundary_nodes_type amalgBoundaryNodes("amalgBoundaryNodes", numRows);
933 Kokkos::deep_copy(amalgBoundaryNodes, false);
934
935 // Detect and record rows that correspond to Dirichlet boundary conditions
936 // TODO If we use ArrayRCP<LO>, then we can record boundary nodes as usual. Size
937 // TODO the array one bigger than the number of local rows, and the last entry can
938 // TODO hold the actual number of boundary nodes. Clever, huh?
939 auto pointBoundaryNodes = MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows_kokkos_host(*A, dirichletThreshold);
940 if (rowSumTol > 0.)
941 Utilities::ApplyRowSumCriterionHost(*A, rowSumTol, pointBoundaryNodes);
942
943 // extract striding information
944 LO blkSize = A->GetFixedBlockSize(); //< the full block size (number of dofs per node in strided map)
945 LO blkId = -1; //< the block id within the strided map (or -1 if it is a full block map)
946 LO blkPartSize = A->GetFixedBlockSize(); //< stores the size of the block within the strided map
947 if (A->IsView("stridedMaps") == true) {
948 Teuchos::RCP<const Map> myMap = A->getRowMap("stridedMaps");
949 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
950 TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
951 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
952 blkId = strMap->getStridedBlockId();
953 if (blkId > -1)
954 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
955 }
956
957 // extract diagonal data for dropping strategy
958 RCP<Vector> ghostedDiag = MueLu::Utilities<SC, LO, GO, NO>::GetMatrixOverlappedDiagonal(*A);
959 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
960
961 // loop over all local nodes
962 LO realnnz = 0;
963 rows[0] = 0;
964 Array<LO> indicesExtra;
965 for (LO row = 0; row < numRows; row++) {
966 ArrayView<const LO> indices;
967 indicesExtra.resize(0);
968
969 // The amalgamated row is marked as Dirichlet iff all point rows are Dirichlet
970 // Note, that pointBoundaryNodes lives on the dofmap (and not the node map).
971 // Therefore, looping over all dofs is fine here. We use blkPartSize as we work
972 // with local ids.
973 // TODO: Here we have different options of how to define a node to be a boundary (or Dirichlet)
974 // node.
975 bool isBoundary = false;
976 if (pL.get<bool>("aggregation: greedy Dirichlet") == true) {
977 for (LO j = 0; j < blkPartSize; j++) {
978 if (pointBoundaryNodes[row * blkPartSize + j]) {
979 isBoundary = true;
980 break;
981 }
982 }
983 } else {
984 isBoundary = true;
985 for (LO j = 0; j < blkPartSize; j++) {
986 if (!pointBoundaryNodes[row * blkPartSize + j]) {
987 isBoundary = false;
988 break;
989 }
990 }
991 }
992
993 // Merge rows of A
994 // The array indicesExtra contains local column node ids for the current local node "row"
995 if (!isBoundary)
996 MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
997 else
998 indicesExtra.push_back(row);
999 indices = indicesExtra;
1000 numTotal += indices.size();
1001
1002 // add the local column node ids to the full columns array which
1003 // contains the local column node ids for all local node rows
1004 LO nnz = indices.size(), rownnz = 0;
1005 for (LO colID = 0; colID < nnz; colID++) {
1006 LO col = indices[colID];
1007 columns[realnnz++] = col;
1008 rownnz++;
1009 }
1010
1011 if (rownnz == 1) {
1012 // If the only element remaining after filtering is diagonal, mark node as boundary
1013 // FIXME: this should really be replaced by the following
1014 // if (indices.size() == 1 && indices[0] == row)
1015 // boundaryNodes[row] = true;
1016 // We do not do it this way now because there is no framework for distinguishing isolated
1017 // and boundary nodes in the aggregation algorithms
1018 amalgBoundaryNodes[row] = true;
1019 }
1020 rows[row + 1] = realnnz;
1021 } // for (LO row = 0; row < numRows; row++)
1022 // columns.resize(realnnz);
1023
1024 RCP<LWGraph> graph = rcp(new LWGraph(rows, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), uniqueMap, nonUniqueMap, "amalgamated graph of A"));
1025 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1026
1027 if (GetVerbLevel() & Statistics1) {
1028 GO numLocalBoundaryNodes = 0;
1029 GO numGlobalBoundaryNodes = 0;
1030
1031 for (size_t i = 0; i < amalgBoundaryNodes.size(); ++i)
1032 if (amalgBoundaryNodes(i))
1033 numLocalBoundaryNodes++;
1034
1035 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
1036 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1037 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes
1038 << " agglomerated Dirichlet nodes" << std::endl;
1039 }
1040
1041 Set(currentLevel, "Graph", graph);
1042 Set(currentLevel, "DofsPerNode", blkSize); // full block size
1043 }
1044
1045 } else if (algo == "distance laplacian") {
1046 LO blkSize = A->GetFixedBlockSize();
1047 GO indexBase = A->getRowMap()->getIndexBase();
1048 // [*0*] : FIXME
1049 // ap: somehow, if I move this line to [*1*], Belos throws an error
1050 // I'm not sure what's going on. Do we always have to Get data, if we did
1051 // DeclareInput for it?
1052 // RCP<RealValuedMultiVector> Coords = Get< RCP<RealValuedMultiVector > >(currentLevel, "Coordinates");
1053
1054 // Detect and record rows that correspond to Dirichlet boundary conditions
1055 // TODO If we use ArrayRCP<LO>, then we can record boundary nodes as usual. Size
1056 // TODO the array one bigger than the number of local rows, and the last entry can
1057 // TODO hold the actual number of boundary nodes. Clever, huh?
1058 auto pointBoundaryNodes = MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows_kokkos_host(*A, dirichletThreshold);
1059 if (rowSumTol > 0.)
1060 Utilities::ApplyRowSumCriterionHost(*A, rowSumTol, pointBoundaryNodes);
1061
1062 if ((blkSize == 1) && (threshold == STS::zero())) {
1063 // Trivial case: scalar problem, no dropping. Can return original graph
1064 RCP<LWGraph> graph = rcp(new LWGraph(A->getCrsGraph(), "graph of A"));
1065 graph->SetBoundaryNodeMap(pointBoundaryNodes);
1066 graphType = "unamalgamated";
1067 numTotal = A->getLocalNumEntries();
1068
1069 if (GetVerbLevel() & Statistics1) {
1070 GO numLocalBoundaryNodes = 0;
1071 GO numGlobalBoundaryNodes = 0;
1072 for (size_t i = 0; i < pointBoundaryNodes.size(); ++i)
1073 if (pointBoundaryNodes(i))
1074 numLocalBoundaryNodes++;
1075 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
1076 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1077 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
1078 }
1079
1080 Set(currentLevel, "DofsPerNode", blkSize);
1081 Set(currentLevel, "Graph", graph);
1082
1083 } else {
1084 // ap: We make quite a few assumptions here; general case may be a lot different,
1085 // but much much harder to implement. We assume that:
1086 // 1) all maps are standard maps, not strided maps
1087 // 2) global indices of dofs in A are related to dofs in coordinates in a simple arithmetic
1088 // way: rows i*blkSize, i*blkSize+1, ..., i*blkSize + (blkSize-1) correspond to node i
1089 //
1090 // NOTE: Potentially, some of the code below could be simplified with UnAmalgamationInfo,
1091 // but as I totally don't understand that code, here is my solution
1092
1093 // [*1*]: see [*0*]
1094
1095 // Check that the number of local coordinates is consistent with the #rows in A
1096 TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getLocalNumElements() / blkSize != Coords->getLocalLength(), Exceptions::Incompatible,
1097 "Coordinate vector length (" << Coords->getLocalLength() << ") is incompatible with number of rows in A (" << A->getRowMap()->getLocalNumElements() << ") by modulo block size (" << blkSize << ").");
1098
1099 const RCP<const Map> colMap = A->getColMap();
1100 RCP<const Map> uniqueMap, nonUniqueMap;
1101 Array<LO> colTranslation;
1102 if (blkSize == 1) {
1103 uniqueMap = A->getRowMap();
1104 nonUniqueMap = A->getColMap();
1105 graphType = "unamalgamated";
1106
1107 } else {
1108 uniqueMap = Coords->getMap();
1109 TEUCHOS_TEST_FOR_EXCEPTION(uniqueMap->getIndexBase() != indexBase, Exceptions::Incompatible,
1110 "Different index bases for matrix and coordinates");
1111
1112 AmalgamationFactory::AmalgamateMap(*(A->getColMap()), *A, nonUniqueMap, colTranslation);
1113
1114 graphType = "amalgamated";
1115 }
1116 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
1117
1118 RCP<RealValuedMultiVector> ghostedCoords;
1119 RCP<Vector> ghostedLaplDiag;
1120 Teuchos::ArrayRCP<SC> ghostedLaplDiagData;
1121 if (threshold != STS::zero()) {
1122 // Get ghost coordinates
1123 RCP<const Import> importer;
1124 {
1125 SubFactoryMonitor m1(*this, "Import construction", currentLevel);
1126 if (blkSize == 1 && realA->getCrsGraph()->getImporter() != Teuchos::null) {
1127 GetOStream(Warnings1) << "Using existing importer from matrix graph" << std::endl;
1128 importer = realA->getCrsGraph()->getImporter();
1129 } else {
1130 GetOStream(Warnings0) << "Constructing new importer instance" << std::endl;
1131 importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
1132 }
1133 } // subtimer
1134 ghostedCoords = Xpetra::MultiVectorFactory<real_type, LO, GO, NO>::Build(nonUniqueMap, Coords->getNumVectors());
1135 {
1136 SubFactoryMonitor m1(*this, "Coordinate import", currentLevel);
1137 ghostedCoords->doImport(*Coords, *importer, Xpetra::INSERT);
1138 } // subtimer
1139
1140 // Construct Distance Laplacian diagonal
1141 RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
1142 Array<LO> indicesExtra;
1143 Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
1144 if (threshold != STS::zero()) {
1145 const size_t numVectors = ghostedCoords->getNumVectors();
1146 coordData.reserve(numVectors);
1147 for (size_t j = 0; j < numVectors; j++) {
1148 Teuchos::ArrayRCP<const real_type> tmpData = ghostedCoords->getData(j);
1149 coordData.push_back(tmpData);
1150 }
1151 }
1152 {
1153 SubFactoryMonitor m1(*this, "Laplacian local diagonal", currentLevel);
1154 ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
1155 for (LO row = 0; row < numRows; row++) {
1156 ArrayView<const LO> indices;
1157
1158 if (blkSize == 1) {
1159 ArrayView<const SC> vals;
1160 A->getLocalRowView(row, indices, vals);
1161
1162 } else {
1163 // Merge rows of A
1164 indicesExtra.resize(0);
1165 MergeRows(*A, row, indicesExtra, colTranslation);
1166 indices = indicesExtra;
1167 }
1168
1169 LO nnz = indices.size();
1170 bool haveAddedToDiag = false;
1171 for (LO colID = 0; colID < nnz; colID++) {
1172 const LO col = indices[colID];
1173
1174 if (row != col) {
1175 if (use_dlap_weights == SINGLE_WEIGHTS) {
1176 /*printf("[%d,%d] Unweighted Distance = %6.4e Weighted Distance = %6.4e\n",row,col,
1177 MueLu::Utilities<real_type,LO,GO,NO>::Distance2(coordData, row, col),
1178 MueLu::Utilities<real_type,LO,GO,NO>::Distance2(dlap_weights(),coordData, row, col));*/
1179 localLaplDiagData[row] += STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(dlap_weights(), coordData, row, col);
1180 } else if (use_dlap_weights == BLOCK_WEIGHTS) {
1181 int block_id = row % interleaved_blocksize;
1182 int block_start = block_id * interleaved_blocksize;
1183 localLaplDiagData[row] += STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(dlap_weights(block_start, interleaved_blocksize), coordData, row, col);
1184 } else {
1185 // printf("[%d,%d] Unweighted Distance = %6.4e\n",row,col,MueLu::Utilities<real_type,LO,GO,NO>::Distance2(coordData, row, col));
1186 localLaplDiagData[row] += STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(coordData, row, col);
1187 }
1188 haveAddedToDiag = true;
1189 }
1190 }
1191 // Deal with the situation where boundary conditions have only been enforced on rows, but not on columns.
1192 // We enforce dropping of these entries by assigning a very large number to the diagonal entries corresponding to BCs.
1193 if (!haveAddedToDiag)
1194 localLaplDiagData[row] = STS::squareroot(STS::rmax());
1195 }
1196 } // subtimer
1197 {
1198 SubFactoryMonitor m1(*this, "Laplacian distributed diagonal", currentLevel);
1199 ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
1200 ghostedLaplDiag->doImport(*localLaplDiag, *importer, Xpetra::INSERT);
1201 ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
1202 } // subtimer
1203
1204 } else {
1205 GetOStream(Runtime0) << "Skipping distance laplacian construction due to 0 threshold" << std::endl;
1206 }
1207
1208 // NOTE: ghostedLaplDiagData might be zero if we don't actually calculate the laplacian
1209
1210 // allocate space for the local graph
1211 typename LWGraph::row_type::non_const_type rows("rows", numRows + 1);
1212 typename LWGraph::entries_type::non_const_type columns("columns", A->getLocalNumEntries());
1213
1214#ifdef HAVE_MUELU_DEBUG
1215 // DEBUGGING
1216 for (LO i = 0; i < (LO)columns.size(); i++) columns[i] = -666;
1217#endif
1218
1219 // Extra array for if we're allowing symmetrization with cutting
1220 ArrayRCP<LO> rows_stop;
1221 bool use_stop_array = threshold != STS::zero() && distanceLaplacianAlgo == scaled_cut_symmetric;
1222 if (use_stop_array)
1223 // rows_stop = typename LWGraph::row_type::non_const_type("rows_stop", numRows);
1224 rows_stop.resize(numRows);
1225
1226 typename LWGraph::boundary_nodes_type amalgBoundaryNodes("amalgBoundaryNodes", numRows);
1227 Kokkos::deep_copy(amalgBoundaryNodes, false);
1228
1229 LO realnnz = 0;
1230 rows(0) = 0;
1231
1232 Array<LO> indicesExtra;
1233 {
1234 SubFactoryMonitor m1(*this, "Laplacian dropping", currentLevel);
1235 Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
1236 if (threshold != STS::zero()) {
1237 const size_t numVectors = ghostedCoords->getNumVectors();
1238 coordData.reserve(numVectors);
1239 for (size_t j = 0; j < numVectors; j++) {
1240 Teuchos::ArrayRCP<const real_type> tmpData = ghostedCoords->getData(j);
1241 coordData.push_back(tmpData);
1242 }
1243 }
1244
1245 ArrayView<const SC> vals; // CMS hackery
1246 for (LO row = 0; row < numRows; row++) {
1247 ArrayView<const LO> indices;
1248 indicesExtra.resize(0);
1249 bool isBoundary = false;
1250
1251 if (blkSize == 1) {
1252 // ArrayView<const SC> vals;//CMS uncomment
1253 A->getLocalRowView(row, indices, vals);
1254 isBoundary = pointBoundaryNodes[row];
1255 } else {
1256 // The amalgamated row is marked as Dirichlet iff all point rows are Dirichlet
1257 isBoundary = true;
1258 for (LO j = 0; j < blkSize; j++) {
1259 if (!pointBoundaryNodes[row * blkSize + j]) {
1260 isBoundary = false;
1261 break;
1262 }
1263 }
1264
1265 // Merge rows of A
1266 if (!isBoundary)
1267 MergeRows(*A, row, indicesExtra, colTranslation);
1268 else
1269 indicesExtra.push_back(row);
1270 indices = indicesExtra;
1271 }
1272 numTotal += indices.size();
1273
1274 LO nnz = indices.size(), rownnz = 0;
1275
1276 if (use_stop_array) {
1277 rows(row + 1) = rows(row) + nnz;
1278 realnnz = rows(row);
1279 }
1280
1281 if (threshold != STS::zero()) {
1282 // default
1283 if (distanceLaplacianAlgo == defaultAlgo) {
1284 /* Standard Distance Laplacian */
1285 for (LO colID = 0; colID < nnz; colID++) {
1286 LO col = indices[colID];
1287
1288 if (row == col) {
1289 columns(realnnz++) = col;
1290 rownnz++;
1291 continue;
1292 }
1293
1294 // We do not want the distance Laplacian aggregating boundary nodes
1295 if (isBoundary) continue;
1296
1297 SC laplVal;
1298 if (use_dlap_weights == SINGLE_WEIGHTS) {
1299 laplVal = STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(dlap_weights(), coordData, row, col);
1300 } else if (use_dlap_weights == BLOCK_WEIGHTS) {
1301 int block_id = row % interleaved_blocksize;
1302 int block_start = block_id * interleaved_blocksize;
1303 laplVal = STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(dlap_weights(block_start, interleaved_blocksize), coordData, row, col);
1304 } else {
1305 laplVal = STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(coordData, row, col);
1306 }
1307 real_type aiiajj = STS::magnitude(realThreshold * realThreshold * ghostedLaplDiagData[row] * ghostedLaplDiagData[col]);
1308 real_type aij = STS::magnitude(laplVal * laplVal);
1309
1310 if (aij > aiiajj) {
1311 columns(realnnz++) = col;
1312 rownnz++;
1313 } else {
1314 numDropped++;
1315 }
1316 }
1317 } else {
1318 /* Cut Algorithm */
1319 using DropTol = Details::DropTol<real_type, LO>;
1320 std::vector<DropTol> drop_vec;
1321 drop_vec.reserve(nnz);
1322 const real_type zero = Teuchos::ScalarTraits<real_type>::zero();
1323 const real_type one = Teuchos::ScalarTraits<real_type>::one();
1324
1325 // find magnitudes
1326 for (LO colID = 0; colID < nnz; colID++) {
1327 LO col = indices[colID];
1328
1329 if (row == col) {
1330 drop_vec.emplace_back(zero, one, colID, false);
1331 continue;
1332 }
1333 // We do not want the distance Laplacian aggregating boundary nodes
1334 if (isBoundary) continue;
1335
1336 SC laplVal;
1337 if (use_dlap_weights == SINGLE_WEIGHTS) {
1338 laplVal = STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(dlap_weights(), coordData, row, col);
1339 } else if (use_dlap_weights == BLOCK_WEIGHTS) {
1340 int block_id = row % interleaved_blocksize;
1341 int block_start = block_id * interleaved_blocksize;
1342 laplVal = STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(dlap_weights(block_start, interleaved_blocksize), coordData, row, col);
1343 } else {
1344 laplVal = STS::one() / MueLu::Utilities<real_type, LO, GO, NO>::Distance2(coordData, row, col);
1345 }
1346
1347 real_type aiiajj = STS::magnitude(ghostedLaplDiagData[row] * ghostedLaplDiagData[col]);
1348 real_type aij = STS::magnitude(laplVal * laplVal);
1349
1350 drop_vec.emplace_back(aij, aiiajj, colID, false);
1351 }
1352
1353 const size_t n = drop_vec.size();
1354
1355 if (distanceLaplacianAlgo == unscaled_cut) {
1356 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol const& a, DropTol const& b) {
1357 return a.val > b.val;
1358 });
1359
1360 bool drop = false;
1361 for (size_t i = 1; i < n; ++i) {
1362 if (!drop) {
1363 auto const& x = drop_vec[i - 1];
1364 auto const& y = drop_vec[i];
1365 auto a = x.val;
1366 auto b = y.val;
1367 if (realThreshold * realThreshold * a > b) {
1368 drop = true;
1369#ifdef HAVE_MUELU_DEBUG
1370 if (distanceLaplacianCutVerbose) {
1371 std::cout << "DJS: KEEP, N, ROW: " << i + 1 << ", " << n << ", " << row << std::endl;
1372 }
1373#endif
1374 }
1375 }
1376 drop_vec[i].drop = drop;
1377 }
1378 } else if (distanceLaplacianAlgo == scaled_cut || distanceLaplacianAlgo == scaled_cut_symmetric) {
1379 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol const& a, DropTol const& b) {
1380 return a.val / a.diag > b.val / b.diag;
1381 });
1382
1383 bool drop = false;
1384 for (size_t i = 1; i < n; ++i) {
1385 if (!drop) {
1386 auto const& x = drop_vec[i - 1];
1387 auto const& y = drop_vec[i];
1388 auto a = x.val / x.diag;
1389 auto b = y.val / y.diag;
1390 if (realThreshold * realThreshold * a > b) {
1391 drop = true;
1392#ifdef HAVE_MUELU_DEBUG
1393 if (distanceLaplacianCutVerbose) {
1394 std::cout << "DJS: KEEP, N, ROW: " << i + 1 << ", " << n << ", " << row << std::endl;
1395 }
1396#endif
1397 }
1398 }
1399 drop_vec[i].drop = drop;
1400 }
1401 }
1402
1403 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol const& a, DropTol const& b) {
1404 return a.col < b.col;
1405 });
1406
1407 for (LO idxID = 0; idxID < (LO)drop_vec.size(); idxID++) {
1408 LO col = indices[drop_vec[idxID].col];
1409
1410 // don't drop diagonal
1411 if (row == col) {
1412 columns(realnnz++) = col;
1413 rownnz++;
1414 // printf("(%d,%d) KEEP %13s matrix = %6.4e\n",row,row,"DIAGONAL",drop_vec[idxID].aux_val);
1415 continue;
1416 }
1417
1418 if (!drop_vec[idxID].drop) {
1419 columns(realnnz++) = col;
1420 // printf("(%d,%d) KEEP dlap = %6.4e matrix = %6.4e\n",row,col,drop_vec[idxID].val/drop_vec[idxID].diag,drop_vec[idxID].aux_val);
1421 rownnz++;
1422 } else {
1423 // printf("(%d,%d) DROP dlap = %6.4e matrix = %6.4e\n",row,col,drop_vec[idxID].val/drop_vec[idxID].diag,drop_vec[idxID].aux_val);
1424 numDropped++;
1425 }
1426 }
1427 }
1428 } else {
1429 // Skip laplace calculation and threshold comparison for zero threshold
1430 for (LO colID = 0; colID < nnz; colID++) {
1431 LO col = indices[colID];
1432 columns(realnnz++) = col;
1433 rownnz++;
1434 }
1435 }
1436
1437 if (rownnz == 1) {
1438 // If the only element remaining after filtering is diagonal, mark node as boundary
1439 // FIXME: this should really be replaced by the following
1440 // if (indices.size() == 1 && indices[0] == row)
1441 // boundaryNodes[row] = true;
1442 // We do not do it this way now because there is no framework for distinguishing isolated
1443 // and boundary nodes in the aggregation algorithms
1444 amalgBoundaryNodes[row] = true;
1445 }
1446
1447 if (use_stop_array)
1448 rows_stop[row] = rownnz + rows[row];
1449 else
1450 rows[row + 1] = realnnz;
1451 } // for (LO row = 0; row < numRows; row++)
1452
1453 } // subtimer
1454
1455 if (use_stop_array) {
1456 // Do symmetrization of the cut matrix
1457 // NOTE: We assume nested row/column maps here
1458 for (LO row = 0; row < numRows; row++) {
1459 for (LO colidx = rows[row]; colidx < rows_stop[row]; colidx++) {
1460 LO col = columns[colidx];
1461 if (col >= numRows) continue;
1462
1463 bool found = false;
1464 for (LO t_col = rows(col); !found && t_col < rows_stop[col]; t_col++) {
1465 if (columns[t_col] == row)
1466 found = true;
1467 }
1468 // We didn't find the transpose buddy, so let's symmetrize, unless we'd be symmetrizing
1469 // into a Dirichlet unknown. In that case don't.
1470 if (!found && !pointBoundaryNodes[col] && Teuchos::as<typename LWGraph::row_type::value_type>(rows_stop[col]) < rows[col + 1]) {
1471 LO new_idx = rows_stop[col];
1472 // printf("(%d,%d) SYMADD entry\n",col,row);
1473 columns[new_idx] = row;
1474 rows_stop[col]++;
1475 numDropped--;
1476 }
1477 }
1478 }
1479
1480 // Condense everything down
1481 LO current_start = 0;
1482 for (LO row = 0; row < numRows; row++) {
1483 LO old_start = current_start;
1484 for (LO col = rows(row); col < rows_stop[row]; col++) {
1485 if (current_start != col) {
1486 columns(current_start) = columns(col);
1487 }
1488 current_start++;
1489 }
1490 rows[row] = old_start;
1491 }
1492 rows(numRows) = realnnz = current_start;
1493 }
1494
1495 RCP<LWGraph> graph;
1496 {
1497 SubFactoryMonitor m1(*this, "Build amalgamated graph", currentLevel);
1498 graph = rcp(new LWGraph(rows, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), uniqueMap, nonUniqueMap, "amalgamated graph of A"));
1499 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1500 } // subtimer
1501
1502 if (GetVerbLevel() & Statistics1) {
1503 GO numLocalBoundaryNodes = 0;
1504 GO numGlobalBoundaryNodes = 0;
1505
1506 for (size_t i = 0; i < amalgBoundaryNodes.size(); ++i)
1507 if (amalgBoundaryNodes(i))
1508 numLocalBoundaryNodes++;
1509
1510 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
1511 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1512 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " agglomerated Dirichlet nodes"
1513 << " using threshold " << dirichletThreshold << std::endl;
1514 }
1515
1516 Set(currentLevel, "Graph", graph);
1517 Set(currentLevel, "DofsPerNode", blkSize);
1518 }
1519 }
1520
1521 if (GetVerbLevel() & Statistics1) {
1522 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
1523 GO numGlobalTotal, numGlobalDropped;
1524 MueLu_sumAll(comm, numTotal, numGlobalTotal);
1525 MueLu_sumAll(comm, numDropped, numGlobalDropped);
1526 GetOStream(Statistics1) << "Number of dropped entries in " << graphType << " matrix graph: " << numGlobalDropped << "/" << numGlobalTotal;
1527 if (numGlobalTotal != 0)
1528 GetOStream(Statistics1) << " (" << 100 * Teuchos::as<double>(numGlobalDropped) / Teuchos::as<double>(numGlobalTotal) << "%)";
1529 GetOStream(Statistics1) << std::endl;
1530 }
1531
1532 } else {
1533 // what Tobias has implemented
1534
1535 SC threshold = as<SC>(pL.get<double>("aggregation: drop tol"));
1536 // GetOStream(Runtime0) << "algorithm = \"" << algo << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
1537 GetOStream(Runtime0) << "algorithm = \""
1538 << "failsafe"
1539 << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
1540 Set<bool>(currentLevel, "Filtering", (threshold != STS::zero()));
1541
1542 RCP<const Map> rowMap = A->getRowMap();
1543 RCP<const Map> colMap = A->getColMap();
1544
1545 LO blockdim = 1; // block dim for fixed size blocks
1546 GO indexBase = rowMap->getIndexBase(); // index base of maps
1547 GO offset = 0;
1548
1549 // 1) check for blocking/striding information
1550 if (A->IsView("stridedMaps") &&
1551 Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
1552 Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
1553 RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
1554 TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null, Exceptions::BadCast, "MueLu::CoalesceFactory::Build: cast to strided row map failed.");
1555 blockdim = strMap->getFixedBlockSize();
1556 offset = strMap->getOffset();
1557 oldView = A->SwitchToView(oldView);
1558 GetOStream(Statistics1) << "CoalesceDropFactory::Build():"
1559 << " found blockdim=" << blockdim << " from strided maps. offset=" << offset << std::endl;
1560 } else
1561 GetOStream(Statistics1) << "CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
1562
1563 // 2) get row map for amalgamated matrix (graph of A)
1564 // with same distribution over all procs as row map of A
1565 RCP<const Map> nodeMap = amalInfo->getNodeRowMap();
1566 GetOStream(Statistics1) << "CoalesceDropFactory: nodeMap " << nodeMap->getLocalNumElements() << "/" << nodeMap->getGlobalNumElements() << " elements" << std::endl;
1567
1568 // 3) create graph of amalgamated matrix
1569 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getLocalMaxNumRowEntries() * blockdim);
1570
1571 LO numRows = A->getRowMap()->getLocalNumElements();
1572 LO numNodes = nodeMap->getLocalNumElements();
1573 typename LWGraph::boundary_nodes_type amalgBoundaryNodes("amalgBoundaryNodes", numNodes);
1574 Kokkos::deep_copy(amalgBoundaryNodes, false);
1575 const ArrayRCP<int> numberDirichletRowsPerNode(numNodes, 0); // helper array counting the number of Dirichlet nodes associated with node
1576 bool bIsDiagonalEntry = false; // boolean flag stating that grid==gcid
1577
1578 // 4) do amalgamation. generate graph of amalgamated matrix
1579 // Note, this code is much more inefficient than the leightwight implementation
1580 // Most of the work has already been done in the AmalgamationFactory
1581 for (LO row = 0; row < numRows; row++) {
1582 // get global DOF id
1583 GO grid = rowMap->getGlobalElement(row);
1584
1585 // reinitialize boolean helper variable
1586 bIsDiagonalEntry = false;
1587
1588 // translate grid to nodeid
1589 GO nodeId = AmalgamationFactory::DOFGid2NodeId(grid, blockdim, offset, indexBase);
1590
1591 size_t nnz = A->getNumEntriesInLocalRow(row);
1592 Teuchos::ArrayView<const LO> indices;
1593 Teuchos::ArrayView<const SC> vals;
1594 A->getLocalRowView(row, indices, vals);
1595
1596 RCP<std::vector<GO>> cnodeIds = Teuchos::rcp(new std::vector<GO>); // global column block ids
1597 LO realnnz = 0;
1598 for (LO col = 0; col < Teuchos::as<LO>(nnz); col++) {
1599 GO gcid = colMap->getGlobalElement(indices[col]); // global column id
1600
1601 if (vals[col] != STS::zero()) {
1602 GO cnodeId = AmalgamationFactory::DOFGid2NodeId(gcid, blockdim, offset, indexBase);
1603 cnodeIds->push_back(cnodeId);
1604 realnnz++; // increment number of nnz in matrix row
1605 if (grid == gcid) bIsDiagonalEntry = true;
1606 }
1607 }
1608
1609 if (realnnz == 1 && bIsDiagonalEntry == true) {
1610 LO lNodeId = nodeMap->getLocalElement(nodeId);
1611 numberDirichletRowsPerNode[lNodeId] += 1; // increment Dirichlet row counter associated with lNodeId
1612 if (numberDirichletRowsPerNode[lNodeId] == blockdim) // mark full Dirichlet nodes
1613 amalgBoundaryNodes[lNodeId] = true;
1614 }
1615
1616 Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp(cnodeIds);
1617
1618 if (arr_cnodeIds.size() > 0)
1619 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
1620 }
1621 // fill matrix graph
1622 crsGraph->fillComplete(nodeMap, nodeMap);
1623
1624 // 5) create MueLu Graph object
1625 RCP<LWGraph> graph = rcp(new LWGraph(crsGraph, "amalgamated graph of A"));
1626
1627 // Detect and record rows that correspond to Dirichlet boundary conditions
1628 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1629
1630 if (GetVerbLevel() & Statistics1) {
1631 GO numLocalBoundaryNodes = 0;
1632 GO numGlobalBoundaryNodes = 0;
1633 for (size_t i = 0; i < amalgBoundaryNodes.size(); ++i)
1634 if (amalgBoundaryNodes(i))
1635 numLocalBoundaryNodes++;
1636 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
1637 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1638 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
1639 }
1640
1641 // 6) store results in Level
1642 // graph->SetBoundaryNodeMap(gBoundaryNodeMap);
1643 Set(currentLevel, "DofsPerNode", blockdim);
1644 Set(currentLevel, "Graph", graph);
1645
1646 } // if (doExperimentalWrap) ... else ...
1647
1648} // Build
1649
1650template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1651void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MergeRows(const Matrix& A, const LO row, Array<LO>& cols, const Array<LO>& translation) const {
1652 typedef typename ArrayView<const LO>::size_type size_type;
1653
1654 // extract striding information
1655 LO blkSize = A.GetFixedBlockSize(); //< stores the size of the block within the strided map
1656 if (A.IsView("stridedMaps") == true) {
1657 Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
1658 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
1659 TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
1660 if (strMap->getStridedBlockId() > -1)
1661 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1662 }
1663
1664 // count nonzero entries in all dof rows associated with node row
1665 size_t nnz = 0, pos = 0;
1666 for (LO j = 0; j < blkSize; j++)
1667 nnz += A.getNumEntriesInLocalRow(row * blkSize + j);
1668
1669 if (nnz == 0) {
1670 cols.resize(0);
1671 return;
1672 }
1673
1674 cols.resize(nnz);
1675
1676 // loop over all local dof rows associated with local node "row"
1677 ArrayView<const LO> inds;
1678 ArrayView<const SC> vals;
1679 for (LO j = 0; j < blkSize; j++) {
1680 A.getLocalRowView(row * blkSize + j, inds, vals);
1681 size_type numIndices = inds.size();
1682
1683 if (numIndices == 0) // skip empty dof rows
1684 continue;
1685
1686 // cols: stores all local node ids for current local node id "row"
1687 cols[pos++] = translation[inds[0]];
1688 for (size_type k = 1; k < numIndices; k++) {
1689 LO nodeID = translation[inds[k]];
1690 // Here we try to speed up the process by reducing the size of an array
1691 // to sort. This works if the column nonzeros belonging to the same
1692 // node are stored consequently.
1693 if (nodeID != cols[pos - 1])
1694 cols[pos++] = nodeID;
1695 }
1696 }
1697 cols.resize(pos);
1698 nnz = pos;
1699
1700 // Sort and remove duplicates
1701 std::sort(cols.begin(), cols.end());
1702 pos = 0;
1703 for (size_t j = 1; j < nnz; j++)
1704 if (cols[j] != cols[pos])
1705 cols[++pos] = cols[j];
1706 cols.resize(pos + 1);
1707}
1708
1709template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1710void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MergeRowsWithDropping(const Matrix& A, const LO row, const ArrayRCP<const SC>& ghostedDiagVals, SC threshold, Array<LO>& cols, const Array<LO>& translation) const {
1711 typedef typename ArrayView<const LO>::size_type size_type;
1712 typedef Teuchos::ScalarTraits<SC> STS;
1713
1714 // extract striding information
1715 LO blkSize = A.GetFixedBlockSize(); //< stores the size of the block within the strided map
1716 if (A.IsView("stridedMaps") == true) {
1717 Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
1718 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
1719 TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
1720 if (strMap->getStridedBlockId() > -1)
1721 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1722 }
1723
1724 // count nonzero entries in all dof rows associated with node row
1725 size_t nnz = 0, pos = 0;
1726 for (LO j = 0; j < blkSize; j++)
1727 nnz += A.getNumEntriesInLocalRow(row * blkSize + j);
1728
1729 if (nnz == 0) {
1730 cols.resize(0);
1731 return;
1732 }
1733
1734 cols.resize(nnz);
1735
1736 // loop over all local dof rows associated with local node "row"
1737 ArrayView<const LO> inds;
1738 ArrayView<const SC> vals;
1739 for (LO j = 0; j < blkSize; j++) {
1740 A.getLocalRowView(row * blkSize + j, inds, vals);
1741 size_type numIndices = inds.size();
1742
1743 if (numIndices == 0) // skip empty dof rows
1744 continue;
1745
1746 // cols: stores all local node ids for current local node id "row"
1747 LO prevNodeID = -1;
1748 for (size_type k = 0; k < numIndices; k++) {
1749 LO dofID = inds[k];
1750 LO nodeID = translation[inds[k]];
1751
1752 // we avoid a square root by using squared values
1753 typename STS::magnitudeType aiiajj = STS::magnitude(threshold * threshold * ghostedDiagVals[dofID] * ghostedDiagVals[row * blkSize + j]); // eps^2 * |a_ii| * |a_jj|
1754 typename STS::magnitudeType aij = STS::magnitude(vals[k] * vals[k]);
1755
1756 // check dropping criterion
1757 if (aij > aiiajj || (row * blkSize + j == dofID)) {
1758 // accept entry in graph
1759
1760 // Here we try to speed up the process by reducing the size of an array
1761 // to sort. This works if the column nonzeros belonging to the same
1762 // node are stored consequently.
1763 if (nodeID != prevNodeID) {
1764 cols[pos++] = nodeID;
1765 prevNodeID = nodeID;
1766 }
1767 }
1768 }
1769 }
1770 cols.resize(pos);
1771 nnz = pos;
1772
1773 // Sort and remove duplicates
1774 std::sort(cols.begin(), cols.end());
1775 pos = 0;
1776 for (size_t j = 1; j < nnz; j++)
1777 if (cols[j] != cols[pos])
1778 cols[++pos] = cols[j];
1779 cols.resize(pos + 1);
1780
1781 return;
1782}
1783
1784template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1785Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BlockDiagonalize(Level& currentLevel, const RCP<Matrix>& A, bool generate_matrix) const {
1786 typedef Teuchos::ScalarTraits<SC> STS;
1787
1788 const ParameterList& pL = GetParameterList();
1789 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<double>("aggregation: Dirichlet threshold")));
1790 const typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<double>("aggregation: row sum drop tol"));
1791
1792 RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector>>(currentLevel, "BlockNumber");
1793 RCP<LocalOrdinalVector> ghostedBlockNumber;
1794 GetOStream(Statistics1) << "Using BlockDiagonal Graph before dropping (with provided blocking)" << std::endl;
1795
1796 // Ghost the column block numbers if we need to
1797 RCP<const Import> importer = A->getCrsGraph()->getImporter();
1798 if (!importer.is_null()) {
1799 SubFactoryMonitor m1(*this, "Block Number import", currentLevel);
1800 ghostedBlockNumber = Xpetra::VectorFactory<LO, LO, GO, NO>::Build(importer->getTargetMap());
1801 ghostedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
1802 } else {
1803 ghostedBlockNumber = BlockNumber;
1804 }
1805
1806 // Accessors for block numbers
1807 Teuchos::ArrayRCP<const LO> row_block_number = BlockNumber->getData(0);
1808 Teuchos::ArrayRCP<const LO> col_block_number = ghostedBlockNumber->getData(0);
1809
1810 // allocate space for the local graph
1811 typename CrsMatrix::local_matrix_device_type::row_map_type::host_mirror_type::non_const_type rows_mat;
1812 typename LWGraph::row_type::non_const_type rows_graph;
1813 typename LWGraph::entries_type::non_const_type columns;
1814 typename CrsMatrix::local_matrix_device_type::values_type::host_mirror_type::non_const_type values;
1815 RCP<CrsMatrixWrap> crs_matrix_wrap;
1816
1817 if (generate_matrix) {
1818 crs_matrix_wrap = rcp(new CrsMatrixWrap(A->getRowMap(), A->getColMap(), 0));
1819 rows_mat = typename CrsMatrix::local_matrix_device_type::row_map_type::host_mirror_type::non_const_type("rows_mat", A->getLocalNumRows() + 1);
1820 } else {
1821 rows_graph = typename LWGraph::row_type::non_const_type("rows_graph", A->getLocalNumRows() + 1);
1822 }
1823 columns = typename LWGraph::entries_type::non_const_type("columns", A->getLocalNumEntries());
1824 values = typename CrsMatrix::local_matrix_device_type::values_type::host_mirror_type::non_const_type("values", A->getLocalNumEntries());
1825
1826 LO realnnz = 0;
1827 GO numDropped = 0, numTotal = 0;
1828 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
1829 LO row_block = row_block_number[row];
1830 size_t nnz = A->getNumEntriesInLocalRow(row);
1831 ArrayView<const LO> indices;
1832 ArrayView<const SC> vals;
1833 A->getLocalRowView(row, indices, vals);
1834
1835 LO rownnz = 0;
1836 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
1837 LO col = indices[colID];
1838 LO col_block = col_block_number[col];
1839
1840 if (row_block == col_block) {
1841 if (generate_matrix) values[realnnz] = vals[colID];
1842 columns[realnnz++] = col;
1843 rownnz++;
1844 } else
1845 numDropped++;
1846 }
1847 if (generate_matrix)
1848 rows_mat[row + 1] = realnnz;
1849 else
1850 rows_graph[row + 1] = realnnz;
1851 }
1852
1853 auto boundaryNodes = MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows_kokkos_host(*A, dirichletThreshold);
1854 if (rowSumTol > 0.)
1855 Utilities::ApplyRowSumCriterionHost(*A, rowSumTol, boundaryNodes);
1856
1857 numTotal = A->getLocalNumEntries();
1858
1859 if (GetVerbLevel() & Statistics1) {
1860 GO numLocalBoundaryNodes = 0;
1861 GO numGlobalBoundaryNodes = 0;
1862 for (size_t i = 0; i < boundaryNodes.size(); ++i)
1863 if (boundaryNodes(i))
1864 numLocalBoundaryNodes++;
1865 RCP<const Teuchos::Comm<int>> comm = A->getRowMap()->getComm();
1866 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1867 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
1868
1869 GO numGlobalTotal, numGlobalDropped;
1870 MueLu_sumAll(comm, numTotal, numGlobalTotal);
1871 MueLu_sumAll(comm, numDropped, numGlobalDropped);
1872 GetOStream(Statistics1) << "Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped << "/" << numGlobalTotal;
1873 if (numGlobalTotal != 0)
1874 GetOStream(Statistics1) << " (" << 100 * Teuchos::as<double>(numGlobalDropped) / Teuchos::as<double>(numGlobalTotal) << "%)";
1875 GetOStream(Statistics1) << std::endl;
1876 }
1877
1878 Set(currentLevel, "Filtering", true);
1879
1880 if (generate_matrix) {
1881 // NOTE: Trying to use A's Import/Export objects will cause the code to segfault back in Build() with errors on the Import
1882 // if you're using Epetra. I'm not really sure why. By using the Col==Domain and Row==Range maps, we get null Import/Export objects
1883 // here, which is legit, because we never use them anyway.
1884 if constexpr (std::is_same<typename LWGraph::row_type,
1885 typename CrsMatrix::local_matrix_device_type::row_map_type>::value) {
1886 crs_matrix_wrap->getCrsMatrix()->setAllValues(rows_mat, columns, values);
1887 } else {
1888 auto rows_mat2 = typename CrsMatrix::local_matrix_device_type::row_map_type::non_const_type("rows_mat2", rows_mat.extent(0));
1889 Kokkos::deep_copy(rows_mat2, rows_mat);
1890 auto columns2 = typename CrsMatrix::local_graph_device_type::entries_type::non_const_type("columns2", columns.extent(0));
1891 Kokkos::deep_copy(columns2, columns);
1892 auto values2 = typename CrsMatrix::local_matrix_device_type::values_type::non_const_type("values2", values.extent(0));
1893 Kokkos::deep_copy(values2, values);
1894 crs_matrix_wrap->getCrsMatrix()->setAllValues(rows_mat2, columns2, values2);
1895 }
1896 crs_matrix_wrap->getCrsMatrix()->expertStaticFillComplete(A->getColMap(), A->getRowMap());
1897 } else {
1898 RCP<LWGraph> graph = rcp(new LWGraph(rows_graph, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), A->getRowMap(), A->getColMap(), "block-diagonalized graph of A"));
1899 graph->SetBoundaryNodeMap(boundaryNodes);
1900 Set(currentLevel, "Graph", graph);
1901 }
1902
1903 Set(currentLevel, "DofsPerNode", 1);
1904 return crs_matrix_wrap;
1905}
1906
1907template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1908void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BlockDiagonalizeGraph(const RCP<LWGraph>& inputGraph, const RCP<LocalOrdinalVector>& ghostedBlockNumber, RCP<LWGraph>& outputGraph, RCP<const Import>& importer) const {
1909 TEUCHOS_TEST_FOR_EXCEPTION(ghostedBlockNumber.is_null(), Exceptions::RuntimeError, "BlockDiagonalizeGraph(): ghostedBlockNumber is null.");
1910 const ParameterList& pL = GetParameterList();
1911
1912 const bool localizeColoringGraph = pL.get<bool>("aggregation: coloring: localize color graph");
1913
1914 GetOStream(Statistics1) << "Using BlockDiagonal Graph after Dropping (with provided blocking)";
1915 if (localizeColoringGraph)
1916 GetOStream(Statistics1) << ", with localization" << std::endl;
1917 else
1918 GetOStream(Statistics1) << ", without localization" << std::endl;
1919
1920 // Accessors for block numbers
1921 Teuchos::ArrayRCP<const LO> row_block_number = ghostedBlockNumber->getData(0);
1922 Teuchos::ArrayRCP<const LO> col_block_number = ghostedBlockNumber->getData(0);
1923
1924 // allocate space for the local graph
1925 ArrayRCP<size_t> rows_mat;
1926 typename LWGraph::row_type::non_const_type rows_graph("rows_graph", inputGraph->GetNodeNumVertices() + 1);
1927 typename LWGraph::entries_type::non_const_type columns("columns", inputGraph->GetNodeNumEdges());
1928
1929 LO realnnz = 0;
1930 GO numDropped = 0, numTotal = 0;
1931 const LO numRows = Teuchos::as<LO>(inputGraph->GetDomainMap()->getLocalNumElements());
1932 if (localizeColoringGraph) {
1933 for (LO row = 0; row < numRows; ++row) {
1934 LO row_block = row_block_number[row];
1935 auto indices = inputGraph->getNeighborVertices(row);
1936
1937 LO rownnz = 0;
1938 for (LO colID = 0; colID < Teuchos::as<LO>(indices.length); colID++) {
1939 LO col = indices(colID);
1940 LO col_block = col_block_number[col];
1941
1942 if ((row_block == col_block) && (col < numRows)) {
1943 columns(realnnz++) = col;
1944 rownnz++;
1945 } else
1946 numDropped++;
1947 }
1948 rows_graph(row + 1) = realnnz;
1949 }
1950 } else {
1951 // ghosting of boundary node map
1952 auto boundaryNodes = inputGraph->GetBoundaryNodeMap();
1953 auto boundaryNodesVector = Xpetra::VectorFactory<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>::Build(inputGraph->GetDomainMap());
1954 for (size_t i = 0; i < inputGraph->GetNodeNumVertices(); i++)
1955 boundaryNodesVector->getDataNonConst(0)[i] = boundaryNodes[i];
1956 // Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write("boundary",*boundaryNodesVector);
1957 auto boundaryColumnVector = Xpetra::VectorFactory<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>::Build(inputGraph->GetImportMap());
1958 boundaryColumnVector->doImport(*boundaryNodesVector, *importer, Xpetra::INSERT);
1959 auto boundaryColumn = boundaryColumnVector->getData(0);
1960
1961 for (LO row = 0; row < numRows; ++row) {
1962 LO row_block = row_block_number[row];
1963 auto indices = inputGraph->getNeighborVertices(row);
1964
1965 LO rownnz = 0;
1966 for (LO colID = 0; colID < Teuchos::as<LO>(indices.length); colID++) {
1967 LO col = indices(colID);
1968 LO col_block = col_block_number[col];
1969
1970 if ((row_block == col_block) && ((row == col) || (boundaryColumn[col] == 0))) {
1971 columns(realnnz++) = col;
1972 rownnz++;
1973 } else
1974 numDropped++;
1975 }
1976 rows_graph(row + 1) = realnnz;
1977 }
1978 }
1979
1980 numTotal = inputGraph->GetNodeNumEdges();
1981
1982 if (GetVerbLevel() & Statistics1) {
1983 RCP<const Teuchos::Comm<int>> comm = inputGraph->GetDomainMap()->getComm();
1984 GO numGlobalTotal, numGlobalDropped;
1985 MueLu_sumAll(comm, numTotal, numGlobalTotal);
1986 MueLu_sumAll(comm, numDropped, numGlobalDropped);
1987 GetOStream(Statistics1) << "Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped << "/" << numGlobalTotal;
1988 if (numGlobalTotal != 0)
1989 GetOStream(Statistics1) << " (" << 100 * Teuchos::as<double>(numGlobalDropped) / Teuchos::as<double>(numGlobalTotal) << "%)";
1990 GetOStream(Statistics1) << std::endl;
1991 }
1992
1993 if (localizeColoringGraph) {
1994 outputGraph = rcp(new LWGraph(rows_graph, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), inputGraph->GetDomainMap(), inputGraph->GetImportMap(), "block-diagonalized graph of A"));
1995 outputGraph->SetBoundaryNodeMap(inputGraph->GetBoundaryNodeMap());
1996 } else {
1997 TEUCHOS_ASSERT(inputGraph->GetDomainMap()->lib() == Xpetra::UseTpetra);
1998 auto outputGraph2 = rcp(new LWGraph(rows_graph, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), inputGraph->GetDomainMap(), inputGraph->GetImportMap(), "block-diagonalized graph of A"));
1999
2000 auto tpGraph = Xpetra::toTpetra(rcp_const_cast<const CrsGraph>(outputGraph2->GetCrsGraph()));
2001 auto sym = rcp(new Tpetra::CrsGraphTransposer<LocalOrdinal, GlobalOrdinal, Node>(tpGraph));
2002 auto tpGraphSym = sym->symmetrize();
2003 auto lclGraphSym = tpGraphSym->getLocalGraphHost();
2004 auto colIndsSym = lclGraphSym.entries;
2005
2006 auto rowsSym = tpGraphSym->getLocalRowPtrsHost();
2007 typename LWGraph::row_type::non_const_type rows_graphSym("rows_graphSym", rowsSym.size());
2008 for (size_t row = 0; row < rowsSym.size(); row++)
2009 rows_graphSym(row) = rowsSym(row);
2010 outputGraph = rcp(new LWGraph(rows_graphSym, colIndsSym, inputGraph->GetDomainMap(), Xpetra::toXpetra(tpGraphSym->getColMap()), "block-diagonalized graph of A"));
2011 outputGraph->SetBoundaryNodeMap(inputGraph->GetBoundaryNodeMap());
2012 }
2013}
2014
2015} // namespace MueLu
2016
2017#endif // MUELU_COALESCEDROPFACTORY_DEF_HPP
SET_VALID_ENTRY
#define SET_VALID_ENTRY(name)
MueLu_CoalesceDropFactory_decl.hpp
MueLu_Exceptions.hpp
MueLu_Level.hpp
MueLu_MasterList.hpp
MueLu_Monitor.hpp
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition MueLu_PerfUtils_decl.hpp:25
rows
rowsType rows
Definition MueLu_TentativePFactory_kokkos_def.hpp:84
MueLu::AmalgamationFactory::AmalgamateMap
static void AmalgamateMap(const Map &sourceMap, const Matrix &A, RCP< const Map > &amalgamatedMap, Array< LO > &translation)
Method to create merged map for systems of PDEs.
Definition MueLu_AmalgamationFactory_def.hpp:153
MueLu::AmalgamationFactory::DOFGid2NodeId
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
Translate global (row/column) id to global amalgamation block id.
Definition MueLu_AmalgamationFactory_def.hpp:224
MueLu::CoalesceDropFactory::MergeRows
void MergeRows(const Matrix &A, const LO row, Array< LO > &cols, const Array< LO > &translation) const
Method to merge rows of matrix for systems of PDEs.
Definition MueLu_CoalesceDropFactory_def.hpp:1651
MueLu::CoalesceDropFactory::BlockDiagonalize
Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > BlockDiagonalize(Level &currentLevel, const RCP< Matrix > &A, bool generate_matrix) const
Definition MueLu_CoalesceDropFactory_def.hpp:1785
MueLu::CoalesceDropFactory::DeclareInput
void DeclareInput(Level &currentLevel) const
Input.
Definition MueLu_CoalesceDropFactory_def.hpp:124
MueLu::CoalesceDropFactory::BlockDiagonalizeGraph
void BlockDiagonalizeGraph(const RCP< LWGraph > &inputGraph, const RCP< LocalOrdinalVector > &ghostedBlockNumber, RCP< LWGraph > &outputGraph, RCP< const Import > &importer) const
Definition MueLu_CoalesceDropFactory_def.hpp:1908
MueLu::CoalesceDropFactory::MergeRowsWithDropping
void MergeRowsWithDropping(const Matrix &A, const LO row, const ArrayRCP< const SC > &ghostedDiagVals, SC threshold, Array< LO > &cols, const Array< LO > &translation) const
Definition MueLu_CoalesceDropFactory_def.hpp:1710
MueLu::CoalesceDropFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition MueLu_CoalesceDropFactory_def.hpp:87
MueLu::CoalesceDropFactory::Build
void Build(Level &currentLevel) const
Build an object with this factory.
Definition MueLu_CoalesceDropFactory_def.hpp:143
MueLu::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition MueLu_Exceptions.hpp:21
MueLu::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition MueLu_Exceptions.hpp:50
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition MueLu_Exceptions.hpp:36
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition MueLu_Monitor.hpp:165
MueLu::LWGraphBase::boundary_nodes_type
Kokkos::View< bool *, memory_space > boundary_nodes_type
Definition MueLu_LWGraphBase.hpp:91
MueLu::LWGraphBase::row_type
typename local_graph_type::row_map_type row_type
Definition MueLu_LWGraphBase.hpp:92
MueLu::LWGraph
Lightweight MueLu representation of a compressed row storage graph.
Definition MueLu_LWGraph_decl.hpp:33
MueLu::Level
Class that holds all level-specific information.
Definition MueLu_Level.hpp:63
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition MueLu_Level.cpp:52
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition MueLu_Monitor.hpp:214
MueLu::UtilitiesBase::DetectDirichletRows_kokkos_host
static Kokkos::View< bool *, typename Kokkos::HostSpace > DetectDirichletRows_kokkos_host(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
Definition MueLu_UtilitiesBase_def.hpp:1203
MueLu::UtilitiesBase::GetMatrixMaxMinusOffDiagonal
static Teuchos::RCP< Vector > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i\not=k}(-a_ik), for each for i in the matrix.
Definition MueLu_UtilitiesBase_def.hpp:622
MueLu::UtilitiesBase::ApplyRowSumCriterion
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Apply Rowsum Criterion.
MueLu::UtilitiesBase::Distance2
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar > > &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
Definition MueLu_UtilitiesBase_def.hpp:1019
MueLu::UtilitiesBase::ApplyRowSumCriterionHost
static void ApplyRowSumCriterionHost(const Matrix &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, Kokkos::HostSpace > &dirichletRows)
MueLu::UtilitiesBase::GetMatrixOverlappedDiagonal
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
Definition MueLu_UtilitiesBase_def.hpp:770
MueLu::Utilities
MueLu utility class.
Definition MueLu_Utilities_decl.hpp:70
MueLu
Namespace for MueLu classes and methods.
Definition MueLu_BrickAggregationFactory_decl.hpp:41
MueLu::Warnings0
@ Warnings0
Important warning messages (one line)
Definition MueLu_VerbosityLevel.hpp:20
MueLu::Statistics1
@ Statistics1
Print more statistics.
Definition MueLu_VerbosityLevel.hpp:34
MueLu::Runtime0
@ Runtime0
One-liner description of what is happening.
Definition MueLu_VerbosityLevel.hpp:25
MueLu::Warnings1
@ Warnings1
Additional warnings.
Definition MueLu_VerbosityLevel.hpp:22
MueLu::Parameters0
@ Parameters0
Print class parameters.
Definition MueLu_VerbosityLevel.hpp:30
MueLu::Details::DropTol::operator=
DropTol & operator=(DropTol const &)=default
MueLu::Details::DropTol::DropTol
DropTol(DropTol &&)=default
MueLu::Details::DropTol::DropTol
DropTol(real_type val_, real_type diag_, LO col_, bool drop_)
Definition MueLu_CoalesceDropFactory_def.hpp:65
MueLu::Details::DropTol::operator=
DropTol & operator=(DropTol &&)=default
MueLu::Details::DropTol::DropTol
DropTol(DropTol const &)=default