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> 52 template <
class BCRSMatrix>
55 using ParentType = BCRSMatrix;
61 using Entries = std::vector<std::set<Index> >;
64 using ColIterator =
typename ParentType::ColIterator;
65 using ConstColIterator =
typename ParentType::ConstColIterator;
66 using block_type =
typename ParentType::block_type;
67 using field_type =
typename ParentType::field_type;
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];
123 ParentType& asParent()
126 const ParentType& asParent()
const 133 {
return *overlap_; }
141 assignFromNative(nativeMatrix);
153 template <
class NativeBCRSMatrix>
157 assignFromNative(nativeMatrix);
169 block_type idMatrix(0.0);
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>
220 void assignFromNative(
const NativeBCRSMatrix& nativeMatrix)
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,
451 [[maybe_unused]] ProcessRank peerRank)
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_;
A set of process ranks.
Definition: overlaptypes.hh:148
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...
void assignAdd(const ParentType &nativeMatrix)
Assign and syncronize the overlapping matrix from a non-overlapping one.
Definition: overlappingbcrsmatrix.hh:138
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: blackoilbioeffectsmodules.hh:45
This class creates and manages the foreign overlap given an initial list of border indices and a BCRS...
Definition: domesticoverlapfrombcrsmatrix.hh:52
const Overlap & overlap() const
Returns the domestic overlap for the process.
Definition: overlappingbcrsmatrix.hh:132
void resetFront()
Set the identity matrix on the main diagonal of front indices.
Definition: overlappingbcrsmatrix.hh:166
unsigned ProcessRank
The type of the rank of a process.
Definition: overlaptypes.hh:49
Expresses which degrees of freedom are blacklisted for the parallel linear solvers and which domestic...
Definition: blacklist.hh:48
This class creates and manages the foreign overlap given an initial list of border indices and a BCRS...
void assignCopy(const NativeBCRSMatrix &nativeMatrix)
Assign and syncronize the overlapping matrix from a non-overlapping one.
Definition: overlappingbcrsmatrix.hh:154
An overlap aware block-compressed row storage (BCRS) matrix.
Definition: overlappingbcrsmatrix.hh:53
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
int Index
The type of an index of a degree of freedom.
Definition: overlaptypes.hh:44