MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_BrickAggregationFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_
47 #define MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_
48 
50 #ifdef HAVE_MPI
52 #include <Teuchos_CommHelpers.hpp>
53 #endif
55 
56 #include <Xpetra_Import.hpp>
57 #include <Xpetra_ImportFactory.hpp>
58 #include <Xpetra_Map.hpp>
59 #include <Xpetra_MapFactory.hpp>
60 #include <Xpetra_Matrix.hpp>
61 #include <Xpetra_MultiVector.hpp>
62 #include <Xpetra_MultiVectorFactory.hpp>
63 
64 #include "MueLu_Aggregates.hpp"
65 #include "MueLu_Level.hpp"
66 #include "MueLu_MasterList.hpp"
67 #include "MueLu_Monitor.hpp"
68 #include "MueLu_Utilities.hpp"
69 #include "MueLu_LWGraph.hpp"
70 
71 #include "MueLu_LWGraph.hpp"
72 
73 namespace MueLu {
74 
75 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77  RCP<ParameterList> validParamList = rcp(new ParameterList());
78 
79 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
80  SET_VALID_ENTRY("aggregation: brick x size");
81  SET_VALID_ENTRY("aggregation: brick y size");
82  SET_VALID_ENTRY("aggregation: brick z size");
83  SET_VALID_ENTRY("aggregation: brick x Dirichlet");
84  SET_VALID_ENTRY("aggregation: brick y Dirichlet");
85  SET_VALID_ENTRY("aggregation: brick z Dirichlet");
86 #undef SET_VALID_ENTRY
87 
88  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory for matrix");
89  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for coordinates");
90  return validParamList;
91 }
92 
93 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
95  Input(currentLevel, "A");
96  Input(currentLevel, "Coordinates");
97 }
98 
99 // The current implementation cannot deal with bricks larger than 3x3(x3) in
100 // parallel. The reason is that aggregation infrastructure in place has
101 // major drawbacks.
102 //
103 // Aggregates class is constructed with a help of a provided map, either
104 // taken from a graph, or provided directly. This map is usually taken to be
105 // a column map of a matrix. The reason for that is that if we have an
106 // overlapped aggregation, we want the processor owning aggregates to store
107 // agg id for all nodes in this aggregate. If we used row map, there would
108 // be no way for the processor to know whether there are some other nodes on
109 // a different processor which belong to its aggregate. On the other hand,
110 // using column map allows both vertex2AggId and procWinner arrays in
111 // Aggregates class to store some extra data, such as whether nodes belonging
112 // to a different processor belong to this processor aggregate.
113 //
114 // The drawback of this is that it stores only overlap=1 data. For aggressive
115 // coarsening, such a brick aggregation with a large single dimension of
116 // brick, it could happen that we need to know depth two or more extra nodes
117 // in the other processor subdomain.
118 //
119 // Another issue is that we may have some implicit connection between
120 // aggregate map and maps of A used in the construction of a tentative
121 // prolongator.
122 //
123 // Another issue is that it seems that some info is unused or not required.
124 // Specifically, it seems that if a node belongs to an aggregate on a
125 // different processor, we don't actually need to set vertex2AggId and
126 // procWinner, despite the following comment in
127 // Aggregates decl:
128 // vertex2AggId[k] gives a local id
129 // corresponding to the aggregate to which
130 // local id k has been assigned. While k
131 // is the local id on my processor (MyPID)
132 // vertex2AggId[k] is the local id on the
133 // processor which actually owns the
134 // aggregate. This owning processor has id
135 // given by procWinner[k].
136 // It is possible that that info is only used during arbitration in
137 // CoupledAggregationFactory.
138 //
139 // The steps that we need to do to resolve this issue:
140 // - Break the link between maps in TentativePFactory, allowing any maps in Aggregates
141 // - Allow Aggregates to construct their own maps, if necessary, OR
142 // - construct aggregates based on row map
143 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
145  FactoryMonitor m(*this, "Build", currentLevel);
146 
148 
149  const ParameterList& pL = GetParameterList();
150  RCP<MultiVector_d> coords = Get<RCP<MultiVector_d> >(currentLevel, "Coordinates");
151  RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel, "A");
152  RCP<const Map> rowMap = A->getRowMap();
153  RCP<const Map> colMap = A->getColMap();
154  GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
155 
156  RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
157  int numProcs = comm->getSize();
158  int myRank = comm->getRank();
159 
160  int numPoints = colMap->getLocalNumElements();
161 
162  bx_ = pL.get<int>("aggregation: brick x size");
163  by_ = pL.get<int>("aggregation: brick y size");
164  bz_ = pL.get<int>("aggregation: brick z size");
165 
166  dirichletX_ = pL.get<bool>("aggregation: brick x Dirichlet");
167  dirichletY_ = pL.get<bool>("aggregation: brick y Dirichlet");
168  dirichletZ_ = pL.get<bool>("aggregation: brick z Dirichlet");
169  if (dirichletX_) GetOStream(Runtime0) << "Dirichlet boundaries in the x direction" << std::endl;
170  if (dirichletY_) GetOStream(Runtime0) << "Dirichlet boundaries in the y direction" << std::endl;
171  if (dirichletZ_) GetOStream(Runtime0) << "Dirichlet boundaries in the z direction" << std::endl;
172 
173  if (numProcs > 1) {
174  // TODO: deal with block size > 1 (see comments above)
175  // TEUCHOS_TEST_FOR_EXCEPTION(bx_ > 3 || by_ > 3 || bz_ > 3, Exceptions::RuntimeError, "Currently cannot deal with brick size > 3");
176  }
177 
178  RCP<MultiVector_d> overlappedCoords = coords;
179  RCP<const Import> importer = ImportFactory::Build(coords->getMap(), colMap);
180  if (!importer.is_null()) {
181  overlappedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO>::Build(colMap, coords->getNumVectors());
182  overlappedCoords->doImport(*coords, *importer, Xpetra::INSERT);
183  }
184 
185  // Setup misc structures
186  // Logically, we construct enough data to query topological information of a rectangular grid
187  Setup(comm, overlappedCoords, colMap);
188 
189  GetOStream(Runtime0) << "Using brick size: " << bx_
190  << (nDim_ > 1 ? "x " + toString(by_) : "")
191  << (nDim_ > 2 ? "x " + toString(bz_) : "") << std::endl;
192 
193  // Build the graph
194  BuildGraph(currentLevel, A);
195 
196  // Construct aggregates
197  RCP<Aggregates> aggregates = rcp(new Aggregates(colMap));
198  aggregates->setObjectLabel("Brick");
199 
200  ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
201  ArrayRCP<LO> procWinner = aggregates->GetProcWinner()->getDataNonConst(0);
202 
203  // In the first pass, we set a mapping from a vertex to aggregate global id. We deal with a structured
204  // rectangular mesh, therefore we know the structure of aggregates. For each vertex we can tell exactly
205  // which aggregate it belongs to.
206  // If we determine that the aggregate does not belong to us (i.e. the root vertex does not belong to this
207  // processor, or is outside and we lost "" arbitration), we record the global aggregate id in order to
208  // fetch the local info from the processor owning the aggregate. This is required for aggregates, as it
209  // uses the local aggregate ids of the owning processor.
210  std::set<GO> myAggGIDs, remoteAggGIDs;
211  for (LO LID = 0; LID < numPoints; LID++) {
212  GO aggGID = getAggGID(LID);
213  // printf("[%d] (%d,%d,%d) => agg %d\n",LID,(int)(*xMap_)[x_[LID]],nDim_ > 1 ? (int)(*yMap_)[y_[LID]] : -1,nDim_ > 2 ? (int)(*zMap_)[z_[LID]] : -1,(int)aggGID);
214  if (aggGID == GO_INVALID) continue;
215  // printf("[%d] getRoot = %d\n",(int)LID,(int)getRoot(LID));
216 
217  if ((revMap_.find(getRoot(LID)) != revMap_.end()) && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
218  // Root of the brick aggregate containing GID (<- LID) belongs to us
219  vertex2AggId[LID] = aggGID;
220  myAggGIDs.insert(aggGID);
221 
222  if (isRoot(LID))
223  aggregates->SetIsRoot(LID);
224  // printf("[%d] initial vertex2AggId = %d\n",(int)LID,(int)vertex2AggId[LID]);
225  } else {
226  remoteAggGIDs.insert(aggGID);
227  }
228  }
229  size_t numAggregates = myAggGIDs.size();
230  size_t numRemote = remoteAggGIDs.size();
231  aggregates->SetNumAggregates(numAggregates);
232 
233  std::map<GO, LO> AggG2L; // Map: Agg GID -> Agg LID (possibly on a different processor)
234  std::map<GO, int> AggG2R; // Map: Agg GID -> processor rank owning aggregate
235 
236  Array<GO> myAggGIDsArray(numAggregates), remoteAggGIDsArray(numRemote);
237 
238  // Fill in the maps for aggregates that we own
239  size_t ind = 0;
240  for (typename std::set<GO>::const_iterator it = myAggGIDs.begin(); it != myAggGIDs.end(); it++) {
241  AggG2L[*it] = ind;
242  AggG2R[*it] = myRank;
243 
244  myAggGIDsArray[ind++] = *it;
245  }
246 
247  // The map is a convenient way to fetch remote local indices from global indices.
248  RCP<Map> aggMap = MapFactory::Build(rowMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
249  myAggGIDsArray, 0, comm);
250 
251  ind = 0;
252  for (typename std::set<GO>::const_iterator it = remoteAggGIDs.begin(); it != remoteAggGIDs.end(); it++)
253  remoteAggGIDsArray[ind++] = *it;
254 
255  // Fetch the required aggregate local ids and ranks
256  Array<int> remoteProcIDs(numRemote);
257  Array<LO> remoteLIDs(numRemote);
258  aggMap->getRemoteIndexList(remoteAggGIDsArray, remoteProcIDs, remoteLIDs);
259 
260  // Fill in the maps for aggregates that we don't own but which have some of our vertices
261  for (size_t i = 0; i < numRemote; i++) {
262  AggG2L[remoteAggGIDsArray[i]] = remoteLIDs[i];
263  AggG2R[remoteAggGIDsArray[i]] = remoteProcIDs[i];
264  }
265 
266  // Remap aggregate GIDs to LIDs and set up owning processors
267  for (LO LID = 0; LID < numPoints; LID++) {
268  if (revMap_.find(getRoot(LID)) != revMap_.end() && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
269  GO aggGID = vertex2AggId[LID];
270  if (aggGID != MUELU_UNAGGREGATED) {
271  vertex2AggId[LID] = AggG2L[aggGID];
272  procWinner[LID] = AggG2R[aggGID];
273  }
274  }
275  }
276 
277  GO numGlobalRemote;
278  MueLu_sumAll(comm, as<GO>(numRemote), numGlobalRemote);
279  aggregates->AggregatesCrossProcessors(numGlobalRemote);
280 
281  Set(currentLevel, "Aggregates", aggregates);
282 
283  GetOStream(Statistics1) << aggregates->description() << std::endl;
284 }
285 
286 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
288  Setup(const RCP<const Teuchos::Comm<int> >& comm, const RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO> >& coords, const RCP<const Map>& /* map */) const {
289  nDim_ = coords->getNumVectors();
290 
291  x_ = coords->getData(0);
292  xMap_ = Construct1DMap(comm, x_);
293  nx_ = xMap_->size();
294 
295  ny_ = 1;
296  if (nDim_ > 1) {
297  y_ = coords->getData(1);
298  yMap_ = Construct1DMap(comm, y_);
299  ny_ = yMap_->size();
300  }
301 
302  nz_ = 1;
303  if (nDim_ > 2) {
304  z_ = coords->getData(2);
305  zMap_ = Construct1DMap(comm, z_);
306  nz_ = zMap_->size();
307  }
308 
309  for (size_t ind = 0; ind < coords->getLocalLength(); ind++) {
310  GO i = (*xMap_)[(coords->getData(0))[ind]], j = 0, k = 0;
311  if (nDim_ > 1)
312  j = (*yMap_)[(coords->getData(1))[ind]];
313  if (nDim_ > 2)
314  k = (*zMap_)[(coords->getData(2))[ind]];
315 
316  revMap_[k * ny_ * nx_ + j * nx_ + i] = ind;
317  }
318 
319  // Get the number of aggregates in each direction, correcting for Dirichlet
320  int xboost = dirichletX_ ? 1 : 0;
321  int yboost = dirichletY_ ? 1 : 0;
322  int zboost = dirichletZ_ ? 1 : 0;
323  naggx_ = (nx_ - 2 * xboost) / bx_ + ((nx_ - 2 * xboost) % bx_ ? 1 : 0);
324 
325  if (nDim_ > 1)
326  naggy_ = (ny_ - 2 * yboost) / by_ + ((ny_ - 2 * yboost) % by_ ? 1 : 0);
327  else
328  naggy_ = 1;
329 
330  if (nDim_ > 2)
331  naggz_ = (nz_ - 2 * zboost) / bz_ + ((nz_ - 2 * zboost) % bz_ ? 1 : 0);
332  else
333  naggz_ = 1;
334 }
335 
336 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
340  const ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& x) const {
341  int n = x.size();
342 
343  // Step 1: Create a local vector with unique coordinate points
344  RCP<container> gMap = rcp(new container);
345  for (int i = 0; i < n; i++)
346  (*gMap)[x[i]] = 0;
347 
348 #ifdef HAVE_MPI
349  // Step 2: exchange coordinates
350  // NOTE: we assume the coordinates are double, or double compatible
351  // That means that for complex case, we assume that all imaginary parts are zeros
352  int numProcs = comm->getSize();
353  if (numProcs > 1) {
354  RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm->duplicate());
355 
356  MPI_Comm rawComm = (*dupMpiComm->getRawMpiComm())();
357 
358  int sendCnt = gMap->size(), cnt = 0, recvSize;
359  Array<int> recvCnt(numProcs), Displs(numProcs);
360  Array<double> sendBuf, recvBuf;
361 
362  sendBuf.resize(sendCnt);
363  for (typename container::const_iterator cit = gMap->begin(); cit != gMap->end(); cit++)
364  sendBuf[cnt++] = Teuchos::as<double>(STS::real(cit->first));
365 
366  MPI_Allgather(&sendCnt, 1, MPI_INT, recvCnt.getRawPtr(), 1, MPI_INT, rawComm);
367  Displs[0] = 0;
368  for (int i = 0; i < numProcs - 1; i++)
369  Displs[i + 1] = Displs[i] + recvCnt[i];
370  recvSize = Displs[numProcs - 1] + recvCnt[numProcs - 1];
371  recvBuf.resize(recvSize);
372  MPI_Allgatherv(sendBuf.getRawPtr(), sendCnt, MPI_DOUBLE, recvBuf.getRawPtr(), recvCnt.getRawPtr(), Displs.getRawPtr(), MPI_DOUBLE, rawComm);
373 
374  for (int i = 0; i < recvSize; i++)
375  (*gMap)[as<SC>(recvBuf[i])] = 0;
376  }
377 #endif
378 
379  GO cnt = 0;
380  for (typename container::iterator it = gMap->begin(); it != gMap->end(); it++)
381  it->second = cnt++;
382 
383  return gMap;
384 }
385 
386 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
388  int i, j, k;
389  getIJK(LID, i, j, k);
390 
391  return (k * ny_ * nx_ + j * nx_ + i) == getRoot(LID);
392 }
393 
394 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
396  bool boundary = false;
397  int i, j, k;
398  getIJK(LID, i, j, k);
399  if (dirichletX_ && (i == 0 || i == nx_ - 1))
400  boundary = true;
401  if (nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_ - 1))
402  boundary = true;
403  if (nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_ - 1))
404  boundary = true;
405 
406  return boundary;
407 }
408 
409 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
411  if (isDirichlet(LID))
413 
414  int aggI, aggJ, aggK;
415  getAggIJK(LID, aggI, aggJ, aggK);
416  int xboost = dirichletX_ ? 1 : 0;
417  int yboost = dirichletY_ ? 1 : 0;
418  int zboost = dirichletZ_ ? 1 : 0;
419 
420  int i = xboost + aggI * bx_ + (bx_ - 1) / 2;
421  int j = (nDim_ > 1) ? yboost + aggJ * by_ + (by_ - 1) / 2 : 0;
422  int k = (nDim_ > 2) ? zboost + aggK * bz_ + (bz_ - 1) / 2 : 0;
423 
424  return k * ny_ * nx_ + j * nx_ + i;
425 }
426 
427 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
429  i = (*xMap_)[x_[LID]];
430  j = (nDim_ > 1) ? (*yMap_)[y_[LID]] : 0;
431  k = (nDim_ > 2) ? (*zMap_)[z_[LID]] : 0;
432 }
433 
434 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
436  int xboost = dirichletX_ ? 1 : 0;
437  int yboost = dirichletY_ ? 1 : 0;
438  int zboost = dirichletZ_ ? 1 : 0;
439  int pointI, pointJ, pointK;
440  getIJK(LID, pointI, pointJ, pointK);
441  i = (pointI - xboost) / bx_;
442 
443  if (nDim_ > 1)
444  j = (pointJ - yboost) / by_;
445  else
446  j = 0;
447 
448  if (nDim_ > 2)
449  k = (pointK - zboost) / bz_;
450  else
451  k = 0;
452 }
453 
454 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
456  bool boundary = false;
457 
458  int i, j, k;
459  getIJK(LID, i, j, k);
460  int ii, jj, kk;
461  getAggIJK(LID, ii, jj, kk);
462 
463  if (dirichletX_ && (i == 0 || i == nx_ - 1)) boundary = true;
464  if (nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_ - 1)) boundary = true;
465  if (nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_ - 1)) boundary = true;
466 
467  /*
468  if(boundary)
469  printf("[%d] coord = (%d,%d,%d) {%d,%d,%d} agg = (%d,%d,%d) {%d,%d,%d} => agg %s\n",LID,i,j,k,nx_,ny_,nz_,ii,jj,kk,naggx_,naggy_,naggz_,"BOUNDARY");
470  else
471  printf("[%d] coord = (%d,%d,%d) {%d,%d,%d} agg = (%d,%d,%d) {%d,%d,%d} => agg %d\n",LID,i,j,k,nx_,ny_,nz_,ii,jj,kk,naggx_,naggy_,naggz_,kk*naggy_*naggx_ + jj*naggx_ + ii);
472  */
473 
474  if (boundary)
476  else
477  return Teuchos::as<GlobalOrdinal>(kk * naggy_ * naggx_) + Teuchos::as<GlobalOrdinal>(jj * naggx_) + ii;
478 }
479 
480 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
482  // TODO: Currently only works w/ 1 DOF per node
483  double dirichletThreshold = 0.0;
484 
485  if (bx_ > 1 && (nDim_ <= 1 || by_ > 1) && (nDim_ <= 2 || bz_ > 1)) {
486  FactoryMonitor m(*this, "Generating Graph (trivial)", currentLevel);
487  /*** Case 1: Use the matrix is the graph ***/
488  // Bricks are of non-trivial size in all active dimensions
489  RCP<LWGraph> graph = rcp(new LWGraph(A->getCrsGraph(), "graph of A"));
490  auto boundaryNodes = MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows_kokkos_host(*A, dirichletThreshold);
491  graph->SetBoundaryNodeMap(boundaryNodes);
492 
493  if (GetVerbLevel() & Statistics1) {
494  GO numLocalBoundaryNodes = 0;
495  GO numGlobalBoundaryNodes = 0;
496  for (size_t i = 0; i < boundaryNodes.size(); ++i)
497  if (boundaryNodes(i))
498  numLocalBoundaryNodes++;
499  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
500  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
501  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
502  }
503  Set(currentLevel, "DofsPerNode", 1);
504  Set(currentLevel, "Graph", graph);
505  Set(currentLevel, "Filtering", false);
506  } else {
507  FactoryMonitor m(*this, "Generating Graph", currentLevel);
508  /*** Case 2: Dropping required ***/
509  // There is at least one active dimension in which we are not coarsening.
510  // Those connections need to be dropped
511  bool drop_x = (bx_ == 1);
512  bool drop_y = (nDim_ > 1 && by_ == 1);
513  bool drop_z = (nDim_ > 2 && bz_ == 1);
514 
515  typename LWGraph::row_type::non_const_type rows("rows", A->getLocalNumRows() + 1);
516  typename LWGraph::entries_type::non_const_type columns("columns", A->getLocalNumEntries());
517 
518  size_t N = A->getRowMap()->getLocalNumElements();
519 
520  // FIXME: Do this on the host because indexing functions are host functions
521  auto G = A->getLocalMatrixHost().graph;
522  auto rowptr = G.row_map;
523  auto colind = G.entries;
524 
525  int ct = 0;
526  rows(0) = 0;
527  for (size_t row = 0; row < N; row++) {
528  // NOTE: Assumes that the first part of the colmap is the rowmap
529  int ir, jr, kr;
530  LO row2 = A->getColMap()->getLocalElement(A->getRowMap()->getGlobalElement(row));
531  getIJK(row2, ir, jr, kr);
532 
533  for (size_t cidx = rowptr[row]; cidx < rowptr[row + 1]; cidx++) {
534  int ic, jc, kc;
535  LO col = colind[cidx];
536  getIJK(col, ic, jc, kc);
537 
538  if ((row2 != col) && ((drop_x && ir != ic) || (drop_y && jr != jc) || (drop_z && kr != kc))) {
539  // Drop it
540  // printf("[%4d] DROP row = (%d,%d,%d) col = (%d,%d,%d)\n",(int)row,ir,jr,kr,ic,jc,kc);
541  } else {
542  // Keep it
543  // printf("[%4d] KEEP row = (%d,%d,%d) col = (%d,%d,%d)\n",(int)row,ir,jr,kr,ic,jc,kc);
544  columns(ct) = col;
545  ct++;
546  }
547  }
548  rows(row + 1) = ct;
549  } // end for
550 
551  RCP<LWGraph> graph = rcp(new LWGraph(rows, columns, A->getRowMap(), A->getColMap(), "thresholded graph of A"));
552 
553  auto boundaryNodes = MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows_kokkos_host(*A, dirichletThreshold);
554  graph->SetBoundaryNodeMap(boundaryNodes);
555 
556  if (GetVerbLevel() & Statistics1) {
557  GO numLocalBoundaryNodes = 0;
558  GO numGlobalBoundaryNodes = 0;
559  for (size_t i = 0; i < boundaryNodes.size(); ++i)
560  if (boundaryNodes(i))
561  numLocalBoundaryNodes++;
562  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
563  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
564  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
565  }
566  Set(currentLevel, "DofsPerNode", 1);
567  Set(currentLevel, "Graph", graph);
568  Set(currentLevel, "Filtering", true);
569  } // end else
570 
571 } // end BuildGraph
572 
573 } // namespace MueLu
574 
575 #endif /* MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_ */
#define MueLu_sumAll(rcpComm, in, out)
void getAggIJK(LocalOrdinal LID, int &i, int &j, int &k) const
MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
const RCP< LOVector > & GetProcWinner() const
Returns constant vector that maps local node IDs to owning processor IDs.
KOKKOS_INLINE_FUNCTION void SetBoundaryNodeMap(const boundary_nodes_type bndry)
Set boolean array indicating which rows correspond to Dirichlet boundaries.
Container class for aggregation information.
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
LocalOrdinal LO
One-liner description of what is happening.
void Build(Level &currentLevel) const
Build aggregates.
void SetIsRoot(LO i, bool value=true)
Set root node information.
GlobalOrdinal getAggGID(LocalOrdinal LID) const
std::map< Scalar, GlobalOrdinal, compare > container
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Setup(const RCP< const Teuchos::Comm< int > > &comm, const RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LO, GO, NO > > &coords, const RCP< const Map > &map) const
#define MUELU_UNAGGREGATED
virtual void setObjectLabel(const std::string &objectLabel)
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
void resize(size_type new_size, const value_type &x=value_type())
void getIJK(LocalOrdinal LID, int &i, int &j, int &k) const
Lightweight MueLu representation of a compressed row storage graph.
RCP< container > Construct1DMap(const RCP< const Teuchos::Comm< int > > &comm, const ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &x) const
KOKKOS_INLINE_FUNCTION void AggregatesCrossProcessors(const bool &flag)
Record whether aggregates include DOFs from other processes.
Node NO
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)
std::string description() const
Return a simple one-line description of this object.
void BuildGraph(Level &currentLevel, const RCP< Matrix > &A) const
GlobalOrdinal getRoot(LocalOrdinal LID) const
void DeclareInput(Level &currentLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.
#define SET_VALID_ENTRY(name)
bool is_null() const