46#ifndef MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_
47#define MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_
51#include <Teuchos_DefaultMpiComm.hpp>
52#include <Teuchos_CommHelpers.hpp>
54#include <Teuchos_OrdinalTraits.hpp>
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>
64#include "MueLu_Aggregates.hpp"
68#include "MueLu_Utilities.hpp"
70#include "MueLu_Graph.hpp"
71#include "MueLu_LWGraph.hpp"
76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 RCP<ParameterList> validParamList = rcp(
new ParameterList());
80#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
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;
94 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 Input(currentLevel,
"A");
97 Input(currentLevel,
"Coordinates");
144 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
148 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> MultiVector_d;
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();
157 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
158 int numProcs = comm->getSize();
159 int myRank = comm->getRank();
161 int numPoints = colMap->getLocalNumElements();
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");
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;
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);
188 Setup(comm, overlappedCoords, colMap);
190 GetOStream(
Runtime0) <<
"Using brick size: " << bx_
191 << (nDim_ > 1 ?
"x " +
toString(by_) :
"")
192 << (nDim_ > 2 ?
"x " +
toString(bz_) :
"") << std::endl;
195 BuildGraph(currentLevel,A);
198 RCP<Aggregates> aggregates = rcp(
new Aggregates(colMap));
199 aggregates->setObjectLabel(
"Brick");
201 ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
202 ArrayRCP<LO> procWinner = aggregates->GetProcWinner() ->getDataNonConst(0);
211 std::set<GO> myAggGIDs, remoteAggGIDs;
212 for (LO LID = 0; LID < numPoints; LID++) {
213 GO aggGID = getAggGID(LID);
215 if(aggGID == GO_INVALID)
continue;
218 if ((revMap_.find(getRoot(LID)) != revMap_.end()) && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
220 vertex2AggId[LID] = aggGID;
221 myAggGIDs.insert(aggGID);
224 aggregates->SetIsRoot(LID);
227 remoteAggGIDs.insert(aggGID);
230 size_t numAggregates = myAggGIDs .size();
231 size_t numRemote = remoteAggGIDs.size();
232 aggregates->SetNumAggregates(numAggregates);
234 std::map<GO,LO> AggG2L;
235 std::map<GO,int> AggG2R;
237 Array<GO> myAggGIDsArray(numAggregates), remoteAggGIDsArray(numRemote);
241 for (
typename std::set<GO>::const_iterator it = myAggGIDs.begin(); it != myAggGIDs.end(); it++) {
243 AggG2R[*it] = myRank;
245 myAggGIDsArray[ind++] = *it;
249 RCP<Map> aggMap = MapFactory::Build(rowMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
250 myAggGIDsArray, 0, comm);
253 for (
typename std::set<GO>::const_iterator it = remoteAggGIDs.begin(); it != remoteAggGIDs.end(); it++)
254 remoteAggGIDsArray[ind++] = *it;
257 Array<int> remoteProcIDs(numRemote);
258 Array<LO> remoteLIDs (numRemote);
259 aggMap->getRemoteIndexList(remoteAggGIDsArray, remoteProcIDs, remoteLIDs);
262 for (
size_t i = 0; i < numRemote; i++) {
263 AggG2L[remoteAggGIDsArray[i]] = remoteLIDs [i];
264 AggG2R[remoteAggGIDsArray[i]] = remoteProcIDs[i];
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];
272 vertex2AggId[LID] = AggG2L[aggGID];
273 procWinner [LID] = AggG2R[aggGID];
281 aggregates->AggregatesCrossProcessors(numGlobalRemote);
283 Set(currentLevel,
"Aggregates", aggregates);
285 GetOStream(
Statistics1) << aggregates->description() << std::endl;
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>& )
const {
291 nDim_ = coords->getNumVectors();
293 x_ = coords->getData(0);
294 xMap_ = Construct1DMap(comm, x_);
299 y_ = coords->getData(1);
300 yMap_ = Construct1DMap(comm, y_);
306 z_ = coords->getData(2);
307 zMap_ = Construct1DMap(comm, z_);
311 for (
size_t ind = 0; ind < coords->getLocalLength(); ind++) {
312 GO i = (*xMap_)[(coords->getData(0))[ind]], j = 0, k = 0;
314 j = (*yMap_)[(coords->getData(1))[ind]];
316 k = (*zMap_)[(coords->getData(2))[ind]];
318 revMap_[k*ny_*nx_ + j*nx_ + i] = ind;
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);
329 naggy_ = (ny_-2*yboost)/by_ + ( (ny_-2*yboost) % by_ ? 1 : 0);
334 naggz_ = (nz_-2*zboost)/bz_ + ( (nz_-2*zboost) % bz_ ? 1 : 0);
340 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
341 RCP<typename BrickAggregationFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::container>
344 const ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& x)
const
349 RCP<container> gMap = rcp(
new container);
350 for (
int i = 0; i < n; i++)
357 int numProcs = comm->getSize();
359 RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm->duplicate());
361 MPI_Comm rawComm = (*dupMpiComm->getRawMpiComm())();
363 int sendCnt = gMap->size(), cnt = 0, recvSize;
364 Array<int> recvCnt(numProcs), Displs(numProcs);
365 Array<double> sendBuf, recvBuf;
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));
371 MPI_Allgather(&sendCnt, 1, MPI_INT, recvCnt.getRawPtr(), 1, MPI_INT, rawComm);
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);
379 for (
int i = 0; i < recvSize; i++)
380 (*gMap)[as<SC>(recvBuf[i])] = 0;
385 for (
typename container::iterator it = gMap->begin(); it != gMap->end(); it++)
391 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
396 return (k*ny_*nx_ + j*nx_ + i) == getRoot(LID);
399 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
401 bool boundary =
false;
404 if( dirichletX_ && (i == 0 || i == nx_-1) )
406 if(nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_-1) )
408 if(nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_-1) )
415 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
418 return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
421 getAggIJK(LID,aggI,aggJ,aggK);
422 int xboost = dirichletX_ ? 1 : 0;
423 int yboost = dirichletY_ ? 1 : 0;
424 int zboost = dirichletZ_ ? 1 : 0;
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;
430 return k*ny_*nx_ + j*nx_ + i;
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;
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_;
450 if (nDim_ > 1) j = (pointJ-yboost)/by_;
453 if (nDim_ > 2) k = (pointK-zboost)/bz_;
457 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
459 bool boundary =
false;
464 getAggIJK(LID,ii,jj,kk);
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;
478 return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
480 return Teuchos::as<GlobalOrdinal>(kk*naggy_*naggx_) + Teuchos::as<GlobalOrdinal>(jj*naggx_) + ii;
485 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
488 double dirichletThreshold = 0.0;
490 if(bx_ > 1 && (nDim_ <= 1 || by_ > 1) && (nDim_ <=2 || bz_>1) ) {
491 FactoryMonitor m(*
this,
"Generating Graph (trivial)", currentLevel);
494 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
496 graph->SetBoundaryNodeMap(boundaryNodes);
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;
508 Set(currentLevel,
"DofsPerNode", 1);
509 Set(currentLevel,
"Graph", graph);
510 Set(currentLevel,
"Filtering",
false);
517 bool drop_x = (bx_ == 1);
518 bool drop_y = (nDim_> 1 && by_ == 1);
519 bool drop_z = (nDim_> 2 && bz_ == 1);
521 ArrayRCP<LO>
rows (A->getLocalNumRows()+1);
522 ArrayRCP<LO> columns(A->getLocalNumEntries());
524 size_t N = A->getRowMap()->getLocalNumElements();
527 auto G = A->getLocalMatrixHost().graph;
528 auto rowptr = G.row_map;
529 auto colind = G.entries;
533 for(
size_t row=0; row<N; row++) {
536 LO row2 = A->getColMap()->getLocalElement(A->getRowMap()->getGlobalElement(row));
537 getIJK(row2,ir,jr,kr);
539 for(
size_t cidx=rowptr[row]; cidx<rowptr[row+1]; cidx++) {
541 LO col = colind[cidx];
542 getIJK(col,ic,jc,kc);
544 if( (row2 !=col) && ((drop_x && ir != ic) || (drop_y && jr != jc) || (drop_z && kr != kc) )) {
558 RCP<GraphBase> graph = rcp(
new LWGraph(
rows, columns, A->getRowMap(), A->getColMap(),
"thresholded graph of A"));
562 graph->SetBoundaryNodeMap(boundaryNodes);
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;
574 Set(currentLevel,
"DofsPerNode", 1);
575 Set(currentLevel,
"Graph", graph);
576 Set(currentLevel,
"Filtering",
true);
#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
bool isDirichlet(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 ¤tLevel) const
Input.
void Build(Level ¤tLevel) const
Build aggregates.
void BuildGraph(Level ¤tLevel, 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
bool isRoot(LocalOrdinal LID) 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.
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.