MueLu Version of the Day
Loading...
Searching...
No Matches
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
51#include <Teuchos_DefaultMpiComm.hpp>
52#include <Teuchos_CommHelpers.hpp>
53#endif
54#include <Teuchos_OrdinalTraits.hpp>
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_GraphBase.hpp"
70#include "MueLu_Graph.hpp"
71#include "MueLu_LWGraph.hpp"
72
73
74namespace MueLu {
75
76 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
78 RCP<ParameterList> validParamList = rcp(new ParameterList());
79
80#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
81 SET_VALID_ENTRY("aggregation: brick x size");
82 SET_VALID_ENTRY("aggregation: brick y size");
83 SET_VALID_ENTRY("aggregation: brick z size");
84 SET_VALID_ENTRY("aggregation: brick x Dirichlet");
85 SET_VALID_ENTRY("aggregation: brick y Dirichlet");
86 SET_VALID_ENTRY("aggregation: brick z Dirichlet");
87#undef SET_VALID_ENTRY
88
89 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory for matrix");
90 validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for coordinates");
91 return validParamList;
92 }
93
94 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
96 Input(currentLevel, "A");
97 Input(currentLevel, "Coordinates");
98 }
99
100 // The current implementation cannot deal with bricks larger than 3x3(x3) in
101 // parallel. The reason is that aggregation infrastructure in place has
102 // major drawbacks.
103 //
104 // Aggregates class is constructed with a help of a provided map, either
105 // taken from a graph, or provided directly. This map is usually taken to be
106 // a column map of a matrix. The reason for that is that if we have an
107 // overlapped aggregation, we want the processor owning aggregates to store
108 // agg id for all nodes in this aggregate. If we used row map, there would
109 // be no way for the processor to know whether there are some other nodes on
110 // a different processor which belong to its aggregate. On the other hand,
111 // using column map allows both vertex2AggId and procWinner arrays in
112 // Aggregates class to store some extra data, such as whether nodes belonging
113 // to a different processor belong to this processor aggregate.
114 //
115 // The drawback of this is that it stores only overlap=1 data. For aggressive
116 // coarsening, such a brick aggregation with a large single dimension of
117 // brick, it could happen that we need to know depth two or more extra nodes
118 // in the other processor subdomain.
119 //
120 // Another issue is that we may have some implicit connection between
121 // aggregate map and maps of A used in the construction of a tentative
122 // prolongator.
123 //
124 // Another issue is that it seems that some info is unused or not required.
125 // Specifically, it seems that if a node belongs to an aggregate on a
126 // different processor, we don't actually need to set vertex2AggId and
127 // procWinner, despite the following comment in
128 // Aggregates decl:
129 // vertex2AggId[k] gives a local id
130 // corresponding to the aggregate to which
131 // local id k has been assigned. While k
132 // is the local id on my processor (MyPID)
133 // vertex2AggId[k] is the local id on the
134 // processor which actually owns the
135 // aggregate. This owning processor has id
136 // given by procWinner[k].
137 // It is possible that that info is only used during arbitration in
138 // CoupledAggregationFactory.
139 //
140 // The steps that we need to do to resolve this issue:
141 // - Break the link between maps in TentativePFactory, allowing any maps in Aggregates
142 // - Allow Aggregates to construct their own maps, if necessary, OR
143 // - construct aggregates based on row map
144 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
146 FactoryMonitor m(*this, "Build", currentLevel);
147
148 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> MultiVector_d;
149
150 const ParameterList& pL = GetParameterList();
151 RCP<MultiVector_d> coords = Get<RCP<MultiVector_d> >(currentLevel, "Coordinates");
152 RCP<Matrix> A = Get< RCP<Matrix> > (currentLevel, "A");
153 RCP<const Map> rowMap = A->getRowMap();
154 RCP<const Map> colMap = A->getColMap();
155 GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
156
157 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
158 int numProcs = comm->getSize();
159 int myRank = comm->getRank();
160
161 int numPoints = colMap->getLocalNumElements();
162
163 bx_ = pL.get<int>("aggregation: brick x size");
164 by_ = pL.get<int>("aggregation: brick y size");
165 bz_ = pL.get<int>("aggregation: brick z size");
166
167 dirichletX_ = pL.get<bool>("aggregation: brick x Dirichlet");
168 dirichletY_ = pL.get<bool>("aggregation: brick y Dirichlet");
169 dirichletZ_ = pL.get<bool>("aggregation: brick z Dirichlet");
170 if(dirichletX_) GetOStream(Runtime0) << "Dirichlet boundaries in the x direction"<<std::endl;
171 if(dirichletY_) GetOStream(Runtime0) << "Dirichlet boundaries in the y direction"<<std::endl;
172 if(dirichletZ_) GetOStream(Runtime0) << "Dirichlet boundaries in the z direction"<<std::endl;
173
174 if (numProcs > 1) {
175 // TODO: deal with block size > 1 (see comments above)
176 //TEUCHOS_TEST_FOR_EXCEPTION(bx_ > 3 || by_ > 3 || bz_ > 3, Exceptions::RuntimeError, "Currently cannot deal with brick size > 3");
177 }
178
179 RCP<MultiVector_d> overlappedCoords = coords;
180 RCP<const Import> importer = ImportFactory::Build(coords->getMap(), colMap);
181 if (!importer.is_null()) {
182 overlappedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(colMap, coords->getNumVectors());
183 overlappedCoords->doImport(*coords, *importer, Xpetra::INSERT);
184 }
185
186 // Setup misc structures
187 // Logically, we construct enough data to query topological information of a rectangular grid
188 Setup(comm, overlappedCoords, colMap);
189
190 GetOStream(Runtime0) << "Using brick size: " << bx_
191 << (nDim_ > 1 ? "x " + toString(by_) : "")
192 << (nDim_ > 2 ? "x " + toString(bz_) : "") << std::endl;
193
194 // Build the graph
195 BuildGraph(currentLevel,A);
196
197 // Construct aggregates
198 RCP<Aggregates> aggregates = rcp(new Aggregates(colMap));
199 aggregates->setObjectLabel("Brick");
200
201 ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
202 ArrayRCP<LO> procWinner = aggregates->GetProcWinner() ->getDataNonConst(0);
203
204 // In the first pass, we set a mapping from a vertex to aggregate global id. We deal with a structured
205 // rectangular mesh, therefore we know the structure of aggregates. For each vertex we can tell exactly
206 // which aggregate it belongs to.
207 // If we determine that the aggregate does not belong to us (i.e. the root vertex does not belong to this
208 // processor, or is outside and we lost "" arbitration), we record the global aggregate id in order to
209 // fetch the local info from the processor owning the aggregate. This is required for aggregates, as it
210 // uses the local aggregate ids of the owning processor.
211 std::set<GO> myAggGIDs, remoteAggGIDs;
212 for (LO LID = 0; LID < numPoints; LID++) {
213 GO aggGID = getAggGID(LID);
214 // 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);
215 if(aggGID == GO_INVALID) continue;
216 // printf("[%d] getRoot = %d\n",(int)LID,(int)getRoot(LID));
217
218 if ((revMap_.find(getRoot(LID)) != revMap_.end()) && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
219 // Root of the brick aggregate containing GID (<- LID) belongs to us
220 vertex2AggId[LID] = aggGID;
221 myAggGIDs.insert(aggGID);
222
223 if (isRoot(LID))
224 aggregates->SetIsRoot(LID);
225 // printf("[%d] initial vertex2AggId = %d\n",(int)LID,(int)vertex2AggId[LID]);
226 } else {
227 remoteAggGIDs.insert(aggGID);
228 }
229 }
230 size_t numAggregates = myAggGIDs .size();
231 size_t numRemote = remoteAggGIDs.size();
232 aggregates->SetNumAggregates(numAggregates);
233
234 std::map<GO,LO> AggG2L; // Map: Agg GID -> Agg LID (possibly on a different processor)
235 std::map<GO,int> AggG2R; // Map: Agg GID -> processor rank owning aggregate
236
237 Array<GO> myAggGIDsArray(numAggregates), remoteAggGIDsArray(numRemote);
238
239 // Fill in the maps for aggregates that we own
240 size_t ind = 0;
241 for (typename std::set<GO>::const_iterator it = myAggGIDs.begin(); it != myAggGIDs.end(); it++) {
242 AggG2L[*it] = ind;
243 AggG2R[*it] = myRank;
244
245 myAggGIDsArray[ind++] = *it;
246 }
247
248 // The map is a convenient way to fetch remote local indices from global indices.
249 RCP<Map> aggMap = MapFactory::Build(rowMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
250 myAggGIDsArray, 0, comm);
251
252 ind = 0;
253 for (typename std::set<GO>::const_iterator it = remoteAggGIDs.begin(); it != remoteAggGIDs.end(); it++)
254 remoteAggGIDsArray[ind++] = *it;
255
256 // Fetch the required aggregate local ids and ranks
257 Array<int> remoteProcIDs(numRemote);
258 Array<LO> remoteLIDs (numRemote);
259 aggMap->getRemoteIndexList(remoteAggGIDsArray, remoteProcIDs, remoteLIDs);
260
261 // Fill in the maps for aggregates that we don't own but which have some of our vertices
262 for (size_t i = 0; i < numRemote; i++) {
263 AggG2L[remoteAggGIDsArray[i]] = remoteLIDs [i];
264 AggG2R[remoteAggGIDsArray[i]] = remoteProcIDs[i];
265 }
266
267 // Remap aggregate GIDs to LIDs and set up owning processors
268 for (LO LID = 0; LID < numPoints; LID++) {
269 if (revMap_.find(getRoot(LID)) != revMap_.end() && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
270 GO aggGID = vertex2AggId[LID];
271 if(aggGID != MUELU_UNAGGREGATED) {
272 vertex2AggId[LID] = AggG2L[aggGID];
273 procWinner [LID] = AggG2R[aggGID];
274 }
275 }
276 }
277
278
279 GO numGlobalRemote;
280 MueLu_sumAll(comm, as<GO>(numRemote), numGlobalRemote);
281 aggregates->AggregatesCrossProcessors(numGlobalRemote);
282
283 Set(currentLevel, "Aggregates", aggregates);
284
285 GetOStream(Statistics1) << aggregates->description() << std::endl;
286 }
287
288 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
290 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 {
291 nDim_ = coords->getNumVectors();
292
293 x_ = coords->getData(0);
294 xMap_ = Construct1DMap(comm, x_);
295 nx_ = xMap_->size();
296
297 ny_ = 1;
298 if (nDim_ > 1) {
299 y_ = coords->getData(1);
300 yMap_ = Construct1DMap(comm, y_);
301 ny_ = yMap_->size();
302 }
303
304 nz_ = 1;
305 if (nDim_ > 2) {
306 z_ = coords->getData(2);
307 zMap_ = Construct1DMap(comm, z_);
308 nz_ = zMap_->size();
309 }
310
311 for (size_t ind = 0; ind < coords->getLocalLength(); ind++) {
312 GO i = (*xMap_)[(coords->getData(0))[ind]], j = 0, k = 0;
313 if (nDim_ > 1)
314 j = (*yMap_)[(coords->getData(1))[ind]];
315 if (nDim_ > 2)
316 k = (*zMap_)[(coords->getData(2))[ind]];
317
318 revMap_[k*ny_*nx_ + j*nx_ + i] = ind;
319 }
320
321
322 // Get the number of aggregates in each direction, correcting for Dirichlet
323 int xboost = dirichletX_ ? 1 : 0;
324 int yboost = dirichletY_ ? 1 : 0;
325 int zboost = dirichletZ_ ? 1 : 0;
326 naggx_ = (nx_-2*xboost)/bx_ + ((nx_-2*xboost) % bx_ ? 1 : 0);
327
328 if(nDim_ > 1)
329 naggy_ = (ny_-2*yboost)/by_ + ( (ny_-2*yboost) % by_ ? 1 : 0);
330 else
331 naggy_ = 1;
332
333 if(nDim_ > 2)
334 naggz_ = (nz_-2*zboost)/bz_ + ( (nz_-2*zboost) % bz_ ? 1 : 0);
335 else
336 naggz_ = 1;
337
338 }
339
340 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
341 RCP<typename BrickAggregationFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::container>
343 Construct1DMap (const RCP<const Teuchos::Comm<int> >& comm,
344 const ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& x) const
345 {
346 int n = x.size();
347
348 // Step 1: Create a local vector with unique coordinate points
349 RCP<container> gMap = rcp(new container);
350 for (int i = 0; i < n; i++)
351 (*gMap)[x[i]] = 0;
352
353#ifdef HAVE_MPI
354 // Step 2: exchange coordinates
355 // NOTE: we assume the coordinates are double, or double compatible
356 // That means that for complex case, we assume that all imaginary parts are zeros
357 int numProcs = comm->getSize();
358 if (numProcs > 1) {
359 RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm->duplicate());
360
361 MPI_Comm rawComm = (*dupMpiComm->getRawMpiComm())();
362
363 int sendCnt = gMap->size(), cnt = 0, recvSize;
364 Array<int> recvCnt(numProcs), Displs(numProcs);
365 Array<double> sendBuf, recvBuf;
366
367 sendBuf.resize(sendCnt);
368 for (typename container::const_iterator cit = gMap->begin(); cit != gMap->end(); cit++)
369 sendBuf[cnt++] = Teuchos::as<double>(STS::real(cit->first));
370
371 MPI_Allgather(&sendCnt, 1, MPI_INT, recvCnt.getRawPtr(), 1, MPI_INT, rawComm);
372 Displs[0] = 0;
373 for (int i = 0; i < numProcs-1; i++)
374 Displs[i+1] = Displs[i] + recvCnt[i];
375 recvSize = Displs[numProcs-1] + recvCnt[numProcs-1];
376 recvBuf.resize(recvSize);
377 MPI_Allgatherv(sendBuf.getRawPtr(), sendCnt, MPI_DOUBLE, recvBuf.getRawPtr(), recvCnt.getRawPtr(), Displs.getRawPtr(), MPI_DOUBLE, rawComm);
378
379 for (int i = 0; i < recvSize; i++)
380 (*gMap)[as<SC>(recvBuf[i])] = 0;
381 }
382#endif
383
384 GO cnt = 0;
385 for (typename container::iterator it = gMap->begin(); it != gMap->end(); it++)
386 it->second = cnt++;
387
388 return gMap;
389 }
390
391 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
393 int i,j,k;
394 getIJK(LID,i,j,k);
395
396 return (k*ny_*nx_ + j*nx_ + i) == getRoot(LID);
397 }
398
399 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
401 bool boundary = false;
402 int i,j,k;
403 getIJK(LID,i,j,k);
404 if( dirichletX_ && (i == 0 || i == nx_-1) )
405 boundary = true;
406 if(nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_-1) )
407 boundary = true;
408 if(nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_-1) )
409 boundary = true;
410
411 return boundary;
412 }
413
414
415 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
417 if(isDirichlet(LID))
418 return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
419
420 int aggI,aggJ,aggK;
421 getAggIJK(LID,aggI,aggJ,aggK);
422 int xboost = dirichletX_ ? 1 : 0;
423 int yboost = dirichletY_ ? 1 : 0;
424 int zboost = dirichletZ_ ? 1 : 0;
425
426 int i = xboost + aggI*bx_ + (bx_-1)/2;
427 int j = (nDim_>1) ? yboost + aggJ*by_ + (by_-1)/2 : 0;
428 int k = (nDim_>2) ? zboost + aggK*bz_ + (bz_-1)/2 : 0;
429
430 return k*ny_*nx_ + j*nx_ + i;
431 }
432
433 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
435 i = (*xMap_)[x_[LID]];
436 j = (nDim_>1) ? (*yMap_)[y_[LID]] : 0;
437 k = (nDim_>2) ? (*zMap_)[z_[LID]] : 0;
438 }
439
440
441 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
443 int xboost = dirichletX_ ? 1 : 0;
444 int yboost = dirichletY_ ? 1 : 0;
445 int zboost = dirichletZ_ ? 1 : 0;
446 int pointI, pointJ, pointK;
447 getIJK(LID,pointI,pointJ,pointK);
448 i = (pointI-xboost)/bx_;
449
450 if (nDim_ > 1) j = (pointJ-yboost)/by_;
451 else j = 0;
452
453 if (nDim_ > 2) k = (pointK-zboost)/bz_;
454 else k = 0;
455 }
456
457 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
459 bool boundary = false;
460
461 int i, j, k;
462 getIJK(LID,i,j,k);
463 int ii , jj, kk;
464 getAggIJK(LID,ii,jj,kk);
465
466 if( dirichletX_ && (i == 0 || i == nx_ - 1)) boundary = true;
467 if (nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_ - 1)) boundary = true;
468 if (nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_ - 1)) boundary = true;
469
470 /*
471 if(boundary)
472 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");
473 else
474 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);
475 */
476
477 if (boundary)
478 return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
479 else
480 return Teuchos::as<GlobalOrdinal>(kk*naggy_*naggx_) + Teuchos::as<GlobalOrdinal>(jj*naggx_) + ii;
481
482 }
483
484
485 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
487 // TODO: Currently only works w/ 1 DOF per node
488 double dirichletThreshold = 0.0;
489
490 if(bx_ > 1 && (nDim_ <= 1 || by_ > 1) && (nDim_ <=2 || bz_>1) ) {
491 FactoryMonitor m(*this, "Generating Graph (trivial)", currentLevel);
492 /*** Case 1: Use the matrix is the graph ***/
493 // Bricks are of non-trivial size in all active dimensions
494 RCP<GraphBase> graph = rcp(new Graph(A->getCrsGraph(), "graph of A"));
495 ArrayRCP<bool > boundaryNodes = Teuchos::arcp_const_cast<bool>(MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold));
496 graph->SetBoundaryNodeMap(boundaryNodes);
497
498 if (GetVerbLevel() & Statistics1) {
499 GO numLocalBoundaryNodes = 0;
500 GO numGlobalBoundaryNodes = 0;
501 for (LO i = 0; i < boundaryNodes.size(); ++i)
502 if (boundaryNodes[i])
503 numLocalBoundaryNodes++;
504 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
505 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
506 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
507 }
508 Set(currentLevel, "DofsPerNode", 1);
509 Set(currentLevel, "Graph", graph);
510 Set(currentLevel, "Filtering",false);
511 }
512 else {
513 FactoryMonitor m(*this, "Generating Graph", currentLevel);
514 /*** Case 2: Dropping required ***/
515 // There is at least one active dimension in which we are not coarsening.
516 // Those connections need to be dropped
517 bool drop_x = (bx_ == 1);
518 bool drop_y = (nDim_> 1 && by_ == 1);
519 bool drop_z = (nDim_> 2 && bz_ == 1);
520
521 ArrayRCP<LO> rows (A->getLocalNumRows()+1);
522 ArrayRCP<LO> columns(A->getLocalNumEntries());
523
524 size_t N = A->getRowMap()->getLocalNumElements();
525
526 // FIXME: Do this on the host because indexing functions are host functions
527 auto G = A->getLocalMatrixHost().graph;
528 auto rowptr = G.row_map;
529 auto colind = G.entries;
530
531 int ct=0;
532 rows[0] = 0;
533 for(size_t row=0; row<N; row++) {
534 // NOTE: Assumes that the first part of the colmap is the rowmap
535 int ir,jr,kr;
536 LO row2 = A->getColMap()->getLocalElement(A->getRowMap()->getGlobalElement(row));
537 getIJK(row2,ir,jr,kr);
538
539 for(size_t cidx=rowptr[row]; cidx<rowptr[row+1]; cidx++) {
540 int ic,jc,kc;
541 LO col = colind[cidx];
542 getIJK(col,ic,jc,kc);
543
544 if( (row2 !=col) && ((drop_x && ir != ic) || (drop_y && jr != jc) || (drop_z && kr != kc) )) {
545 // Drop it
546 // printf("[%4d] DROP row = (%d,%d,%d) col = (%d,%d,%d)\n",(int)row,ir,jr,kr,ic,jc,kc);
547 }
548 else {
549 // Keep it
550 // printf("[%4d] KEEP row = (%d,%d,%d) col = (%d,%d,%d)\n",(int)row,ir,jr,kr,ic,jc,kc);
551 columns[ct] = col;
552 ct++;
553 }
554 }
555 rows[row+1] = ct;
556 }//end for
557
558 RCP<GraphBase> graph = rcp(new LWGraph(rows, columns, A->getRowMap(), A->getColMap(), "thresholded graph of A"));
559
560
561 ArrayRCP<bool > boundaryNodes = Teuchos::arcp_const_cast<bool>(MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold));
562 graph->SetBoundaryNodeMap(boundaryNodes);
563
564 if (GetVerbLevel() & Statistics1) {
565 GO numLocalBoundaryNodes = 0;
566 GO numGlobalBoundaryNodes = 0;
567 for (LO i = 0; i < boundaryNodes.size(); ++i)
568 if (boundaryNodes[i])
569 numLocalBoundaryNodes++;
570 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
571 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
572 GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
573 }
574 Set(currentLevel, "DofsPerNode", 1);
575 Set(currentLevel, "Graph", graph);
576 Set(currentLevel, "Filtering",true);
577 }//end else
578
579
580 }//end BuildGraph
581
582
583
584
585} //namespace MueLu
586
587#endif /* MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_ */
#define MUELU_UNAGGREGATED
#define SET_VALID_ENTRY(name)
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Container class for aggregation information.
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
GlobalOrdinal getRoot(LocalOrdinal LID) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
GlobalOrdinal getAggGID(LocalOrdinal LID) const
void DeclareInput(Level &currentLevel) const
Input.
void Build(Level &currentLevel) const
Build aggregates.
void BuildGraph(Level &currentLevel, const RCP< Matrix > &A) const
void getAggIJK(LocalOrdinal LID, int &i, int &j, int &k) const
RCP< container > Construct1DMap(const RCP< const Teuchos::Comm< int > > &comm, const ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &x) const
void getIJK(LocalOrdinal LID, int &i, int &j, int &k) const
std::map< Scalar, GlobalOrdinal, compare > container
Timer to be used in factories. Similar to Monitor but with additional timers.
MueLu representation of a compressed row storage graph.
Lightweight MueLu representation of a compressed row storage graph.
Class that holds all level-specific information.
MueLu utility class.
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.
@ Runtime0
One-liner description of what is happening.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.