My Project
Loading...
Searching...
No Matches
CSRGraphFromCoordinates_impl.hpp
1/*
2 Copyright 2016 SINTEF ICT, Applied Mathematics.
3 Copyright 2016 Statoil ASA.
4 Copyright 2022 Equinor ASA
5
6 This file is part of the Open Porous Media Project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include <algorithm>
23#include <cassert>
24#include <exception>
25#include <iterator>
26#include <optional>
27#include <stdexcept>
28#include <string>
29#include <utility>
30#include <vector>
31
32// ---------------------------------------------------------------------
33// Class Opm::utility::CSRGraphFromCoordinates::Connections
34// ---------------------------------------------------------------------
35
36template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
37void
39Connections::add(const VertexID v1, const VertexID v2)
40{
41 this->i_.push_back(v1);
42 this->j_.push_back(v2);
43
44 this->max_i_ = std::max(this->max_i_.value_or(BaseVertexID{}), this->i_.back());
45 this->max_j_ = std::max(this->max_j_.value_or(BaseVertexID{}), this->j_.back());
46}
47
48template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
49void
51Connections::add(VertexID maxRowIdx,
52 VertexID maxColIdx,
53 const Neighbours& rows,
54 const Neighbours& cols)
55{
56 if (cols.size() != rows.size()) {
57 throw std::invalid_argument {
58 "Coordinate format column index table size does not match "
59 "row index table size"
60 };
61 }
62
63 this->i_.insert(this->i_.end(), rows .begin(), rows .end());
64 this->j_.insert(this->j_.end(), cols .begin(), cols .end());
65
66 this->max_i_ = std::max(this->max_i_.value_or(BaseVertexID{}), maxRowIdx);
67 this->max_j_ = std::max(this->max_j_.value_or(BaseVertexID{}), maxColIdx);
68}
69
70template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
71void
75 this->j_.clear();
76 this->i_.clear();
77
78 this->max_i_.reset();
79 this->max_j_.reset();
80}
81
82template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
83bool
86{
87 return this->i_.empty();
88}
89
90template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
91bool
94{
95 return this->i_.size() == this->j_.size();
96}
97
98template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
99std::optional<typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx, PermitSelfConnections>::BaseVertexID>
102{
103 return this->max_i_;
105
106template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
107std::optional<typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx, PermitSelfConnections>::BaseVertexID>
111 return this->max_j_;
112}
113
114template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
118{
119 return this->i_.size();
120}
121
122template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
126{
127 return this->i_;
128}
129
130template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
134{
135 return this->j_;
136}
137
138// =====================================================================
139
140// ---------------------------------------------------------------------
141// Class Opm::utility::CSRGraphFromCoordinates::CSR
142// ---------------------------------------------------------------------
143
144template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
145void
147CSR::merge(const Connections& conns,
148 const Offset maxNumVertices,
149 const bool expandExistingIdxMap)
150{
151 const auto maxRow = conns.maxRow();
152
153 if (maxRow.has_value() &&
154 (static_cast<Offset>(*maxRow) >= maxNumVertices))
155 {
156 throw std::invalid_argument {
157 "Number of vertices in input graph (" +
158 std::to_string(*maxRow) + ") "
159 "exceeds maximum graph size implied by explicit size of "
160 "adjacency matrix (" + std::to_string(maxNumVertices) + ')'
161 };
162 }
163
164 this->assemble(conns.rowIndices(), conns.columnIndices(),
165 maxRow.value_or(BaseVertexID{0}),
166 conns.maxCol().value_or(BaseVertexID{0}),
167 expandExistingIdxMap);
168
169 this->compress(maxNumVertices);
170}
171
172template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
175CSR::numRows() const
176{
177 return this->startPointers().empty()
178 ? 0 : this->startPointers().size() - 1;
179}
180
181template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
182typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx, PermitSelfConnections>::BaseVertexID
184CSR::maxRowID() const
185{
186 return this->numRows_ - 1;
187}
188
189template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
190typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx, PermitSelfConnections>::BaseVertexID
192CSR::maxColID() const
193{
194 return this->numCols_ - 1;
195}
196
197template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
200CSR::startPointers() const
201{
202 return this->ia_;
203}
204
205template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
208CSR::columnIndices() const
209{
210 return this->ja_;
211}
212
213template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
217{
218 auto rowIdx = Neighbours{};
219
220 if (this->ia_.empty()) {
221 return rowIdx;
222 }
223
224 rowIdx.reserve(this->ia_.back());
225
226 auto row = BaseVertexID{};
227
228 const auto m = this->ia_.size() - 1;
229 for (auto i = 0*m; i < m; ++i, ++row) {
230 const auto n = this->ia_[i + 1] - this->ia_[i + 0];
231
232 rowIdx.insert(rowIdx.end(), n, row);
233 }
234
235 return rowIdx;
236}
237
238template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
239void
242{
243 this->ia_.clear();
244 this->ja_.clear();
245
246 if constexpr (TrackCompressedIdx) {
247 this->compressedIdx_.clear();
248 }
249
250 this->numRows_ = 0;
251 this->numCols_ = 0;
252}
253
254template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
255void
257CSR::assemble(const Neighbours& rows,
258 const Neighbours& cols,
259 const BaseVertexID maxRowIdx,
260 const BaseVertexID maxColIdx,
261 [[maybe_unused]] const bool expandExistingIdxMap)
262{
263 [[maybe_unused]] auto compressedIdx = this->compressedIdx_;
264 [[maybe_unused]] const auto numOrigNNZ = this->ja_.size();
265
266 auto i = this->coordinateFormatRowIndices();
267 i.insert(i.end(), rows.begin(), rows.end());
268
269 auto j = this->ja_;
270 j.insert(j.end(), cols.begin(), cols.end());
271
272 const auto thisNumRows = std::max(this->numRows_, maxRowIdx + 1);
273 const auto thisNumCols = std::max(this->numCols_, maxColIdx + 1);
274
275 this->preparePushbackRowGrouping(thisNumRows, i);
276
277 this->groupAndTrackColumnIndicesByRow(i, j);
278
279 if constexpr (TrackCompressedIdx) {
280 if (expandExistingIdxMap) {
281 this->remapCompressedIndex(std::move(compressedIdx), numOrigNNZ);
282 }
283 }
284
285 this->numRows_ = thisNumRows;
286 this->numCols_ = thisNumCols;
287}
288
289template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
290void
292CSR::compress(const Offset maxNumVertices)
293{
294 if (this->numRows() > maxNumVertices) {
295 throw std::invalid_argument {
296 "Number of vertices in input graph (" +
297 std::to_string(this->numRows()) + ") "
298 "exceeds maximum graph size implied by explicit size of "
299 "adjacency matrix (" + std::to_string(maxNumVertices) + ')'
300 };
301 }
302
303 this->sortColumnIndicesPerRow();
304
305 // Must be called *after* sortColumnIndicesPerRow().
306 this->condenseDuplicates();
307
308 const auto nRows = this->startPointers().size() - 1;
309 if (nRows < maxNumVertices) {
310 this->ia_.insert(this->ia_.end(),
311 maxNumVertices - nRows,
312 this->startPointers().back());
313 }
314}
315
316template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
317void
320{
321 // Transposition is, in this context, effectively a linear time (O(nnz))
322 // bucket insertion procedure. In other words transposing the structure
323 // twice creates a structure with column indices in (ascendingly) sorted
324 // order.
325
326 this->transpose();
327 this->transpose();
328}
329
330template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
331void
334{
335 // Note: Must be called *after* sortColumnIndicesPerRow().
336
337 const auto colIdx = this->ja_;
338 auto end = colIdx.begin();
339
340 this->ja_.clear();
341
342 [[maybe_unused]] auto compressedIdx = this->compressedIdx_;
343 if constexpr (TrackCompressedIdx) {
344 this->compressedIdx_.clear();
345 }
346
347 const auto numRows = this->ia_.size() - 1;
348 for (auto row = 0*numRows; row < numRows; ++row) {
349 auto begin = end;
350
351 std::advance(end, this->ia_[row + 1] - this->ia_[row + 0]);
352
353 const auto q = this->ja_.size();
354
355 this->condenseAndTrackUniqueColumnsForSingleRow(begin, end);
356
357 this->ia_[row + 0] = q;
358 }
359
360 if constexpr (TrackCompressedIdx) {
361 this->remapCompressedIndex(std::move(compressedIdx));
362 }
363
364 // Record final table sizes.
365 this->ia_.back() = this->ja_.size();
366}
367
368template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
369void
371CSR::preparePushbackRowGrouping(const int numRows,
372 const Neighbours& rowIdx)
373{
374 assert (numRows >= 0);
375
376 this->ia_.assign(numRows + 1, 0);
377
378 // Count number of neighbouring vertices for each row. Accumulate in
379 // "next" bin since we're positioning the end pointers.
380 for (const auto& row : rowIdx) {
381 this->ia_[row + 1] += 1;
382 }
383
384 // Position "end" pointers.
385 //
386 // After this loop, ia_[i + 1] points to the *start* of the range of the
387 // column indices/neighbouring vertices of vertex 'i'. This, in turn,
388 // enables using the statement ja_[ia_[i+1]++] = v in groupAndTrack()
389 // to insert vertex 'v' as a neighbour, at the end of the range of known
390 // neighbours, *and* advance the end pointer of row/vertex 'i'. We use
391 // ia_[0] as an accumulator for the total number of neighbouring
392 // vertices in the graph.
393 //
394 // Note index range: 1..numRows inclusive.
395 for (typename Start::size_type i = 1, n = numRows; i <= n; ++i) {
396 this->ia_[0] += this->ia_[i];
397 this->ia_[i] = this->ia_[0] - this->ia_[i];
398 }
399
400 assert (this->ia_[0] == rowIdx.size());
401}
402
403template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
404void
406CSR::groupAndTrackColumnIndicesByRow(const Neighbours& rowIdx,
407 const Neighbours& colIdx)
408{
409 assert (this->ia_[0] == rowIdx.size());
410
411 const auto nnz = rowIdx.size();
412
413 this->ja_.resize(nnz);
414
415 if constexpr (TrackCompressedIdx) {
416 this->compressedIdx_.clear();
417 this->compressedIdx_.reserve(nnz);
418 }
419
420 // Group/insert column indices according to their associate vertex/row
421 // index.
422 //
423 // At the start of the loop the end pointers ia_[i+1], formed in
424 // preparePushback(), are positioned at the *start* of the column index
425 // range associated to vertex 'i'. After this loop all vertices
426 // neighbouring vertex 'i' will be placed consecutively, in order of
427 // appearance, into ja_. Furthermore, the row pointers ia_ will have
428 // their final position.
429 //
430 // The statement ja_[ia_[i+1]++] = v, split into two statements using
431 // the helper object 'k', inserts 'v' as a neighbouring vertex of vertex
432 // 'i' *and* advances the end pointer ia_[i+1] of that vertex. We use
433 // and maintain the invariant that ia_[i+1] at all times records the
434 // insertion point of the next neighbouring vertex of vertex 'i'. When
435 // the list of neighbouring vertices for vertex 'i' has been exhausted,
436 // ia_[i+1] will hold the start position for in ja_ for vertex i+1.
437 for (auto nz = 0*nnz; nz < nnz; ++nz) {
438 const auto k = this->ia_[rowIdx[nz] + 1] ++;
439
440 this->ja_[k] = colIdx[nz];
441
442 if constexpr (TrackCompressedIdx) {
443 this->compressedIdx_.push_back(k);
444 }
445 }
446
447 this->ia_[0] = 0;
448}
449
450template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
451void
454{
455 [[maybe_unused]] auto compressedIdx = this->compressedIdx_;
456
457 {
458 const auto rowIdx = this->coordinateFormatRowIndices();
459 const auto colIdx = this->ja_;
460
461 this->preparePushbackRowGrouping(this->numCols_, colIdx);
462
463 // Note parameter order. Transposition switches role of rows and
464 // columns.
465 this->groupAndTrackColumnIndicesByRow(colIdx, rowIdx);
466 }
467
468 if constexpr (TrackCompressedIdx) {
469 this->remapCompressedIndex(std::move(compressedIdx));
470 }
471
472 std::swap(this->numRows_, this->numCols_);
473}
474
475template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
476void
478CSR::condenseAndTrackUniqueColumnsForSingleRow(typename Neighbours::const_iterator begin,
479 typename Neighbours::const_iterator end)
480{
481 // We assume that we're only called *after* sortColumnIndicesPerRow()
482 // whence duplicate elements appear consecutively in [begin, end).
483 //
484 // Note: This is essentially the same as std::unique(begin, end) save
485 // for the return value and the fact that we additionally record the
486 // 'compressedIdx_' mapping. That mapping enables subsequent, decoupled
487 // accumulation of the 'sa_' contributions.
488
489 while (begin != end) {
490 // Note: Order of ja_ and compressedIdx_ matters here.
491
492 if constexpr (TrackCompressedIdx) {
493 this->compressedIdx_.push_back(this->ja_.size());
494 }
495
496 this->ja_.push_back(*begin);
497
498 auto next_unique =
499 std::find_if(begin, end, [last = this->ja_.back()]
500 (const auto j) { return j != last; });
501
502 if constexpr (TrackCompressedIdx) {
503 // Number of duplicate elements in [begin, next_unique).
504 const auto ndup = std::distance(begin, next_unique);
505
506 if (ndup > 1) {
507 // Insert ndup - 1 copies of .back() to represent the
508 // duplicate pairs in [begin, next_unique). We subtract one
509 // to account for .push_back() above representing *begin.
510 this->compressedIdx_.insert(this->compressedIdx_.end(),
511 ndup - 1,
512 this->compressedIdx_.back());
513 }
514 }
515
516 begin = next_unique;
517 }
518}
519
520template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
521void
523remapCompressedIndex([[maybe_unused]] Start&& compressedIdx,
524 [[maybe_unused]] std::optional<typename Start::size_type> numOrig)
525{
526 if constexpr (TrackCompressedIdx) {
527 for (auto& i : compressedIdx) {
528 i = this->compressedIdx_[i];
529 }
530
531 if (numOrig.has_value() && (*numOrig < this->compressedIdx_.size())) {
532 // Client called add() after compress(). Remap existing portion
533 // of compressedIdx (above), and append new entries (here).
534 compressedIdx
535 .insert(compressedIdx.end(),
536 this->compressedIdx_.begin() + *numOrig,
537 this->compressedIdx_.end());
538 }
539
540 this->compressedIdx_.swap(compressedIdx);
541 }
542}
543
544// =====================================================================
545
546// ---------------------------------------------------------------------
547// Class Opm::utility::CSRGraphFromCoordinates
548// ---------------------------------------------------------------------
549
550template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
556
557template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
558void
560addConnection(const VertexID v1, const VertexID v2)
561{
562 if ((v1 < 0) || (v2 < 0)) {
563 throw std::invalid_argument {
564 "Vertex IDs must be non-negative. Got (v1,v2) = ("
565 + std::to_string(v1) + ", " + std::to_string(v2)
566 + ')'
567 };
568 }
569
570 if constexpr (! PermitSelfConnections) {
571 if (v1 == v2) {
572 // Ignore self connections.
573 return;
574 }
575 }
576
577 this->uncompressed_.add(v1, v2);
578}
579
580template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
581void
583compress(const Offset maxNumVertices, const bool expandExistingIdxMap)
584{
585 if (! this->uncompressed_.isValid()) {
586 throw std::logic_error {
587 "Cannot compress invalid connection list"
588 };
589 }
590
591 this->csr_.merge(this->uncompressed_, maxNumVertices, expandExistingIdxMap);
592
593 this->uncompressed_.clear();
594}
595
596template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
602
603template <typename VertexID, bool TrackCompressedIdx, bool PermitSelfConnections>
606{
607 const auto& ia = this->startPointers();
608
609 return ia.empty() ? 0 : ia.back();
610}
Form CSR adjacency matrix representation of unstructured graphs.
Definition CSRGraphFromCoordinates.hpp:56
Offset numEdges() const
Retrieve number of edges (non-zero matrix elements) in input graph.
Definition CSRGraphFromCoordinates_impl.hpp:605
void clear()
Clear all internal buffers, but preserve allocated capacity.
Definition CSRGraphFromCoordinates_impl.hpp:551
void addConnection(VertexID v1, VertexID v2)
Add flow rate connection between regions.
Definition CSRGraphFromCoordinates_impl.hpp:560
typename Neighbours::size_type Offset
Offset into neighbour array.
Definition CSRGraphFromCoordinates.hpp:68
Offset numVertices() const
Retrieve number of rows (source entities) in input graph.
Definition CSRGraphFromCoordinates_impl.hpp:598
void compress(Offset maxNumVertices, bool expandExistingIdxMap=false)
Form CSR adjacency matrix representation of input graph from connections established in previous call...
Definition CSRGraphFromCoordinates_impl.hpp:583
std::vector< BaseVertexID > Neighbours
Representation of neighbouring regions.
Definition CSRGraphFromCoordinates.hpp:65
std::vector< Offset > Start
CSR start pointers.
Definition CSRGraphFromCoordinates.hpp:71