25 #ifndef EWOMS_OVERLAPPING_BCRS_MATRIX_HH
26 #define EWOMS_OVERLAPPING_BCRS_MATRIX_HH
33 #include <opm/material/common/Valgrind.hpp>
35 #include <dune/istl/scalarproducts.hh>
36 #include <dune/istl/io.hh>
51 template <
class BCRSMatrix>
54 typedef BCRSMatrix ParentType;
60 typedef std::vector<std::set<Index> > Entries;
77 overlap_ = std::shared_ptr<Overlap>(
78 new Overlap(nativeMatrix, borderList, blackList, overlapSize));
81 MPI_Comm_rank(MPI_COMM_WORLD, &myRank_);
91 if (overlap_.use_count() == 0)
95 const PeerSet &peerSet = overlap_->peerSet();
96 typename PeerSet::const_iterator peerIt = peerSet.begin();
97 typename PeerSet::const_iterator peerEndIt = peerSet.end();
98 for (; peerIt != peerEndIt; ++peerIt) {
99 int peerRank = *peerIt;
101 delete rowSizesRecvBuff_[peerRank];
102 delete rowIndicesRecvBuff_[peerRank];
103 delete entryColIndicesRecvBuff_[peerRank];
104 delete entryValuesRecvBuff_[peerRank];
106 delete numRowsSendBuff_[peerRank];
107 delete rowSizesSendBuff_[peerRank];
108 delete rowIndicesSendBuff_[peerRank];
109 delete entryColIndicesSendBuff_[peerRank];
110 delete entryValuesSendBuff_[peerRank];
118 {
return *overlap_; }
126 assignFromNative_(nativeMatrix);
141 assignFromNative_(nativeMatrix);
153 block_type idMatrix(0.0);
154 for (
unsigned i = 0; i < idMatrix.size(); ++i)
155 idMatrix[i][i] = 1.0;
157 int numLocal = overlap_->numLocal();
158 int numDomestic = overlap_->numDomestic();
159 for (
int domRowIdx = numLocal; domRowIdx < numDomestic; ++domRowIdx) {
160 if (overlap_->isFront(domRowIdx)) {
162 (*this)[domRowIdx] = 0.0;
163 (*this)[domRowIdx][domRowIdx] = idMatrix;
172 for (
int i = 0; i < this->N(); ++i) {
173 if (overlap_->isLocal(i))
177 std::cout <<
"row " << i <<
" ";
179 typedef typename BCRSMatrix::ConstColIterator ColIt;
180 ColIt colIt = (*this)[i].begin();
181 ColIt colEndIt = (*this)[i].end();
182 for (
int j = 0; j < this->M(); ++j) {
183 if (colIt != colEndIt && j == colIt.index()) {
185 if (overlap_->isBorder(j))
187 else if (overlap_->isLocal(j))
195 std::cout <<
"\n" << std::flush;
197 Dune::printSparseMatrix(std::cout,
198 *static_cast<const BCRSMatrix *>(
this),
204 void assignFromNative_(
const BCRSMatrix &nativeMatrix)
207 BCRSMatrix::operator=(0.0);
210 for (
unsigned nativeRowIdx = 0; nativeRowIdx < nativeMatrix.N(); ++nativeRowIdx) {
211 int domesticRowIdx = overlap_->nativeToDomestic(nativeRowIdx);
212 if (domesticRowIdx < 0) {
216 ConstColIterator nativeColIt = nativeMatrix[nativeRowIdx].begin();
217 const ConstColIterator &nativeColEndIt = nativeMatrix[nativeRowIdx].end();
218 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
219 int domesticColIdx = overlap_->nativeToDomestic(nativeColIt.index());
224 if (domesticColIdx < 0)
225 domesticColIdx = overlap_->blackList().nativeToDomestic(nativeColIt.index());
227 if (domesticColIdx < 0)
233 (*this)[domesticRowIdx][domesticColIdx] = *nativeColIt;
238 void build_(
const BCRSMatrix &nativeMatrix)
240 int numDomestic = overlap_->numDomestic();
243 this->setSize(numDomestic, numDomestic);
244 this->setBuildMode(ParentType::random);
247 buildIndices_(nativeMatrix);
250 int numDomesticEntriesInRowFor_(
const BCRSMatrix &nativeMatrix,
int peerRank,
int rowIdx)
254 typedef typename BCRSMatrix::ConstColIterator ColIt;
255 ColIt colIt = nativeMatrix[rowIdx].begin();
256 ColIt colEndIt = nativeMatrix[rowIdx].end();
257 for (; colIt != colEndIt; ++colIt) {
258 if (overlap_->isDomesticIndexFor(peerRank, colIt.index()))
265 void buildIndices_(
const BCRSMatrix &nativeMatrix)
270 entries_.resize(overlap_->numDomestic());
271 for (
unsigned nativeRowIdx = 0; nativeRowIdx < nativeMatrix.N(); ++nativeRowIdx) {
272 int domesticRowIdx = overlap_->nativeToDomestic(nativeRowIdx);
273 if (domesticRowIdx < 0)
276 ConstColIterator nativeColIt = nativeMatrix[nativeRowIdx].begin();
277 ConstColIterator nativeColEndIt = nativeMatrix[nativeRowIdx].end();
278 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
279 int domesticColIdx = overlap_->nativeToDomestic(nativeColIt.index());
284 if (domesticColIdx < 0) {
285 domesticColIdx = overlap_->blackList().nativeToDomestic(nativeColIt.index());
288 if (domesticColIdx < 0)
291 entries_[domesticRowIdx].insert(domesticColIdx);
300 const PeerSet &peerSet = overlap_->peerSet();
301 typename PeerSet::const_iterator peerIt = peerSet.begin();
302 typename PeerSet::const_iterator peerEndIt = peerSet.end();
303 for (; peerIt != peerEndIt; ++peerIt) {
304 int peerRank = *peerIt;
305 sendIndices_(nativeMatrix, peerRank);
309 peerIt = peerSet.begin();
310 for (; peerIt != peerEndIt; ++peerIt) {
311 int peerRank = *peerIt;
312 receiveIndices_(peerRank);
316 peerIt = peerSet.begin();
317 for (; peerIt != peerEndIt; ++peerIt) {
318 int peerRank = *peerIt;
320 numRowsSendBuff_[peerRank]->wait();
321 rowSizesSendBuff_[peerRank]->wait();
322 rowIndicesSendBuff_[peerRank]->wait();
323 entryColIndicesSendBuff_[peerRank]->wait();
327 globalToDomesticBuff_(*rowIndicesSendBuff_[peerRank]);
328 globalToDomesticBuff_(*entryColIndicesSendBuff_[peerRank]);
336 int numDomestic = overlap_->numDomestic();
337 for (
int rowIdx = 0; rowIdx < numDomestic; ++rowIdx) {
339 const auto &colIndices = entries_[rowIdx];
340 auto colIdxIt = colIndices.begin();
341 const auto &colIdxEndIt = colIndices.end();
342 for (; colIdxIt != colIdxEndIt; ++colIdxIt) {
350 this->setrowsize(rowIdx, numCols);
355 for (
int rowIdx = 0; rowIdx < numDomestic; ++rowIdx) {
356 const auto &colIndices = entries_[rowIdx];
358 auto colIdxIt = colIndices.begin();
359 const auto &colIdxEndIt = colIndices.end();
360 for (; colIdxIt != colIdxEndIt; ++colIdxIt) {
365 this->addindex(rowIdx, *colIdxIt);
375 void sendIndices_(
const BCRSMatrix &nativeMatrix,
int peerRank)
379 int numOverlapRows = overlap_->foreignOverlapSize(peerRank);
380 numRowsSendBuff_[peerRank] =
new MpiBuffer<Index>(1);
381 (*numRowsSendBuff_[peerRank])[0] = numOverlapRows;
382 numRowsSendBuff_[peerRank]->send(peerRank);
386 rowIndicesSendBuff_[peerRank] =
new MpiBuffer<Index>(numOverlapRows);
387 rowSizesSendBuff_[peerRank] =
new MpiBuffer<Index>(numOverlapRows);
390 typedef std::set<int> ColumnIndexSet;
391 typedef std::map<int, ColumnIndexSet> EntryTuples;
393 EntryTuples entryIndices;
395 for (
int overlapOffset = 0; overlapOffset < numOverlapRows; ++overlapOffset) {
396 int domesticRowIdx = overlap_->foreignOverlapOffsetToDomesticIdx(peerRank, overlapOffset);
397 int nativeRowIdx = overlap_->domesticToNative(domesticRowIdx);
398 int globalRowIdx = overlap_->domesticToGlobal(domesticRowIdx);
400 ColumnIndexSet &colIndices = entryIndices[globalRowIdx];
402 typedef typename BCRSMatrix::ConstColIterator ColIt;
403 ColIt nativeColIt = nativeMatrix[nativeRowIdx].begin();
404 ColIt nativeColEndIt = nativeMatrix[nativeRowIdx].end();
405 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
406 int nativeColIdx = nativeColIt.index();
407 int domesticColIdx = overlap_->nativeToDomestic(nativeColIdx);
409 if (domesticColIdx < 0)
412 domesticColIdx = overlap_->blackList().nativeToDomestic(nativeColIdx);
414 if (domesticColIdx < 0)
420 int globalColIdx = overlap_->domesticToGlobal(domesticColIdx);
421 colIndices.insert(globalColIdx);
427 entryColIndicesSendBuff_[peerRank] =
new MpiBuffer<Index>(numEntries);
428 int overlapEntryIdx = 0;
429 for (
int overlapOffset = 0; overlapOffset < numOverlapRows; ++overlapOffset) {
430 int domesticRowIdx = overlap_->foreignOverlapOffsetToDomesticIdx(peerRank, overlapOffset);
431 int globalRowIdx = overlap_->domesticToGlobal(domesticRowIdx);
433 (*rowIndicesSendBuff_[peerRank])[overlapOffset] = globalRowIdx;
435 const ColumnIndexSet &colIndexSet = entryIndices[globalRowIdx];
436 auto* rssb = rowSizesSendBuff_[peerRank];
437 (*rssb)[overlapOffset] = colIndexSet.size();
438 for (
auto it = colIndexSet.begin(); it != colIndexSet.end(); ++it) {
439 int globalColIdx = *it;
441 (*entryColIndicesSendBuff_[peerRank])[overlapEntryIdx] = globalColIdx;
447 rowSizesSendBuff_[peerRank]->send(peerRank);
448 rowIndicesSendBuff_[peerRank]->send(peerRank);
449 entryColIndicesSendBuff_[peerRank]->send(peerRank);
453 entryValuesSendBuff_[peerRank] =
new MpiBuffer<block_type>(numEntries);
458 void receiveIndices_(
int peerRank)
462 Index numOverlapRows;
463 auto &numRowsRecvBuff = numRowsRecvBuff_[peerRank];
464 numRowsRecvBuff.resize(1);
465 numRowsRecvBuff.receive(peerRank);
466 numOverlapRows = numRowsRecvBuff[0];
470 rowSizesRecvBuff_[peerRank] =
new MpiBuffer<Index>(numOverlapRows);
471 rowIndicesRecvBuff_[peerRank] =
new MpiBuffer<Index>(numOverlapRows);
472 rowSizesRecvBuff_[peerRank]->receive(peerRank);
473 rowIndicesRecvBuff_[peerRank]->receive(peerRank);
477 int totalIndices = 0;
478 for (
Index i = 0; i < numOverlapRows; ++i) {
479 totalIndices += (*rowSizesRecvBuff_[peerRank])[i];
483 entryColIndicesRecvBuff_[peerRank] =
new MpiBuffer<Index>(totalIndices);
484 entryValuesRecvBuff_[peerRank] =
new MpiBuffer<block_type>(totalIndices);
487 entryColIndicesRecvBuff_[peerRank]->receive(peerRank);
491 globalToDomesticBuff_(*rowIndicesRecvBuff_[peerRank]);
492 globalToDomesticBuff_(*entryColIndicesRecvBuff_[peerRank]);
496 for (
Index i = 0; i < numOverlapRows; ++i) {
497 int domRowIdx = (*rowIndicesRecvBuff_[peerRank])[i];
498 for (
Index j = 0; j < (*rowSizesRecvBuff_[peerRank])[i]; ++j) {
499 int domColIdx = (*entryColIndicesRecvBuff_[peerRank])[k];
500 entries_[domRowIdx].insert(domColIdx);
511 const PeerSet &peerSet = overlap_->peerSet();
512 typename PeerSet::const_iterator peerIt = peerSet.begin();
513 typename PeerSet::const_iterator peerEndIt = peerSet.end();
514 for (; peerIt != peerEndIt; ++peerIt) {
515 int peerRank = *peerIt;
517 sendEntries_(peerRank);
521 peerIt = peerSet.begin();
522 for (; peerIt != peerEndIt; ++peerIt) {
523 int peerRank = *peerIt;
525 receiveAddEntries_(peerRank);
530 peerIt = peerSet.begin();
531 for (; peerIt != peerEndIt; ++peerIt) {
532 int peerRank = *peerIt;
533 entryValuesSendBuff_[peerRank]->wait();
542 const PeerSet &peerSet = overlap_->peerSet();
543 typename PeerSet::const_iterator peerIt = peerSet.begin();
544 typename PeerSet::const_iterator peerEndIt = peerSet.end();
545 for (; peerIt != peerEndIt; ++peerIt) {
546 int peerRank = *peerIt;
548 sendEntries_(peerRank);
552 peerIt = peerSet.begin();
553 for (; peerIt != peerEndIt; ++peerIt) {
554 int peerRank = *peerIt;
556 receiveCopyEntries_(peerRank);
561 peerIt = peerSet.begin();
562 for (; peerIt != peerEndIt; ++peerIt) {
563 int peerRank = *peerIt;
564 entryValuesSendBuff_[peerRank]->wait();
568 void sendEntries_(
int peerRank)
571 auto &mpiSendBuff = *entryValuesSendBuff_[peerRank];
573 auto &mpiRowIndicesSendBuff = *rowIndicesSendBuff_[peerRank];
574 auto &mpiRowSizesSendBuff = *rowSizesSendBuff_[peerRank];
575 auto &mpiColIndicesSendBuff = *entryColIndicesSendBuff_[peerRank];
579 for (
unsigned i = 0; i < mpiRowIndicesSendBuff.size(); ++i) {
580 Index domRowIdx = mpiRowIndicesSendBuff[i];
582 for (
Index j = 0; j < static_cast<Index>(mpiRowSizesSendBuff[i]); ++j)
585 Index domColIdx = mpiColIndicesSendBuff[k];
588 mpiSendBuff[k] = (*this)[domRowIdx][domColIdx];
593 mpiSendBuff.send(peerRank);
597 void receiveAddEntries_(
int peerRank)
600 auto &mpiRecvBuff = *entryValuesRecvBuff_[peerRank];
602 auto &mpiRowIndicesRecvBuff = *rowIndicesRecvBuff_[peerRank];
603 auto &mpiRowSizesRecvBuff = *rowSizesRecvBuff_[peerRank];
604 auto &mpiColIndicesRecvBuff = *entryColIndicesRecvBuff_[peerRank];
606 mpiRecvBuff.receive(peerRank);
610 for (
Index i = 0; i < static_cast<Index>(mpiRowIndicesRecvBuff.size()); ++i) {
611 Index domRowIdx = mpiRowIndicesRecvBuff[i];
612 for (
Index j = 0; j < static_cast<Index>(mpiRowSizesRecvBuff[i]); ++j, ++k) {
613 Index domColIdx = mpiColIndicesRecvBuff[k];
619 (*this)[domRowIdx][domColIdx] += mpiRecvBuff[k];
625 void receiveCopyEntries_(
int peerRank)
628 MpiBuffer<block_type> &mpiRecvBuff = *entryValuesRecvBuff_[peerRank];
630 MpiBuffer<Index> &mpiRowIndicesRecvBuff = *rowIndicesRecvBuff_[peerRank];
631 MpiBuffer<Index> &mpiRowSizesRecvBuff = *rowSizesRecvBuff_[peerRank];
632 MpiBuffer<Index> &mpiColIndicesRecvBuff = *entryColIndicesRecvBuff_[peerRank];
634 mpiRecvBuff.receive(peerRank);
638 for (
int i = 0; i < mpiRowIndicesRecvBuff.size(); ++i) {
639 Index domRowIdx = mpiRowIndicesRecvBuff[i];
640 for (
int j = 0; j < mpiRowSizesRecvBuff[i]; ++j, ++k) {
641 Index domColIdx = mpiColIndicesRecvBuff[k];
647 (*this)[domRowIdx][domColIdx] = mpiRecvBuff[k];
653 void globalToDomesticBuff_(MpiBuffer<Index> &idxBuff)
655 for (
unsigned i = 0; i < idxBuff.size(); ++i)
656 idxBuff[i] = overlap_->globalToDomestic(idxBuff[i]);
661 std::shared_ptr<Overlap> overlap_;
663 std::map<ProcessRank, MpiBuffer<Index> *> numRowsSendBuff_;
664 std::map<ProcessRank, MpiBuffer<Index> *> rowSizesSendBuff_;
665 std::map<ProcessRank, MpiBuffer<Index> *> rowIndicesSendBuff_;
666 std::map<ProcessRank, MpiBuffer<Index> *> entryColIndicesSendBuff_;
667 std::map<ProcessRank, MpiBuffer<block_type> *> entryValuesSendBuff_;
669 std::map<ProcessRank, MpiBuffer<Index> > numRowsRecvBuff_;
670 std::map<ProcessRank, MpiBuffer<Index> *> rowSizesRecvBuff_;
671 std::map<ProcessRank, MpiBuffer<Index> *> rowIndicesRecvBuff_;
672 std::map<ProcessRank, MpiBuffer<Index> *> entryColIndicesRecvBuff_;
673 std::map<ProcessRank, MpiBuffer<block_type> *> entryValuesRecvBuff_;
A set of process ranks.
Definition: overlaptypes.hh:147
Expresses which degrees of freedom are blacklisted for the parallel linear solvers and which domestic...
Definition: blacklist.hh:45
An overlap aware block-compressed row storage (BCRS) matrix.
Definition: overlappingbcrsmatrix.hh:52
Simplifies handling of buffers to be used in conjunction with MPI.
void assignCopy(const BCRSMatrix &nativeMatrix)
Assign and syncronize the overlapping matrix from a non-overlapping one.
Definition: overlappingbcrsmatrix.hh:138
std::list< BorderIndex > BorderList
This class managages a list of indices which are on the border of a process' partition of the grid...
Definition: overlaptypes.hh:119
OverlappingBCRSMatrix(const OverlappingBCRSMatrix &other)
Definition: overlappingbcrsmatrix.hh:68
Definition: baseauxiliarymodule.hh:35
const Overlap & overlap() const
Returns the domestic overlap for the process.
Definition: overlappingbcrsmatrix.hh:117
OverlappingBCRSMatrix(const BCRSMatrix &nativeMatrix, const BorderList &borderList, const BlackList &blackList, int overlapSize)
Definition: overlappingbcrsmatrix.hh:72
void assignAdd(const BCRSMatrix &nativeMatrix)
Assign and syncronize the overlapping matrix from a non-overlapping one.
Definition: overlappingbcrsmatrix.hh:123
~OverlappingBCRSMatrix()
Definition: overlappingbcrsmatrix.hh:89
ParentType::block_type block_type
Definition: overlappingbcrsmatrix.hh:65
This class creates and manages the foreign overlap given an initial list of border indices and a BCRS...
Definition: domesticoverlapfrombcrsmatrix.hh:51
int Index
The type of an index of a degree of freedom.
Definition: overlaptypes.hh:43
This class maps domestic row indices to and from "global" indices which is used to construct an algeb...
void print() const
Definition: overlappingbcrsmatrix.hh:168
Ewoms::Linear::DomesticOverlapFromBCRSMatrix< BCRSMatrix > Overlap
Definition: overlappingbcrsmatrix.hh:57
ParentType::ColIterator ColIterator
Definition: overlappingbcrsmatrix.hh:63
ParentType::ConstColIterator ConstColIterator
Definition: overlappingbcrsmatrix.hh:64
Expresses which degrees of freedom are blacklisted for the parallel linear solvers and which domestic...
void resetFront()
Set the identity matrix on the main diagonal of front indices.
Definition: overlappingbcrsmatrix.hh:150
This class creates and manages the foreign overlap given an initial list of border indices and a BCRS...