27 #ifndef EWOMS_OVERLAPPING_BCRS_MATRIX_HH
28 #define EWOMS_OVERLAPPING_BCRS_MATRIX_HH
35 #include <opm/common/Valgrind.hpp>
37 #include <dune/istl/scalarproducts.hh>
38 #include <dune/istl/io.hh>
53 template <
class BCRSMatrix>
56 typedef BCRSMatrix ParentType;
62 typedef std::vector<std::set<Index> > Entries;
65 typedef typename ParentType::ColIterator ColIterator;
66 typedef typename ParentType::ConstColIterator ConstColIterator;
67 typedef typename ParentType::block_type block_type;
68 typedef typename ParentType::field_type field_type;
75 template <
class NativeBCRSMatrix>
77 const BorderList& borderList,
81 overlap_ = std::make_shared<Overlap>(nativeMatrix, borderList, blackList, overlapSize);
84 MPI_Comm_rank(MPI_COMM_WORLD, &myRank_);
94 if (overlap_.use_count() == 0)
98 const PeerSet& peerSet = overlap_->peerSet();
99 typename PeerSet::const_iterator peerIt = peerSet.begin();
100 typename PeerSet::const_iterator peerEndIt = peerSet.end();
101 for (; peerIt != peerEndIt; ++peerIt) {
102 ProcessRank peerRank = *peerIt;
104 delete rowSizesRecvBuff_[peerRank];
105 delete rowIndicesRecvBuff_[peerRank];
106 delete entryColIndicesRecvBuff_[peerRank];
107 delete entryValuesRecvBuff_[peerRank];
109 delete numRowsSendBuff_[peerRank];
110 delete rowSizesSendBuff_[peerRank];
111 delete rowIndicesSendBuff_[peerRank];
112 delete entryColIndicesSendBuff_[peerRank];
113 delete entryValuesSendBuff_[peerRank];
117 ParentType& asParent()
120 const ParentType& asParent()
const
127 {
return *overlap_; }
132 template <
class NativeBCRSMatrix>
136 assignFromNative(nativeMatrix);
148 template <
class NativeBCRSMatrix>
152 assignFromNative(nativeMatrix);
164 block_type idMatrix(0.0);
165 for (
unsigned i = 0; i < idMatrix.size(); ++i)
166 idMatrix[i][i] = 1.0;
168 int numLocal = overlap_->numLocal();
169 int numDomestic = overlap_->numDomestic();
170 for (
int domRowIdx = numLocal; domRowIdx < numDomestic; ++domRowIdx) {
171 if (overlap_->isFront(domRowIdx)) {
173 (*this)[domRowIdx] = 0.0;
174 (*this)[domRowIdx][domRowIdx] = idMatrix;
183 for (
int i = 0; i < this->N(); ++i) {
184 if (overlap_->isLocal(i))
188 std::cout <<
"row " << i <<
" ";
190 typedef typename BCRSMatrix::ConstColIterator ColIt;
191 ColIt colIt = (*this)[i].begin();
192 ColIt colEndIt = (*this)[i].end();
193 for (
int j = 0; j < this->M(); ++j) {
194 if (colIt != colEndIt && j == colIt.index()) {
196 if (overlap_->isBorder(j))
198 else if (overlap_->isLocal(j))
206 std::cout <<
"\n" << std::flush;
208 Dune::printSparseMatrix(std::cout,
209 *static_cast<const BCRSMatrix *>(
this),
214 template <
class NativeBCRSMatrix>
215 void assignFromNative(
const NativeBCRSMatrix& nativeMatrix)
218 BCRSMatrix::operator=(0.0);
221 for (
unsigned nativeRowIdx = 0; nativeRowIdx < nativeMatrix.N(); ++nativeRowIdx) {
222 Index domesticRowIdx = overlap_->nativeToDomestic(static_cast<Index>(nativeRowIdx));
223 if (domesticRowIdx < 0) {
227 auto nativeColIt = nativeMatrix[nativeRowIdx].begin();
228 const auto& nativeColEndIt = nativeMatrix[nativeRowIdx].end();
229 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
230 Index domesticColIdx = overlap_->nativeToDomestic(static_cast<Index>(nativeColIt.index()));
235 if (domesticColIdx < 0)
236 domesticColIdx = overlap_->blackList().nativeToDomestic(static_cast<Index>(nativeColIt.index()));
238 if (domesticColIdx < 0)
247 const auto& src = *nativeColIt;
248 auto& dest = (*this)[
static_cast<unsigned>(domesticRowIdx)][static_cast<unsigned>(domesticColIdx)];
249 for (
unsigned i = 0; i < src.rows; ++i) {
250 for (
unsigned j = 0; j < src.cols; ++j) {
251 dest[i][j] =
static_cast<field_type
>(src[i][j]);
262 const PeerSet& peerSet = overlap_->peerSet();
263 typename PeerSet::const_iterator peerIt = peerSet.begin();
264 typename PeerSet::const_iterator peerEndIt = peerSet.end();
265 for (; peerIt != peerEndIt; ++peerIt) {
266 ProcessRank peerRank = *peerIt;
268 sendEntries_(peerRank);
272 peerIt = peerSet.begin();
273 for (; peerIt != peerEndIt; ++peerIt) {
274 ProcessRank peerRank = *peerIt;
276 receiveAddEntries_(peerRank);
281 peerIt = peerSet.begin();
282 for (; peerIt != peerEndIt; ++peerIt) {
283 ProcessRank peerRank = *peerIt;
284 entryValuesSendBuff_[peerRank]->wait();
293 const PeerSet& peerSet = overlap_->peerSet();
294 typename PeerSet::const_iterator peerIt = peerSet.begin();
295 typename PeerSet::const_iterator peerEndIt = peerSet.end();
296 for (; peerIt != peerEndIt; ++peerIt) {
297 ProcessRank peerRank = *peerIt;
299 sendEntries_(peerRank);
303 peerIt = peerSet.begin();
304 for (; peerIt != peerEndIt; ++peerIt) {
305 ProcessRank peerRank = *peerIt;
307 receiveCopyEntries_(peerRank);
312 peerIt = peerSet.begin();
313 for (; peerIt != peerEndIt; ++peerIt) {
314 ProcessRank peerRank = *peerIt;
315 entryValuesSendBuff_[peerRank]->wait();
320 template <
class NativeBCRSMatrix>
321 void build_(
const NativeBCRSMatrix& nativeMatrix)
323 size_t numDomestic = overlap_->numDomestic();
326 this->setSize(numDomestic, numDomestic);
327 this->setBuildMode(ParentType::random);
330 buildIndices_(nativeMatrix);
333 template <
class NativeBCRSMatrix>
334 void buildIndices_(
const NativeBCRSMatrix& nativeMatrix)
339 entries_.resize(overlap_->numDomestic());
340 for (
unsigned nativeRowIdx = 0; nativeRowIdx < nativeMatrix.N(); ++nativeRowIdx) {
341 int domesticRowIdx = overlap_->nativeToDomestic(static_cast<Index>(nativeRowIdx));
342 if (domesticRowIdx < 0)
345 auto nativeColIt = nativeMatrix[nativeRowIdx].begin();
346 const auto& nativeColEndIt = nativeMatrix[nativeRowIdx].end();
347 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
348 int domesticColIdx = overlap_->nativeToDomestic(static_cast<Index>(nativeColIt.index()));
353 if (domesticColIdx < 0) {
354 domesticColIdx = overlap_->blackList().nativeToDomestic(static_cast<Index>(nativeColIt.index()));
357 if (domesticColIdx < 0)
360 entries_[
static_cast<unsigned>(domesticRowIdx)].insert(domesticColIdx);
369 const PeerSet& peerSet = overlap_->peerSet();
370 typename PeerSet::const_iterator peerIt = peerSet.begin();
371 typename PeerSet::const_iterator peerEndIt = peerSet.end();
372 for (; peerIt != peerEndIt; ++peerIt) {
373 ProcessRank peerRank = *peerIt;
374 sendIndices_(nativeMatrix, peerRank);
378 peerIt = peerSet.begin();
379 for (; peerIt != peerEndIt; ++peerIt) {
380 ProcessRank peerRank = *peerIt;
381 receiveIndices_(peerRank);
385 peerIt = peerSet.begin();
386 for (; peerIt != peerEndIt; ++peerIt) {
387 ProcessRank peerRank = *peerIt;
389 numRowsSendBuff_[peerRank]->wait();
390 rowSizesSendBuff_[peerRank]->wait();
391 rowIndicesSendBuff_[peerRank]->wait();
392 entryColIndicesSendBuff_[peerRank]->wait();
396 globalToDomesticBuff_(*rowIndicesSendBuff_[peerRank]);
397 globalToDomesticBuff_(*entryColIndicesSendBuff_[peerRank]);
405 size_t numDomestic = overlap_->numDomestic();
406 for (
unsigned rowIdx = 0; rowIdx < numDomestic; ++rowIdx) {
407 unsigned numCols = 0;
408 const auto& colIndices = entries_[rowIdx];
409 auto colIdxIt = colIndices.begin();
410 const auto& colIdxEndIt = colIndices.end();
411 for (; colIdxIt != colIdxEndIt; ++colIdxIt) {
419 this->setrowsize(rowIdx, numCols);
424 for (
unsigned rowIdx = 0; rowIdx < numDomestic; ++rowIdx) {
425 const auto& colIndices = entries_[rowIdx];
427 auto colIdxIt = colIndices.begin();
428 const auto& colIdxEndIt = colIndices.end();
429 for (; colIdxIt != colIdxEndIt; ++colIdxIt) {
434 this->addindex(rowIdx, static_cast<unsigned>(*colIdxIt));
444 template <
class NativeBCRSMatrix>
445 void sendIndices_(
const NativeBCRSMatrix& nativeMatrix, ProcessRank peerRank)
449 size_t numOverlapRows = overlap_->foreignOverlapSize(peerRank);
450 numRowsSendBuff_[peerRank] =
new MpiBuffer<unsigned>(1);
451 (*numRowsSendBuff_[peerRank])[0] = static_cast<unsigned>(numOverlapRows);
452 numRowsSendBuff_[peerRank]->send(peerRank);
456 rowIndicesSendBuff_[peerRank] =
new MpiBuffer<Index>(numOverlapRows);
457 rowSizesSendBuff_[peerRank] =
new MpiBuffer<unsigned>(numOverlapRows);
460 typedef std::set<int> ColumnIndexSet;
461 typedef std::map<int, ColumnIndexSet> EntryTuples;
463 EntryTuples entryIndices;
464 unsigned numEntries = 0;
465 for (
unsigned overlapOffset = 0; overlapOffset < numOverlapRows; ++overlapOffset) {
466 Index domesticRowIdx = overlap_->foreignOverlapOffsetToDomesticIdx(peerRank, overlapOffset);
467 Index nativeRowIdx = overlap_->domesticToNative(domesticRowIdx);
468 Index globalRowIdx = overlap_->domesticToGlobal(domesticRowIdx);
470 ColumnIndexSet& colIndices = entryIndices[globalRowIdx];
472 auto nativeColIt = nativeMatrix[
static_cast<unsigned>(nativeRowIdx)].begin();
473 const auto& nativeColEndIt = nativeMatrix[
static_cast<unsigned>(nativeRowIdx)].end();
474 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
475 unsigned nativeColIdx =
static_cast<unsigned>(nativeColIt.index());
476 Index domesticColIdx = overlap_->nativeToDomestic(static_cast<Index>(nativeColIdx));
478 if (domesticColIdx < 0)
481 domesticColIdx = overlap_->blackList().nativeToDomestic(static_cast<Index>(nativeColIdx));
483 if (domesticColIdx < 0)
489 Index globalColIdx = overlap_->domesticToGlobal(domesticColIdx);
490 colIndices.insert(globalColIdx);
496 entryColIndicesSendBuff_[peerRank] =
new MpiBuffer<Index>(numEntries);
497 Index overlapEntryIdx = 0;
498 for (
unsigned overlapOffset = 0; overlapOffset < numOverlapRows; ++overlapOffset) {
499 Index domesticRowIdx = overlap_->foreignOverlapOffsetToDomesticIdx(peerRank, overlapOffset);
500 Index globalRowIdx = overlap_->domesticToGlobal(domesticRowIdx);
502 (*rowIndicesSendBuff_[peerRank])[overlapOffset] = globalRowIdx;
504 const ColumnIndexSet& colIndexSet = entryIndices[globalRowIdx];
505 auto* rssb = rowSizesSendBuff_[peerRank];
506 (*rssb)[overlapOffset] =
static_cast<unsigned>(colIndexSet.size());
507 for (
auto it = colIndexSet.begin(); it != colIndexSet.end(); ++it) {
508 int globalColIdx = *it;
510 (*entryColIndicesSendBuff_[peerRank])[static_cast<unsigned>(overlapEntryIdx)] = globalColIdx;
516 rowSizesSendBuff_[peerRank]->send(peerRank);
517 rowIndicesSendBuff_[peerRank]->send(peerRank);
518 entryColIndicesSendBuff_[peerRank]->send(peerRank);
522 entryValuesSendBuff_[peerRank] =
new MpiBuffer<block_type>(numEntries);
527 void receiveIndices_(ProcessRank peerRank)
531 unsigned numOverlapRows;
532 auto& numRowsRecvBuff = numRowsRecvBuff_[peerRank];
533 numRowsRecvBuff.resize(1);
534 numRowsRecvBuff.receive(peerRank);
535 numOverlapRows = numRowsRecvBuff[0];
539 rowSizesRecvBuff_[peerRank] =
new MpiBuffer<unsigned>(numOverlapRows);
540 rowIndicesRecvBuff_[peerRank] =
new MpiBuffer<Index>(numOverlapRows);
541 rowSizesRecvBuff_[peerRank]->receive(peerRank);
542 rowIndicesRecvBuff_[peerRank]->receive(peerRank);
546 unsigned totalIndices = 0;
547 for (
unsigned i = 0; i < numOverlapRows; ++i)
548 totalIndices += (*rowSizesRecvBuff_[peerRank])[i];
551 entryColIndicesRecvBuff_[peerRank] =
new MpiBuffer<Index>(totalIndices);
552 entryValuesRecvBuff_[peerRank] =
new MpiBuffer<block_type>(totalIndices);
555 entryColIndicesRecvBuff_[peerRank]->receive(peerRank);
559 globalToDomesticBuff_(*rowIndicesRecvBuff_[peerRank]);
560 globalToDomesticBuff_(*entryColIndicesRecvBuff_[peerRank]);
564 for (
unsigned i = 0; i < numOverlapRows; ++i) {
565 Index domRowIdx = (*rowIndicesRecvBuff_[peerRank])[i];
566 for (
unsigned j = 0; j < (*rowSizesRecvBuff_[peerRank])[i]; ++j) {
567 Index domColIdx = (*entryColIndicesRecvBuff_[peerRank])[k];
568 entries_[
static_cast<unsigned>(domRowIdx)].insert(domColIdx);
575 void sendEntries_(ProcessRank peerRank)
578 auto &mpiSendBuff = *entryValuesSendBuff_[peerRank];
580 auto &mpiRowIndicesSendBuff = *rowIndicesSendBuff_[peerRank];
581 auto &mpiRowSizesSendBuff = *rowSizesSendBuff_[peerRank];
582 auto &mpiColIndicesSendBuff = *entryColIndicesSendBuff_[peerRank];
586 for (
unsigned i = 0; i < mpiRowIndicesSendBuff.size(); ++i) {
587 Index domRowIdx = mpiRowIndicesSendBuff[i];
589 for (Index j = 0; j < static_cast<Index>(mpiRowSizesSendBuff[i]); ++j)
592 Index domColIdx = mpiColIndicesSendBuff[k];
595 mpiSendBuff[k] = (*this)[
static_cast<unsigned>(domRowIdx)][static_cast<unsigned>(domColIdx)];
600 mpiSendBuff.send(peerRank);
604 void receiveAddEntries_(ProcessRank peerRank)
607 auto &mpiRecvBuff = *entryValuesRecvBuff_[peerRank];
609 auto &mpiRowIndicesRecvBuff = *rowIndicesRecvBuff_[peerRank];
610 auto &mpiRowSizesRecvBuff = *rowSizesRecvBuff_[peerRank];
611 auto &mpiColIndicesRecvBuff = *entryColIndicesRecvBuff_[peerRank];
613 mpiRecvBuff.receive(peerRank);
617 for (
unsigned i = 0; i < mpiRowIndicesRecvBuff.size(); ++i) {
618 Index domRowIdx = mpiRowIndicesRecvBuff[i];
619 for (
unsigned j = 0; j < mpiRowSizesRecvBuff[i]; ++j, ++k) {
620 Index domColIdx = mpiColIndicesRecvBuff[k];
626 (*this)[
static_cast<unsigned>(domRowIdx)][static_cast<unsigned>(domColIdx)] += mpiRecvBuff[k];
632 void receiveCopyEntries_(
int peerRank)
635 MpiBuffer<block_type> &mpiRecvBuff = *entryValuesRecvBuff_[peerRank];
637 MpiBuffer<Index> &mpiRowIndicesRecvBuff = *rowIndicesRecvBuff_[peerRank];
638 MpiBuffer<unsigned> &mpiRowSizesRecvBuff = *rowSizesRecvBuff_[peerRank];
639 MpiBuffer<Index> &mpiColIndicesRecvBuff = *entryColIndicesRecvBuff_[peerRank];
641 mpiRecvBuff.receive(peerRank);
645 for (
unsigned i = 0; i < mpiRowIndicesRecvBuff.size(); ++i) {
646 Index domRowIdx = mpiRowIndicesRecvBuff[i];
647 for (
unsigned j = 0; j < mpiRowSizesRecvBuff[i]; ++j, ++k) {
648 Index domColIdx = mpiColIndicesRecvBuff[k];
654 (*this)[
static_cast<unsigned>(domRowIdx)][static_cast<unsigned>(domColIdx)] = mpiRecvBuff[k];
660 void globalToDomesticBuff_(MpiBuffer<Index>& idxBuff)
662 for (
unsigned i = 0; i < idxBuff.size(); ++i)
663 idxBuff[i] = overlap_->globalToDomestic(idxBuff[i]);
668 std::shared_ptr<Overlap> overlap_;
670 std::map<ProcessRank, MpiBuffer<unsigned> *> numRowsSendBuff_;
671 std::map<ProcessRank, MpiBuffer<unsigned> *> rowSizesSendBuff_;
672 std::map<ProcessRank, MpiBuffer<Index> *> rowIndicesSendBuff_;
673 std::map<ProcessRank, MpiBuffer<Index> *> entryColIndicesSendBuff_;
674 std::map<ProcessRank, MpiBuffer<block_type> *> entryValuesSendBuff_;
676 std::map<ProcessRank, MpiBuffer<unsigned> > numRowsRecvBuff_;
677 std::map<ProcessRank, MpiBuffer<unsigned> *> rowSizesRecvBuff_;
678 std::map<ProcessRank, MpiBuffer<Index> *> rowIndicesRecvBuff_;
679 std::map<ProcessRank, MpiBuffer<Index> *> entryColIndicesRecvBuff_;
680 std::map<ProcessRank, MpiBuffer<block_type> *> entryValuesRecvBuff_;
void assignCopy(const NativeBCRSMatrix &nativeMatrix)
Assign and syncronize the overlapping matrix from a non-overlapping one.
Definition: overlappingbcrsmatrix.hh:149
Expresses which degrees of freedom are blacklisted for the parallel linear solvers and which domestic...
Simplifies handling of buffers to be used in conjunction with MPI.
This class maps domestic row indices to and from "global" indices which is used to construct an algeb...
const Overlap & overlap() const
Returns the domestic overlap for the process.
Definition: overlappingbcrsmatrix.hh:126
void assignAdd(const NativeBCRSMatrix &nativeMatrix)
Assign and syncronize the overlapping matrix from a non-overlapping one.
Definition: overlappingbcrsmatrix.hh:133
An overlap aware block-compressed row storage (BCRS) matrix.
Definition: overlappingbcrsmatrix.hh:54
Expresses which degrees of freedom are blacklisted for the parallel linear solvers and which domestic...
Definition: blacklist.hh:47
This class creates and manages the foreign overlap given an initial list of border indices and a BCRS...
This class creates and manages the foreign overlap given an initial list of border indices and a BCRS...
Definition: domesticoverlapfrombcrsmatrix.hh:52
void resetFront()
Set the identity matrix on the main diagonal of front indices.
Definition: overlappingbcrsmatrix.hh:161
A set of process ranks.
Definition: overlaptypes.hh:148