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;
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);
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++)
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);