27#ifndef EWOMS_OVERLAPPING_BCRS_MATRIX_HH
28#define EWOMS_OVERLAPPING_BCRS_MATRIX_HH
35#include <opm/material/common/Valgrind.hpp>
37#include <dune/istl/scalarproducts.hh>
38#include <dune/istl/io.hh>
52template <
class BCRSMatrix>
55 using ParentType = BCRSMatrix;
61 using Entries = std::vector<std::set<Index> >;
74 template <
class NativeBCRSMatrix>
80 overlap_ = std::make_shared<Overlap>(nativeMatrix, borderList, blackList, overlapSize);
83 MPI_Comm_rank(MPI_COMM_WORLD, &myRank_);
95 typename BCRSMatrix::BuildMode)
96 {
throw std::logic_error(
"OverlappingBCRSMatrix objects cannot be build from scratch!"); }
100 if (overlap_.use_count() == 0)
104 const PeerSet& peerSet = overlap_->peerSet();
105 typename PeerSet::const_iterator peerIt = peerSet.begin();
106 typename PeerSet::const_iterator peerEndIt = peerSet.end();
107 for (; peerIt != peerEndIt; ++peerIt) {
110 delete rowSizesRecvBuff_[peerRank];
111 delete rowIndicesRecvBuff_[peerRank];
112 delete entryColIndicesRecvBuff_[peerRank];
113 delete entryValuesRecvBuff_[peerRank];
115 delete numRowsSendBuff_[peerRank];
116 delete rowSizesSendBuff_[peerRank];
117 delete rowIndicesSendBuff_[peerRank];
118 delete entryColIndicesSendBuff_[peerRank];
119 delete entryValuesSendBuff_[peerRank];
133 {
return *overlap_; }
153 template <
class NativeBCRSMatrix>
170 for (
unsigned i = 0; i < idMatrix.size(); ++i)
171 idMatrix[i][i] = 1.0;
173 int numLocal = overlap_->numLocal();
174 int numDomestic = overlap_->numDomestic();
175 for (
int domRowIdx = numLocal; domRowIdx < numDomestic; ++domRowIdx) {
176 if (overlap_->isFront(domRowIdx)) {
178 (*this)[domRowIdx] = 0.0;
179 (*this)[domRowIdx][domRowIdx] = idMatrix;
188 for (
int i = 0; i < this->N(); ++i) {
189 if (overlap_->isLocal(i))
193 std::cout <<
"row " << i <<
" ";
195 using ColIt =
typename BCRSMatrix::ConstColIterator;
196 ColIt colIt = (*this)[i].begin();
197 ColIt colEndIt = (*this)[i].end();
198 for (
int j = 0; j < this->M(); ++j) {
199 if (colIt != colEndIt && j == colIt.index()) {
201 if (overlap_->isBorder(j))
203 else if (overlap_->isLocal(j))
211 std::cout <<
"\n" << std::flush;
213 Dune::printSparseMatrix(std::cout,
214 *
static_cast<const BCRSMatrix *
>(
this),
219 template <
class NativeBCRSMatrix>
223 BCRSMatrix::operator=(0.0);
226 for (
unsigned nativeRowIdx = 0; nativeRowIdx < nativeMatrix.N(); ++nativeRowIdx) {
227 Index domesticRowIdx = overlap_->nativeToDomestic(
static_cast<Index>(nativeRowIdx));
228 if (domesticRowIdx < 0) {
232 auto nativeColIt = nativeMatrix[nativeRowIdx].begin();
233 const auto& nativeColEndIt = nativeMatrix[nativeRowIdx].end();
234 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
235 Index domesticColIdx = overlap_->nativeToDomestic(
static_cast<Index>(nativeColIt.index()));
240 if (domesticColIdx < 0)
241 domesticColIdx = overlap_->blackList().nativeToDomestic(
static_cast<Index>(nativeColIt.index()));
243 if (domesticColIdx < 0)
252 const auto& src = *nativeColIt;
253 auto& dest = (*this)[
static_cast<unsigned>(domesticRowIdx)][
static_cast<unsigned>(domesticColIdx)];
254 for (
unsigned i = 0; i < src.rows; ++i) {
255 for (
unsigned j = 0; j < src.cols; ++j) {
256 dest[i][j] =
static_cast<field_type>(src[i][j]);
267 const PeerSet& peerSet = overlap_->peerSet();
268 typename PeerSet::const_iterator peerIt = peerSet.begin();
269 typename PeerSet::const_iterator peerEndIt = peerSet.end();
270 for (; peerIt != peerEndIt; ++peerIt) {
273 sendEntries_(peerRank);
277 peerIt = peerSet.begin();
278 for (; peerIt != peerEndIt; ++peerIt) {
281 receiveAddEntries_(peerRank);
286 peerIt = peerSet.begin();
287 for (; peerIt != peerEndIt; ++peerIt) {
289 entryValuesSendBuff_[peerRank]->wait();
298 const PeerSet& peerSet = overlap_->peerSet();
299 typename PeerSet::const_iterator peerIt = peerSet.begin();
300 typename PeerSet::const_iterator peerEndIt = peerSet.end();
301 for (; peerIt != peerEndIt; ++peerIt) {
304 sendEntries_(peerRank);
308 peerIt = peerSet.begin();
309 for (; peerIt != peerEndIt; ++peerIt) {
312 receiveCopyEntries_(peerRank);
317 peerIt = peerSet.begin();
318 for (; peerIt != peerEndIt; ++peerIt) {
320 entryValuesSendBuff_[peerRank]->wait();
325 template <
class NativeBCRSMatrix>
326 void build_(
const NativeBCRSMatrix& nativeMatrix)
328 size_t numDomestic = overlap_->numDomestic();
331 this->setSize(numDomestic, numDomestic);
332 this->setBuildMode(ParentType::random);
335 buildIndices_(nativeMatrix);
338 template <
class NativeBCRSMatrix>
339 void buildIndices_(
const NativeBCRSMatrix& nativeMatrix)
344 entries_.resize(overlap_->numDomestic());
345 for (
unsigned nativeRowIdx = 0; nativeRowIdx < nativeMatrix.N(); ++nativeRowIdx) {
346 int domesticRowIdx = overlap_->nativeToDomestic(
static_cast<Index>(nativeRowIdx));
347 if (domesticRowIdx < 0)
350 auto nativeColIt = nativeMatrix[nativeRowIdx].begin();
351 const auto& nativeColEndIt = nativeMatrix[nativeRowIdx].end();
352 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
353 int domesticColIdx = overlap_->nativeToDomestic(
static_cast<Index>(nativeColIt.index()));
358 if (domesticColIdx < 0) {
359 domesticColIdx = overlap_->blackList().nativeToDomestic(
static_cast<Index>(nativeColIt.index()));
362 if (domesticColIdx < 0)
365 entries_[
static_cast<unsigned>(domesticRowIdx)].insert(domesticColIdx);
374 const PeerSet& peerSet = overlap_->peerSet();
375 typename PeerSet::const_iterator peerIt = peerSet.begin();
376 typename PeerSet::const_iterator peerEndIt = peerSet.end();
377 for (; peerIt != peerEndIt; ++peerIt) {
379 sendIndices_(nativeMatrix, peerRank);
383 peerIt = peerSet.begin();
384 for (; peerIt != peerEndIt; ++peerIt) {
386 receiveIndices_(peerRank);
390 peerIt = peerSet.begin();
391 for (; peerIt != peerEndIt; ++peerIt) {
394 numRowsSendBuff_[peerRank]->wait();
395 rowSizesSendBuff_[peerRank]->wait();
396 rowIndicesSendBuff_[peerRank]->wait();
397 entryColIndicesSendBuff_[peerRank]->wait();
401 globalToDomesticBuff_(*rowIndicesSendBuff_[peerRank]);
402 globalToDomesticBuff_(*entryColIndicesSendBuff_[peerRank]);
410 size_t numDomestic = overlap_->numDomestic();
411 for (
unsigned rowIdx = 0; rowIdx < numDomestic; ++rowIdx) {
412 unsigned numCols = 0;
413 const auto& colIndices = entries_[rowIdx];
414 auto colIdxIt = colIndices.begin();
415 const auto& colIdxEndIt = colIndices.end();
416 for (; colIdxIt != colIdxEndIt; ++colIdxIt) {
424 this->setrowsize(rowIdx, numCols);
429 for (
unsigned rowIdx = 0; rowIdx < numDomestic; ++rowIdx) {
430 const auto& colIndices = entries_[rowIdx];
432 auto colIdxIt = colIndices.begin();
433 const auto& colIdxEndIt = colIndices.end();
434 for (; colIdxIt != colIdxEndIt; ++colIdxIt) {
439 this->addindex(rowIdx,
static_cast<unsigned>(*colIdxIt));
449 template <
class NativeBCRSMatrix>
450 void sendIndices_([[maybe_unused]]
const NativeBCRSMatrix& nativeMatrix,
455 size_t numOverlapRows = overlap_->foreignOverlapSize(peerRank);
456 numRowsSendBuff_[peerRank] =
new MpiBuffer<unsigned>(1);
457 (*numRowsSendBuff_[peerRank])[0] =
static_cast<unsigned>(numOverlapRows);
458 numRowsSendBuff_[peerRank]->send(peerRank);
462 rowIndicesSendBuff_[peerRank] =
new MpiBuffer<Index>(numOverlapRows);
463 rowSizesSendBuff_[peerRank] =
new MpiBuffer<unsigned>(numOverlapRows);
466 using ColumnIndexSet = std::set<int>;
467 using EntryTuples = std::map<int, ColumnIndexSet>;
469 EntryTuples entryIndices;
470 unsigned numEntries = 0;
471 for (
unsigned overlapOffset = 0; overlapOffset < numOverlapRows; ++overlapOffset) {
472 Index domesticRowIdx = overlap_->foreignOverlapOffsetToDomesticIdx(peerRank, overlapOffset);
473 Index nativeRowIdx = overlap_->domesticToNative(domesticRowIdx);
474 Index globalRowIdx = overlap_->domesticToGlobal(domesticRowIdx);
476 ColumnIndexSet& colIndices = entryIndices[globalRowIdx];
478 auto nativeColIt = nativeMatrix[
static_cast<unsigned>(nativeRowIdx)].begin();
479 const auto& nativeColEndIt = nativeMatrix[
static_cast<unsigned>(nativeRowIdx)].end();
480 for (; nativeColIt != nativeColEndIt; ++nativeColIt) {
481 unsigned nativeColIdx =
static_cast<unsigned>(nativeColIt.index());
482 Index domesticColIdx = overlap_->nativeToDomestic(
static_cast<Index>(nativeColIdx));
484 if (domesticColIdx < 0)
487 domesticColIdx = overlap_->blackList().nativeToDomestic(
static_cast<Index>(nativeColIdx));
489 if (domesticColIdx < 0)
495 Index globalColIdx = overlap_->domesticToGlobal(domesticColIdx);
496 colIndices.insert(globalColIdx);
502 entryColIndicesSendBuff_[peerRank] =
new MpiBuffer<Index>(numEntries);
503 Index overlapEntryIdx = 0;
504 for (
unsigned overlapOffset = 0; overlapOffset < numOverlapRows; ++overlapOffset) {
505 Index domesticRowIdx = overlap_->foreignOverlapOffsetToDomesticIdx(peerRank, overlapOffset);
506 Index globalRowIdx = overlap_->domesticToGlobal(domesticRowIdx);
508 (*rowIndicesSendBuff_[peerRank])[overlapOffset] = globalRowIdx;
510 const ColumnIndexSet& colIndexSet = entryIndices[globalRowIdx];
511 auto* rssb = rowSizesSendBuff_[peerRank];
512 (*rssb)[overlapOffset] =
static_cast<unsigned>(colIndexSet.size());
513 for (
auto it = colIndexSet.begin(); it != colIndexSet.end(); ++it) {
514 int globalColIdx = *it;
516 (*entryColIndicesSendBuff_[peerRank])[
static_cast<unsigned>(overlapEntryIdx)] = globalColIdx;
522 rowSizesSendBuff_[peerRank]->send(peerRank);
523 rowIndicesSendBuff_[peerRank]->send(peerRank);
524 entryColIndicesSendBuff_[peerRank]->send(peerRank);
528 entryValuesSendBuff_[peerRank] =
new MpiBuffer<block_type>(numEntries);
533 void receiveIndices_([[maybe_unused]]
ProcessRank peerRank)
537 unsigned numOverlapRows;
538 auto& numRowsRecvBuff = numRowsRecvBuff_[peerRank];
539 numRowsRecvBuff.resize(1);
540 numRowsRecvBuff.receive(peerRank);
541 numOverlapRows = numRowsRecvBuff[0];
545 rowSizesRecvBuff_[peerRank] =
new MpiBuffer<unsigned>(numOverlapRows);
546 rowIndicesRecvBuff_[peerRank] =
new MpiBuffer<Index>(numOverlapRows);
547 rowSizesRecvBuff_[peerRank]->receive(peerRank);
548 rowIndicesRecvBuff_[peerRank]->receive(peerRank);
552 unsigned totalIndices = 0;
553 for (
unsigned i = 0; i < numOverlapRows; ++i)
554 totalIndices += (*rowSizesRecvBuff_[peerRank])[i];
557 entryColIndicesRecvBuff_[peerRank] =
new MpiBuffer<Index>(totalIndices);
558 entryValuesRecvBuff_[peerRank] =
new MpiBuffer<block_type>(totalIndices);
561 entryColIndicesRecvBuff_[peerRank]->receive(peerRank);
565 globalToDomesticBuff_(*rowIndicesRecvBuff_[peerRank]);
566 globalToDomesticBuff_(*entryColIndicesRecvBuff_[peerRank]);
570 for (
unsigned i = 0; i < numOverlapRows; ++i) {
571 Index domRowIdx = (*rowIndicesRecvBuff_[peerRank])[i];
572 for (
unsigned j = 0; j < (*rowSizesRecvBuff_[peerRank])[i]; ++j) {
573 Index domColIdx = (*entryColIndicesRecvBuff_[peerRank])[k];
574 entries_[
static_cast<unsigned>(domRowIdx)].insert(domColIdx);
581 void sendEntries_([[maybe_unused]]
ProcessRank peerRank)
584 auto &mpiSendBuff = *entryValuesSendBuff_[peerRank];
586 auto &mpiRowIndicesSendBuff = *rowIndicesSendBuff_[peerRank];
587 auto &mpiRowSizesSendBuff = *rowSizesSendBuff_[peerRank];
588 auto &mpiColIndicesSendBuff = *entryColIndicesSendBuff_[peerRank];
592 for (
unsigned i = 0; i < mpiRowIndicesSendBuff.size(); ++i) {
593 Index domRowIdx = mpiRowIndicesSendBuff[i];
595 for (
Index j = 0; j < static_cast<Index>(mpiRowSizesSendBuff[i]); ++j)
598 Index domColIdx = mpiColIndicesSendBuff[k];
601 mpiSendBuff[k] = (*this)[
static_cast<unsigned>(domRowIdx)][
static_cast<unsigned>(domColIdx)];
606 mpiSendBuff.send(peerRank);
610 void receiveAddEntries_([[maybe_unused]]
ProcessRank peerRank)
613 auto &mpiRecvBuff = *entryValuesRecvBuff_[peerRank];
615 auto &mpiRowIndicesRecvBuff = *rowIndicesRecvBuff_[peerRank];
616 auto &mpiRowSizesRecvBuff = *rowSizesRecvBuff_[peerRank];
617 auto &mpiColIndicesRecvBuff = *entryColIndicesRecvBuff_[peerRank];
619 mpiRecvBuff.receive(peerRank);
623 for (
unsigned i = 0; i < mpiRowIndicesRecvBuff.size(); ++i) {
624 Index domRowIdx = mpiRowIndicesRecvBuff[i];
625 for (
unsigned j = 0; j < mpiRowSizesRecvBuff[i]; ++j, ++k) {
626 Index domColIdx = mpiColIndicesRecvBuff[k];
632 (*this)[
static_cast<unsigned>(domRowIdx)][
static_cast<unsigned>(domColIdx)] += mpiRecvBuff[k];
638 void receiveCopyEntries_([[maybe_unused]]
int peerRank)
641 MpiBuffer<block_type> &mpiRecvBuff = *entryValuesRecvBuff_[peerRank];
643 MpiBuffer<Index> &mpiRowIndicesRecvBuff = *rowIndicesRecvBuff_[peerRank];
644 MpiBuffer<unsigned> &mpiRowSizesRecvBuff = *rowSizesRecvBuff_[peerRank];
645 MpiBuffer<Index> &mpiColIndicesRecvBuff = *entryColIndicesRecvBuff_[peerRank];
647 mpiRecvBuff.receive(peerRank);
651 for (
unsigned i = 0; i < mpiRowIndicesRecvBuff.size(); ++i) {
652 Index domRowIdx = mpiRowIndicesRecvBuff[i];
653 for (
unsigned j = 0; j < mpiRowSizesRecvBuff[i]; ++j, ++k) {
654 Index domColIdx = mpiColIndicesRecvBuff[k];
660 (*this)[
static_cast<unsigned>(domRowIdx)][
static_cast<unsigned>(domColIdx)] = mpiRecvBuff[k];
666 void globalToDomesticBuff_(MpiBuffer<Index>& idxBuff)
668 for (
unsigned i = 0; i < idxBuff.size(); ++i)
669 idxBuff[i] = overlap_->globalToDomestic(idxBuff[i]);
674 std::shared_ptr<Overlap> overlap_;
676 std::map<ProcessRank, MpiBuffer<unsigned> *> numRowsSendBuff_;
677 std::map<ProcessRank, MpiBuffer<unsigned> *> rowSizesSendBuff_;
678 std::map<ProcessRank, MpiBuffer<Index> *> rowIndicesSendBuff_;
679 std::map<ProcessRank, MpiBuffer<Index> *> entryColIndicesSendBuff_;
680 std::map<ProcessRank, MpiBuffer<block_type> *> entryValuesSendBuff_;
682 std::map<ProcessRank, MpiBuffer<unsigned> > numRowsRecvBuff_;
683 std::map<ProcessRank, MpiBuffer<unsigned> *> rowSizesRecvBuff_;
684 std::map<ProcessRank, MpiBuffer<Index> *> rowIndicesRecvBuff_;
685 std::map<ProcessRank, MpiBuffer<Index> *> entryColIndicesRecvBuff_;
686 std::map<ProcessRank, MpiBuffer<block_type> *> entryValuesRecvBuff_;
Expresses which degrees of freedom are blacklisted for the parallel linear solvers and which domestic...
Definition: blacklist.hh:49
This class creates and manages the foreign overlap given an initial list of border indices and a BCRS...
Definition: domesticoverlapfrombcrsmatrix.hh:53
An overlap aware block-compressed row storage (BCRS) matrix.
Definition: overlappingbcrsmatrix.hh:54
typename ParentType::ColIterator ColIterator
Definition: overlappingbcrsmatrix.hh:64
void syncAdd()
Definition: overlappingbcrsmatrix.hh:264
const ParentType & asParent() const
Definition: overlappingbcrsmatrix.hh:126
typename ParentType::block_type block_type
Definition: overlappingbcrsmatrix.hh:66
void assignFromNative(const NativeBCRSMatrix &nativeMatrix)
Definition: overlappingbcrsmatrix.hh:220
typename ParentType::field_type field_type
Definition: overlappingbcrsmatrix.hh:67
OverlappingBCRSMatrix(const NativeBCRSMatrix &nativeMatrix, const BorderList &borderList, const BlackList &blackList, unsigned overlapSize)
Definition: overlappingbcrsmatrix.hh:75
void resetFront()
Set the identity matrix on the main diagonal of front indices.
Definition: overlappingbcrsmatrix.hh:166
void assignAdd(const ParentType &nativeMatrix)
Assign and syncronize the overlapping matrix from a non-overlapping one.
Definition: overlappingbcrsmatrix.hh:138
typename ParentType::ConstColIterator ConstColIterator
Definition: overlappingbcrsmatrix.hh:65
void assignCopy(const NativeBCRSMatrix &nativeMatrix)
Assign and syncronize the overlapping matrix from a non-overlapping one.
Definition: overlappingbcrsmatrix.hh:154
OverlappingBCRSMatrix(size_t, size_t, typename BCRSMatrix::BuildMode)
Definition: overlappingbcrsmatrix.hh:93
~OverlappingBCRSMatrix()
Definition: overlappingbcrsmatrix.hh:98
void syncCopy()
Definition: overlappingbcrsmatrix.hh:295
const Overlap & overlap() const
Returns the domestic overlap for the process.
Definition: overlappingbcrsmatrix.hh:132
ParentType & asParent()
Definition: overlappingbcrsmatrix.hh:123
OverlappingBCRSMatrix(const OverlappingBCRSMatrix &other)
Definition: overlappingbcrsmatrix.hh:70
void print() const
Definition: overlappingbcrsmatrix.hh:184
A set of process ranks.
Definition: overlaptypes.hh:149
unsigned ProcessRank
The type of the rank of a process.
Definition: overlaptypes.hh:49
int Index
The type of an index of a degree of freedom.
Definition: overlaptypes.hh:44
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:120
Definition: blackoilboundaryratevector.hh:37