Extract Facility to Form CSR Representations of Graphs
This commit extracts the internal helpers of class
InterRegFlowMap
out to a new public helper class template, CSRGraphFromCoordinates.
Client code can, at the expense of one additional data member and
some dynamic memory, elect to track the index pairs. This enables
O(1) assembly per element when used as part of a CSR matrix with a
value array, SA.
Class CSRGraphFromCoordinates does not track values. It is purely
for the sake of forming the IA and JA structure arrays. Upon
calling 'compress()', column indices are sorted per row and
duplicate column indices condensed to a single, unique,
representative.
This commit is contained in:
Bård Skaflestad
@@ -463,6 +463,7 @@ if(ENABLE_ECL_INPUT)
|
||||
tests/rst_test.cpp
|
||||
tests/test_ActiveGridCells.cpp
|
||||
tests/test_CopyablePtr.cpp
|
||||
tests/test_CSRGraphFromCoordinates.cpp
|
||||
tests/test_ERsm.cpp
|
||||
tests/test_GuideRate.cpp
|
||||
tests/test_RestartFileView.cpp
|
||||
@@ -756,6 +757,8 @@ list( APPEND PUBLIC_HEADER_FILES
|
||||
opm/common/OpmLog/StreamLog.hpp
|
||||
opm/common/OpmLog/TimerLog.hpp
|
||||
opm/common/utility/ActiveGridCells.hpp
|
||||
opm/common/utility/CSRGraphFromCoordinates.hpp
|
||||
opm/common/utility/CSRGraphFromCoordinates_impl.hpp
|
||||
opm/common/utility/Demangle.hpp
|
||||
opm/common/utility/FileSystem.hpp
|
||||
opm/common/utility/OpmInputError.hpp
|
||||
|
||||
526
opm/common/utility/CSRGraphFromCoordinates.hpp
Normal file
526
opm/common/utility/CSRGraphFromCoordinates.hpp
Normal file
@@ -0,0 +1,526 @@
|
||||
/*
|
||||
Copyright 2016 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2016 Statoil ASA.
|
||||
Copyright 2022 Equinor ASA
|
||||
|
||||
This file is part of the Open Porous Media Project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_UTILITY_CSRGRAPHFROMCOORDINATES_HPP
|
||||
#define OPM_UTILITY_CSRGRAPHFROMCOORDINATES_HPP
|
||||
|
||||
#include <cstddef>
|
||||
#include <optional>
|
||||
#include <type_traits>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
/// \file
|
||||
///
|
||||
/// Facility for converting collection of region ID pairs into a sparse
|
||||
/// (CSR) adjacency matrix representation of a graph. Supports O(nnz)
|
||||
/// compression.
|
||||
|
||||
namespace Opm { namespace utility {
|
||||
|
||||
/// Form CSR adjacency matrix representation of unstructured graphs.
|
||||
/// Optionally maps vertex pairs to compressed indices to enable O(1)
|
||||
/// per-element lookup in assembly-like operations.
|
||||
///
|
||||
/// \tparam VertexID ID type of an abstract vertex in the graph. Could
|
||||
/// for instance be used to represent a cell index or a region index.
|
||||
/// Must be an integral type.
|
||||
///
|
||||
/// \tparam TrackCompressedIdx Whether or not to form a mapping relation
|
||||
/// for vertex pairs to compressed indices. Default value, false,
|
||||
/// bypasses this mapping relation and conserves memory.
|
||||
template <typename VertexID = int, bool TrackCompressedIdx = false>
|
||||
class CSRGraphFromCoordinates
|
||||
{
|
||||
private:
|
||||
using BaseVertexID = std::remove_cv_t<VertexID>;
|
||||
|
||||
static_assert(std::is_integral_v<BaseVertexID>,
|
||||
"The VertexID must be an integral type");
|
||||
|
||||
public:
|
||||
/// Representation of neighbouring regions.
|
||||
using Neighbours = std::vector<BaseVertexID>;
|
||||
|
||||
/// Offset into neighbour array.
|
||||
using Offset = typename Neighbours::size_type;
|
||||
|
||||
/// CSR start pointers.
|
||||
using Start = std::vector<Offset>;
|
||||
|
||||
/// Clear all internal buffers, but preserve allocated capacity.
|
||||
void clear();
|
||||
|
||||
/// Add flow rate connection between regions.
|
||||
///
|
||||
/// \param[in] v1 First vertex in vertex pair. Used as row index.
|
||||
///
|
||||
/// \param[in] r2 Second vertex in vertex pair. Used as column index.
|
||||
///
|
||||
/// If both vertex IDs are the same then this function does nothing.
|
||||
void addConnection(VertexID v1, VertexID v2);
|
||||
|
||||
/// Form CSR adjacency matrix representation of input graph from
|
||||
/// connections established in previous calls to addConnection().
|
||||
///
|
||||
/// \param[in] maxNumVertices Number of rows in resulting CSR
|
||||
/// matrix. If prior calls to addConnection() supply vertex IDs
|
||||
/// (row indices) greater than or equal to \p maxNumVertices,
|
||||
/// then method compress() will throw \code
|
||||
/// std::invalid_argument \endcode.
|
||||
///
|
||||
/// \param[in] expandExistingIdxMap Whether or not preserve and
|
||||
/// update the existing compressed index map. This is
|
||||
/// potentially useful for the case of adding new connections to
|
||||
/// an already compressed graph. The default setting, \c false,
|
||||
/// will disregard any existing index map and create the index
|
||||
/// map from scratch. This runtime parameter is unused if the
|
||||
/// class is not configured to track compressed indices through
|
||||
/// the TrackCompressedIdx class parameter.
|
||||
void compress(Offset maxNumVertices,
|
||||
bool expandExistingIdxMap = false);
|
||||
|
||||
/// Retrieve number of rows (source entities) in input graph.
|
||||
/// Corresponds to value of argument passed to compress(). Valid
|
||||
/// only after calling compress().
|
||||
Offset numVertices() const;
|
||||
|
||||
/// Retrieve number of edges (non-zero matrix elements) in input
|
||||
/// graph.
|
||||
Offset numEdges() const;
|
||||
|
||||
/// Read-only access to compressed structure's start pointers.
|
||||
const Start& startPointers() const
|
||||
{
|
||||
return this->csr_.startPointers();
|
||||
}
|
||||
|
||||
/// Read-only access to compressed structure's column indices,
|
||||
/// ascendingly sorted per row.
|
||||
const Neighbours& columnIndices() const
|
||||
{
|
||||
return this->csr_.columnIndices();
|
||||
}
|
||||
|
||||
/// Read-only access to mapping from input order vertex pairs to
|
||||
/// compressed structure's edge index (location in ja/sa).
|
||||
///
|
||||
/// Available only if client code sets TrackCompressedIdx=true.
|
||||
/// Compiler diagnostic, typically referring to 'enable_if', if
|
||||
/// client code tries to call this function without setting
|
||||
/// class parameter TrackCompressedIdx=true.
|
||||
template <typename Ret = const Start&>
|
||||
std::enable_if_t<TrackCompressedIdx, Ret> compressedIndexMap() const
|
||||
{
|
||||
return this->csr_.compressedIndexMap();
|
||||
}
|
||||
|
||||
// MessageBufferType API should be similar to Dune::MessageBufferIF
|
||||
template <class MessageBufferType>
|
||||
void write(MessageBufferType& buffer) const
|
||||
{
|
||||
this->csr_.write(buffer);
|
||||
}
|
||||
|
||||
// MessageBufferType API should be similar to Dune::MessageBufferIF
|
||||
template <class MessageBufferType>
|
||||
void read(MessageBufferType& buffer)
|
||||
{
|
||||
auto other = CSR{};
|
||||
other.read(buffer);
|
||||
|
||||
this->uncompressed_
|
||||
.add(other.maxRowID(),
|
||||
other.maxColID(),
|
||||
other.coordinateFormatRowIndices(),
|
||||
other.columnIndices());
|
||||
}
|
||||
|
||||
private:
|
||||
/// Coordinate format representation of individual contributions to
|
||||
/// inter-region flows.
|
||||
class Connections
|
||||
{
|
||||
public:
|
||||
/// Add contributions from a single inter-region connection.
|
||||
///
|
||||
/// \param[in] r1 Source region. Zero-based region index/ID.
|
||||
///
|
||||
/// \param[in] r2 Destination region. Zero-based region index.
|
||||
void add(VertexID v1, VertexID v2);
|
||||
|
||||
/// Add contributions from multiple inter-region connections.
|
||||
///
|
||||
/// \param[in] maxRowIdx Maximum row (source region) index
|
||||
/// across all new inter-region connection contributions.
|
||||
///
|
||||
/// \param[in] maxColIdx Maximum column (destination region)
|
||||
/// index across all new inter-region contributions.
|
||||
///
|
||||
/// \param[in] rows Source region indices for all new
|
||||
/// inter-region connection contributions.
|
||||
///
|
||||
/// \param[in] cols Destination region indices for all new
|
||||
/// inter-region connection contributions.
|
||||
void add(VertexID maxRowIdx,
|
||||
VertexID maxColIdx,
|
||||
const Neighbours& rows,
|
||||
const Neighbours& cols);
|
||||
|
||||
/// Clear internal tables. Preserve allocated capacity.
|
||||
void clear();
|
||||
|
||||
/// Predicate.
|
||||
///
|
||||
/// \return Whether or not internal tables are empty.
|
||||
bool empty() const;
|
||||
|
||||
/// Whether or not internal tables meet size consistency
|
||||
/// requirements.
|
||||
bool isValid() const;
|
||||
|
||||
/// Maximum zero-based row index.
|
||||
std::optional<BaseVertexID> maxRow() const;
|
||||
|
||||
/// Maximum zero-based column index.
|
||||
std::optional<BaseVertexID> maxCol() const;
|
||||
|
||||
/// Number of uncompressed contributions in internal tables.
|
||||
typename Neighbours::size_type numContributions() const;
|
||||
|
||||
/// Read-only access to uncompressed row indices.
|
||||
const Neighbours& rowIndices() const;
|
||||
|
||||
/// Read-only access to uncompressed column indices.
|
||||
const Neighbours& columnIndices() const;
|
||||
|
||||
private:
|
||||
/// Zero-based row/source region indices.
|
||||
Neighbours i_{};
|
||||
|
||||
/// Zero-based column/destination region indices.
|
||||
Neighbours j_{};
|
||||
|
||||
/// Maximum row index in \code this->i_ \endcode.
|
||||
std::optional<VertexID> max_i_{};
|
||||
|
||||
/// Maximum column index in \code this->j_ \endcode.
|
||||
std::optional<VertexID> max_j_{};
|
||||
};
|
||||
|
||||
/// Compressed sparse row representation of inter-region flow rates
|
||||
///
|
||||
/// Row and column indices are zero-based vertex IDs. Column
|
||||
/// indices ascendingly sorted per row.
|
||||
class CSR
|
||||
{
|
||||
public:
|
||||
/// Merge coordinate format into existing CSR map.
|
||||
///
|
||||
/// \param[in] conns Coordinate representation of new
|
||||
/// contributions.
|
||||
///
|
||||
/// \param[in] maxNumVertices Maximum number of vertices.
|
||||
///
|
||||
/// \param[in] expandExistingIdxMap Whether or not preserve and
|
||||
/// update the existing compressed index map. This is
|
||||
/// potentially useful for the case of adding new connections
|
||||
/// to an already compressed graph. The default setting, \c
|
||||
/// false, will disregard any existing index map and create
|
||||
/// the index map from scratch. This runtime parameter is
|
||||
/// unused if the class is not configured to track compressed
|
||||
/// indices through the TrackCompressedIdx class parameter.
|
||||
void merge(const Connections& conns,
|
||||
const Offset maxNumVertices,
|
||||
const bool expandExistingIdxMap);
|
||||
|
||||
/// Total number of rows in compressed map structure.
|
||||
Offset numRows() const;
|
||||
|
||||
/// Maximum zero-based row index encountered in mapped structure.
|
||||
BaseVertexID maxRowID() const;
|
||||
|
||||
/// Maximum zero-based column index encountered in mapped structure.
|
||||
BaseVertexID maxColID() const;
|
||||
|
||||
/// Read-only access to compressed structure's start pointers.
|
||||
const Start& startPointers() const;
|
||||
|
||||
/// Read-only access to compressed structure's column indices,
|
||||
/// ascendingly sorted per rwo.
|
||||
const Neighbours& columnIndices() const;
|
||||
|
||||
/// Coordinate format row index vector. Expanded from \code
|
||||
/// startPointers() \endcode.
|
||||
Neighbours coordinateFormatRowIndices() const;
|
||||
|
||||
template <typename Ret = const Start&>
|
||||
std::enable_if_t<TrackCompressedIdx, Ret> compressedIndexMap() const
|
||||
{
|
||||
return this->compressedIdx_;
|
||||
}
|
||||
|
||||
// MessageBufferType API should be similar to Dune::MessageBufferIF
|
||||
template <class MessageBufferType>
|
||||
void write(MessageBufferType& buffer) const
|
||||
{
|
||||
this->writeVector(this->ia_, buffer);
|
||||
this->writeVector(this->ja_, buffer);
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->writeVector(this->compressedIdx_, buffer);
|
||||
}
|
||||
|
||||
buffer.write(this->numRows_);
|
||||
buffer.write(this->numCols_);
|
||||
}
|
||||
|
||||
// MessageBufferType API should be similar to Dune::MessageBufferIF
|
||||
template <class MessageBufferType>
|
||||
void read(MessageBufferType& buffer)
|
||||
{
|
||||
this->readVector(buffer, this->ia_);
|
||||
this->readVector(buffer, this->ja_);
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->readVector(buffer, this->compressedIdx_);
|
||||
}
|
||||
|
||||
buffer.read(this->numRows_);
|
||||
buffer.read(this->numCols_);
|
||||
}
|
||||
|
||||
/// Clear internal tables. Preserve allocated capacity.
|
||||
void clear();
|
||||
|
||||
private:
|
||||
struct EmptyPlaceHolder {};
|
||||
|
||||
/// Start pointers.
|
||||
Start ia_{};
|
||||
|
||||
/// Column indices. Ascendingly sorted per row once structure
|
||||
/// is fully established.
|
||||
Neighbours ja_{};
|
||||
|
||||
/// Destination index in compressed representation. Vector of
|
||||
/// size equal to number of \c addConnection() calls if client
|
||||
/// code requests that compressed indices be tracked (i.e., when
|
||||
/// parameter TrackCompressedIdx == true); Empty structure
|
||||
/// otherwise (default setting).
|
||||
std::conditional_t<TrackCompressedIdx, Start, EmptyPlaceHolder> compressedIdx_{};
|
||||
|
||||
/// Number of active rows in compressed map structure.
|
||||
BaseVertexID numRows_{};
|
||||
|
||||
/// Number of active columns in compressed map structure.
|
||||
/// Tracked as the maximum column index plus one.
|
||||
BaseVertexID numCols_{};
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// Implementation of read()/write()
|
||||
// ---------------------------------------------------------
|
||||
|
||||
template <typename T, class A, class MessageBufferType>
|
||||
void writeVector(const std::vector<T,A>& vec,
|
||||
MessageBufferType& buffer) const
|
||||
{
|
||||
const auto n = vec.size();
|
||||
buffer.write(n);
|
||||
|
||||
for (const auto& x : vec) {
|
||||
buffer.write(x);
|
||||
}
|
||||
}
|
||||
|
||||
template <class MessageBufferType, typename T, class A>
|
||||
void readVector(MessageBufferType& buffer,
|
||||
std::vector<T,A>& vec)
|
||||
{
|
||||
auto n = 0 * vec.size();
|
||||
buffer.read(n);
|
||||
|
||||
vec.resize(n);
|
||||
|
||||
for (auto& x : vec) {
|
||||
buffer.read(x);
|
||||
}
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// Implementation of merge()
|
||||
// ---------------------------------------------------------
|
||||
|
||||
/// Incorporate new, coordinate format contributions into
|
||||
/// existing, possibly empty, CSR mapping structure.
|
||||
///
|
||||
/// On exit the ia_ array holds the proper start pointers while
|
||||
/// ja_ holds the corresponding column indices albeit possibly
|
||||
/// repeated and unsorted.
|
||||
///
|
||||
/// \param[in] rows Row indices of all, possibly repeated,
|
||||
/// coordinate format input contributions. Start pointers \c
|
||||
/// ia_ updated to account for new entries.
|
||||
///
|
||||
/// \param[in] cols Column index of coordinate format intput
|
||||
/// structure. Inserted into \c ja_ according to its
|
||||
/// corresponding row index.
|
||||
///
|
||||
/// \param[in] maxRowID Maximum index in \p rows. Needed to
|
||||
/// ensure proper size of \c ia_.
|
||||
///
|
||||
/// \param[in] maxColID Maximum index in \p cols.
|
||||
///
|
||||
/// \param[in] expandExistingIdxMap Whether or not preserve and
|
||||
/// update the existing compressed index map. This is
|
||||
/// potentially useful for the case of adding new connections
|
||||
/// to an already compressed graph. The default setting, \c
|
||||
/// false, will disregard any existing index map and create
|
||||
/// the index map from scratch. This runtime parameter is
|
||||
/// unused if the class is not configured to track compressed
|
||||
/// indices through the TrackCompressedIdx class parameter.
|
||||
void assemble(const Neighbours& rows,
|
||||
const Neighbours& cols,
|
||||
BaseVertexID maxRowID,
|
||||
BaseVertexID maxColID,
|
||||
bool expandExistingIdxMap);
|
||||
|
||||
/// Sort column indices per row and compress repeated column
|
||||
/// indices down to a single unique element per row. Sum
|
||||
/// repeated values
|
||||
///
|
||||
/// On exit the \c ia_ and \c ja_ arrays all have their
|
||||
/// expected, canonical structure.
|
||||
///
|
||||
/// \param[in] maxNumVertices Maximum number of vertices
|
||||
/// supported by final compressed mapping structure. Ignored
|
||||
/// if less than active number of rows.
|
||||
void compress(const Offset maxNumVertices);
|
||||
|
||||
/// Sort column indices within each mapped row.
|
||||
///
|
||||
/// On exit \c ja_ has ascendingly sorted column indices, albeit
|
||||
/// possibly with repeated entries. This function also updates
|
||||
/// \c compressedIdx_, if applicable, to account for the new
|
||||
/// locations of the non-zero elements in the grouped structure.
|
||||
void sortColumnIndicesPerRow();
|
||||
|
||||
/// Condense repeated column indices per row down to a single
|
||||
/// unique entry for each.
|
||||
///
|
||||
/// Assumes that each row has ascendingly sorted column indices
|
||||
/// in \c ja_ and must therefore be called after member function
|
||||
/// sortColumnIndicesPerRow(). On exit, \c ja_ has its final
|
||||
/// canonical structure and \c compressedIdx_, if applicable,
|
||||
/// knows the final location of each non-zero contribution in
|
||||
/// the input coordinate format.
|
||||
void condenseDuplicates();
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// Implementation of assemble()
|
||||
// ---------------------------------------------------------
|
||||
|
||||
/// Position end pointers at start of row to prepare for column
|
||||
/// index grouping by corresponding row index.
|
||||
///
|
||||
/// Also counts total number of non-zero elements, possibly
|
||||
/// including repetitions, in \code this->ia_[0] \endcode.
|
||||
///
|
||||
/// \param[in] numRows Number of rows in final compressed
|
||||
/// structure. Used to allocate \code this->ia_ \endcode.
|
||||
///
|
||||
/// \param[in] rowIdx Row indices of all, possibly repeated,
|
||||
/// coordinate format input contributions. Needed to count
|
||||
/// the number of possibly repeated column index entries per
|
||||
/// row.
|
||||
void preparePushbackRowGrouping(const int numRows,
|
||||
const Neighbours& rowIdx);
|
||||
|
||||
/// Group column indices by corresponding row index and track
|
||||
/// grouped location of original coordinate format element
|
||||
///
|
||||
/// Appends grouped location to \c compressedIdx_ if needed.
|
||||
///
|
||||
/// \param[in] rowIdx Row index of coordinate format input
|
||||
/// structure. Used as grouping key.
|
||||
///
|
||||
/// \param[in] colIdx Column index of coordinate format intput
|
||||
/// structure. Inserted into \c ja_ according to its
|
||||
/// corresponding row index.
|
||||
void groupAndTrackColumnIndicesByRow(const Neighbours& rowIdx,
|
||||
const Neighbours& colIdx);
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// General utilities
|
||||
// ---------------------------------------------------------
|
||||
|
||||
/// Transpose connectivity structure.
|
||||
///
|
||||
/// Essentially swaps the roles of rows and columns. Also used
|
||||
/// as a basic building block for sortColumnIndicesPerRow().
|
||||
void transpose();
|
||||
|
||||
/// Condense sequences of repeated column indices in a single
|
||||
/// map row down to a single copy of each unique column index.
|
||||
///
|
||||
/// Appends new unique column indices to \code ja_ \endcode
|
||||
///
|
||||
/// Assumes that the map row has ascendingly sorted column
|
||||
/// indices and therefore has the same requirements as
|
||||
/// std::unique. Will also update the internal compressedIdx_
|
||||
/// mapping, if needed, to record new compressed locations for
|
||||
/// the current, uncompressed, non-zero map elements.
|
||||
///
|
||||
/// \param[in] begin Start of map row that contains possibly
|
||||
/// repeated column indices.
|
||||
///
|
||||
/// \param[in] end One-past-end of map row that contains
|
||||
/// possibly repeated column indices.
|
||||
void condenseAndTrackUniqueColumnsForSingleRow(typename Neighbours::const_iterator begin,
|
||||
typename Neighbours::const_iterator end);
|
||||
|
||||
/// Update \c compressedIdx_ mapping, if needed, to account for
|
||||
/// column index reshuffling.
|
||||
///
|
||||
/// \param[in] compressedIdx New compressed index locations of
|
||||
/// the non-zero map entries.
|
||||
///
|
||||
/// \param[in] numOrigNNZ Number of existing, unique NNZs
|
||||
/// (edges) in graph. Needed to support calling add() after
|
||||
/// compress() when TrackCompressedIdx is true.
|
||||
void remapCompressedIndex(Start&& compressedIdx,
|
||||
std::optional<typename Start::size_type> numOrigNNZ = std::nullopt);
|
||||
};
|
||||
|
||||
/// Accumulated coordinate format contributions that have not yet
|
||||
/// been added to the final CSR structure.
|
||||
Connections uncompressed_;
|
||||
|
||||
/// Canonical representation of unique inter-region flow rates.
|
||||
CSR csr_;
|
||||
};
|
||||
|
||||
}} // namespace Opm::utility
|
||||
|
||||
// Actual implementation of member functions in _impl.hpp file.
|
||||
#include <opm/common/utility/CSRGraphFromCoordinates_impl.hpp>
|
||||
|
||||
#endif // OPM_UTILITY_CSRGRAPHFROMCOORDINATES_HPP
|
||||
608
opm/common/utility/CSRGraphFromCoordinates_impl.hpp
Normal file
608
opm/common/utility/CSRGraphFromCoordinates_impl.hpp
Normal file
@@ -0,0 +1,608 @@
|
||||
/*
|
||||
Copyright 2016 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2016 Statoil ASA.
|
||||
Copyright 2022 Equinor ASA
|
||||
|
||||
This file is part of the Open Porous Media Project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <algorithm>
|
||||
#include <cassert>
|
||||
#include <exception>
|
||||
#include <iterator>
|
||||
#include <optional>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
// ---------------------------------------------------------------------
|
||||
// Class Opm::utility::CSRGraphFromCoordinates::Connections
|
||||
// ---------------------------------------------------------------------
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::add(const VertexID v1, const VertexID v2)
|
||||
{
|
||||
this->i_.push_back(v1);
|
||||
this->j_.push_back(v2);
|
||||
|
||||
this->max_i_ = std::max(this->max_i_.value_or(BaseVertexID{}), this->i_.back());
|
||||
this->max_j_ = std::max(this->max_j_.value_or(BaseVertexID{}), this->j_.back());
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::add(VertexID maxRowIdx,
|
||||
VertexID maxColIdx,
|
||||
const Neighbours& rows,
|
||||
const Neighbours& cols)
|
||||
{
|
||||
if (cols.size() != rows.size()) {
|
||||
throw std::invalid_argument {
|
||||
"Coordinate format column index table size does not match "
|
||||
"row index table size"
|
||||
};
|
||||
}
|
||||
|
||||
this->i_.insert(this->i_.end(), rows .begin(), rows .end());
|
||||
this->j_.insert(this->j_.end(), cols .begin(), cols .end());
|
||||
|
||||
this->max_i_ = std::max(this->max_i_.value_or(BaseVertexID{}), maxRowIdx);
|
||||
this->max_j_ = std::max(this->max_j_.value_or(BaseVertexID{}), maxColIdx);
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::clear()
|
||||
{
|
||||
this->j_.clear();
|
||||
this->i_.clear();
|
||||
|
||||
this->max_i_.reset();
|
||||
this->max_j_.reset();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
bool
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::empty() const
|
||||
{
|
||||
return this->i_.empty();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
bool
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::isValid() const
|
||||
{
|
||||
return this->i_.size() == this->j_.size();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
std::optional<typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::BaseVertexID>
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::maxRow() const
|
||||
{
|
||||
return this->max_i_;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
std::optional<typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::BaseVertexID>
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::maxCol() const
|
||||
{
|
||||
return this->max_j_;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Neighbours::size_type
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::numContributions() const
|
||||
{
|
||||
return this->i_.size();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
const typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Neighbours&
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::rowIndices() const
|
||||
{
|
||||
return this->i_;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
const typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Neighbours&
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
Connections::columnIndices() const
|
||||
{
|
||||
return this->j_;
|
||||
}
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// ---------------------------------------------------------------------
|
||||
// Class Opm::utility::CSRGraphFromCoordinates::CSR
|
||||
// ---------------------------------------------------------------------
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::merge(const Connections& conns,
|
||||
const Offset maxNumVertices,
|
||||
const bool expandExistingIdxMap)
|
||||
{
|
||||
const auto maxRow = conns.maxRow();
|
||||
|
||||
if (maxRow.has_value() &&
|
||||
(static_cast<Offset>(*maxRow) >= maxNumVertices))
|
||||
{
|
||||
throw std::invalid_argument {
|
||||
"Number of vertices in input graph (" +
|
||||
std::to_string(*maxRow) + ") "
|
||||
"exceeds maximum graph size implied by explicit size of "
|
||||
"adjacency matrix (" + std::to_string(maxNumVertices) + ')'
|
||||
};
|
||||
}
|
||||
|
||||
this->assemble(conns.rowIndices(), conns.columnIndices(),
|
||||
maxRow.value_or(BaseVertexID{0}),
|
||||
conns.maxCol().value_or(BaseVertexID{0}),
|
||||
expandExistingIdxMap);
|
||||
|
||||
this->compress(maxNumVertices);
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Offset
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::numRows() const
|
||||
{
|
||||
return this->startPointers().empty()
|
||||
? 0 : this->startPointers().size() - 1;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::BaseVertexID
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::maxRowID() const
|
||||
{
|
||||
return this->numRows_ - 1;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::BaseVertexID
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::maxColID() const
|
||||
{
|
||||
return this->numCols_ - 1;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
const typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Start&
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::startPointers() const
|
||||
{
|
||||
return this->ia_;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
const typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Neighbours&
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::columnIndices() const
|
||||
{
|
||||
return this->ja_;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Neighbours
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::coordinateFormatRowIndices() const
|
||||
{
|
||||
auto rowIdx = Neighbours{};
|
||||
|
||||
if (this->ia_.empty()) {
|
||||
return rowIdx;
|
||||
}
|
||||
|
||||
rowIdx.reserve(this->ia_.back());
|
||||
|
||||
auto row = BaseVertexID{};
|
||||
|
||||
const auto m = this->ia_.size() - 1;
|
||||
for (auto i = 0*m; i < m; ++i, ++row) {
|
||||
const auto n = this->ia_[i + 1] - this->ia_[i + 0];
|
||||
|
||||
rowIdx.insert(rowIdx.end(), n, row);
|
||||
}
|
||||
|
||||
return rowIdx;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::clear()
|
||||
{
|
||||
this->ia_.clear();
|
||||
this->ja_.clear();
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->compressedIdx_.clear();
|
||||
}
|
||||
|
||||
this->numRows_ = 0;
|
||||
this->numCols_ = 0;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::assemble(const Neighbours& rows,
|
||||
const Neighbours& cols,
|
||||
const BaseVertexID maxRowIdx,
|
||||
const BaseVertexID maxColIdx,
|
||||
[[maybe_unused]] const bool expandExistingIdxMap)
|
||||
{
|
||||
[[maybe_unused]] auto compressedIdx = this->compressedIdx_;
|
||||
[[maybe_unused]] const auto numOrigNNZ = this->ja_.size();
|
||||
|
||||
auto i = this->coordinateFormatRowIndices();
|
||||
i.insert(i.end(), rows.begin(), rows.end());
|
||||
|
||||
auto j = this->ja_;
|
||||
j.insert(j.end(), cols.begin(), cols.end());
|
||||
|
||||
const auto thisNumRows = std::max(this->numRows_, maxRowIdx + 1);
|
||||
const auto thisNumCols = std::max(this->numCols_, maxColIdx + 1);
|
||||
|
||||
this->preparePushbackRowGrouping(thisNumRows, i);
|
||||
|
||||
this->groupAndTrackColumnIndicesByRow(i, j);
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
if (expandExistingIdxMap) {
|
||||
this->remapCompressedIndex(std::move(compressedIdx), numOrigNNZ);
|
||||
}
|
||||
}
|
||||
|
||||
this->numRows_ = thisNumRows;
|
||||
this->numCols_ = thisNumCols;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::compress(const Offset maxNumVertices)
|
||||
{
|
||||
if (this->numRows() > maxNumVertices) {
|
||||
throw std::invalid_argument {
|
||||
"Number of vertices in input graph (" +
|
||||
std::to_string(this->numRows()) + ") "
|
||||
"exceeds maximum graph size implied by explicit size of "
|
||||
"adjacency matrix (" + std::to_string(maxNumVertices) + ')'
|
||||
};
|
||||
}
|
||||
|
||||
this->sortColumnIndicesPerRow();
|
||||
|
||||
// Must be called *after* sortColumnIndicesPerRow().
|
||||
this->condenseDuplicates();
|
||||
|
||||
const auto nRows = this->startPointers().size() - 1;
|
||||
if (nRows < maxNumVertices) {
|
||||
this->ia_.insert(this->ia_.end(),
|
||||
maxNumVertices - nRows,
|
||||
this->startPointers().back());
|
||||
}
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::sortColumnIndicesPerRow()
|
||||
{
|
||||
// Transposition is, in this context, effectively a linear time (O(nnz))
|
||||
// bucket insertion procedure. In other words transposing the structure
|
||||
// twice creates a structure with column indices in (ascendingly) sorted
|
||||
// order.
|
||||
|
||||
this->transpose();
|
||||
this->transpose();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::condenseDuplicates()
|
||||
{
|
||||
// Note: Must be called *after* sortColumnIndicesPerRow().
|
||||
|
||||
const auto colIdx = this->ja_;
|
||||
auto end = colIdx.begin();
|
||||
|
||||
this->ja_.clear();
|
||||
|
||||
[[maybe_unused]] auto compressedIdx = this->compressedIdx_;
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->compressedIdx_.clear();
|
||||
}
|
||||
|
||||
const auto numRows = this->ia_.size() - 1;
|
||||
for (auto row = 0*numRows; row < numRows; ++row) {
|
||||
auto begin = end;
|
||||
|
||||
std::advance(end, this->ia_[row + 1] - this->ia_[row + 0]);
|
||||
|
||||
const auto q = this->ja_.size();
|
||||
|
||||
this->condenseAndTrackUniqueColumnsForSingleRow(begin, end);
|
||||
|
||||
this->ia_[row + 0] = q;
|
||||
}
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->remapCompressedIndex(std::move(compressedIdx));
|
||||
}
|
||||
|
||||
// Record final table sizes.
|
||||
this->ia_.back() = this->ja_.size();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::preparePushbackRowGrouping(const int numRows,
|
||||
const Neighbours& rowIdx)
|
||||
{
|
||||
assert (numRows >= 0);
|
||||
|
||||
this->ia_.assign(numRows + 1, 0);
|
||||
|
||||
// Count number of neighbouring vertices for each row. Accumulate in
|
||||
// "next" bin since we're positioning the end pointers.
|
||||
for (const auto& row : rowIdx) {
|
||||
this->ia_[row + 1] += 1;
|
||||
}
|
||||
|
||||
// Position "end" pointers.
|
||||
//
|
||||
// After this loop, ia_[i + 1] points to the *start* of the range of the
|
||||
// column indices/neighbouring vertices of vertex 'i'. This, in turn,
|
||||
// enables using the statement ja_[ia_[i+1]++] = v in groupAndTrack()
|
||||
// to insert vertex 'v' as a neighbour, at the end of the range of known
|
||||
// neighbours, *and* advance the end pointer of row/vertex 'i'. We use
|
||||
// ia_[0] as an accumulator for the total number of neighbouring
|
||||
// vertices in the graph.
|
||||
//
|
||||
// Note index range: 1..numRows inclusive.
|
||||
for (typename Start::size_type i = 1, n = numRows; i <= n; ++i) {
|
||||
this->ia_[0] += this->ia_[i];
|
||||
this->ia_[i] = this->ia_[0] - this->ia_[i];
|
||||
}
|
||||
|
||||
assert (this->ia_[0] == rowIdx.size());
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::groupAndTrackColumnIndicesByRow(const Neighbours& rowIdx,
|
||||
const Neighbours& colIdx)
|
||||
{
|
||||
assert (this->ia_[0] == rowIdx.size());
|
||||
|
||||
const auto nnz = rowIdx.size();
|
||||
|
||||
this->ja_.resize(nnz);
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->compressedIdx_.clear();
|
||||
this->compressedIdx_.reserve(nnz);
|
||||
}
|
||||
|
||||
// Group/insert column indices according to their associate vertex/row
|
||||
// index.
|
||||
//
|
||||
// At the start of the loop the end pointers ia_[i+1], formed in
|
||||
// preparePushback(), are positioned at the *start* of the column index
|
||||
// range associated to vertex 'i'. After this loop all vertices
|
||||
// neighbouring vertex 'i' will be placed consecutively, in order of
|
||||
// appearance, into ja_. Furthermore, the row pointers ia_ will have
|
||||
// their final position.
|
||||
//
|
||||
// The statement ja_[ia_[i+1]++] = v, split into two statements using
|
||||
// the helper object 'k', inserts 'v' as a neighbouring vertex of vertex
|
||||
// 'i' *and* advances the end pointer ia_[i+1] of that vertex. We use
|
||||
// and maintain the invariant that ia_[i+1] at all times records the
|
||||
// insertion point of the next neighbouring vertex of vertex 'i'. When
|
||||
// the list of neighbouring vertices for vertex 'i' has been exhausted,
|
||||
// ia_[i+1] will hold the start position for in ja_ for vertex i+1.
|
||||
for (auto nz = 0*nnz; nz < nnz; ++nz) {
|
||||
const auto k = this->ia_[rowIdx[nz] + 1] ++;
|
||||
|
||||
this->ja_[k] = colIdx[nz];
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->compressedIdx_.push_back(k);
|
||||
}
|
||||
}
|
||||
|
||||
this->ia_[0] = 0;
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::transpose()
|
||||
{
|
||||
[[maybe_unused]] auto compressedIdx = this->compressedIdx_;
|
||||
|
||||
{
|
||||
const auto rowIdx = this->coordinateFormatRowIndices();
|
||||
const auto colIdx = this->ja_;
|
||||
|
||||
this->preparePushbackRowGrouping(this->numCols_, colIdx);
|
||||
|
||||
// Note parameter order. Transposition switches role of rows and
|
||||
// columns.
|
||||
this->groupAndTrackColumnIndicesByRow(colIdx, rowIdx);
|
||||
}
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->remapCompressedIndex(std::move(compressedIdx));
|
||||
}
|
||||
|
||||
std::swap(this->numRows_, this->numCols_);
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
CSR::condenseAndTrackUniqueColumnsForSingleRow(typename Neighbours::const_iterator begin,
|
||||
typename Neighbours::const_iterator end)
|
||||
{
|
||||
// We assume that we're only called *after* sortColumnIndicesPerRow()
|
||||
// whence duplicate elements appear consecutively in [begin, end).
|
||||
//
|
||||
// Note: This is essentially the same as std::unique(begin, end) save
|
||||
// for the return value and the fact that we additionally record the
|
||||
// 'compressedIdx_' mapping. That mapping enables subsequent, decoupled
|
||||
// accumulation of the 'sa_' contributions.
|
||||
|
||||
while (begin != end) {
|
||||
// Note: Order of ja_ and compressedIdx_ matters here.
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
this->compressedIdx_.push_back(this->ja_.size());
|
||||
}
|
||||
|
||||
this->ja_.push_back(*begin);
|
||||
|
||||
auto next_unique =
|
||||
std::find_if(begin, end, [last = this->ja_.back()]
|
||||
(const auto j) { return j != last; });
|
||||
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
// Number of duplicate elements in [begin, next_unique).
|
||||
const auto ndup = std::distance(begin, next_unique);
|
||||
|
||||
if (ndup > 1) {
|
||||
// Insert ndup - 1 copies of .back() to represent the
|
||||
// duplicate pairs in [begin, next_unique). We subtract one
|
||||
// to account for .push_back() above representing *begin.
|
||||
this->compressedIdx_.insert(this->compressedIdx_.end(),
|
||||
ndup - 1,
|
||||
this->compressedIdx_.back());
|
||||
}
|
||||
}
|
||||
|
||||
begin = next_unique;
|
||||
}
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::CSR::
|
||||
remapCompressedIndex([[maybe_unused]] Start&& compressedIdx,
|
||||
[[maybe_unused]] std::optional<typename Start::size_type> numOrig)
|
||||
{
|
||||
if constexpr (TrackCompressedIdx) {
|
||||
for (auto& i : compressedIdx) {
|
||||
i = this->compressedIdx_[i];
|
||||
}
|
||||
|
||||
if (numOrig.has_value() && (*numOrig < this->compressedIdx_.size())) {
|
||||
// Client called add() after compress(). Remap existing portion
|
||||
// of compressedIdx (above), and append new entries (here).
|
||||
compressedIdx
|
||||
.insert(compressedIdx.end(),
|
||||
this->compressedIdx_.begin() + *numOrig,
|
||||
this->compressedIdx_.end());
|
||||
}
|
||||
|
||||
this->compressedIdx_.swap(compressedIdx);
|
||||
}
|
||||
}
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// ---------------------------------------------------------------------
|
||||
// Class Opm::utility::CSRGraphFromCoordinates
|
||||
// ---------------------------------------------------------------------
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::clear()
|
||||
{
|
||||
this->uncompressed_.clear();
|
||||
this->csr_.clear();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
addConnection(const VertexID v1, const VertexID v2)
|
||||
{
|
||||
if ((v1 < 0) || (v2 < 0)) {
|
||||
throw std::invalid_argument {
|
||||
"Vertex IDs must be non-negative. Got (v1,v2) = ("
|
||||
+ std::to_string(v1) + ", " + std::to_string(v2)
|
||||
+ ')'
|
||||
};
|
||||
}
|
||||
|
||||
if (v1 == v2) {
|
||||
// Ignore self connections.
|
||||
return;
|
||||
}
|
||||
|
||||
this->uncompressed_.add(v1, v2);
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
void
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::
|
||||
compress(const Offset maxNumVertices, const bool expandExistingIdxMap)
|
||||
{
|
||||
if (! this->uncompressed_.isValid()) {
|
||||
throw std::logic_error {
|
||||
"Cannot compress invalid connection list"
|
||||
};
|
||||
}
|
||||
|
||||
this->csr_.merge(this->uncompressed_, maxNumVertices, expandExistingIdxMap);
|
||||
|
||||
this->uncompressed_.clear();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Offset
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::numVertices() const
|
||||
{
|
||||
return this->csr_.numRows();
|
||||
}
|
||||
|
||||
template <typename VertexID, bool TrackCompressedIdx>
|
||||
typename Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::Offset
|
||||
Opm::utility::CSRGraphFromCoordinates<VertexID, TrackCompressedIdx>::numEdges() const
|
||||
{
|
||||
const auto& ia = this->startPointers();
|
||||
|
||||
return ia.empty() ? 0 : ia.back();
|
||||
}
|
||||
@@ -24,6 +24,8 @@
|
||||
|
||||
#include <opm/output/data/InterRegFlow.hpp>
|
||||
|
||||
#include <opm/common/utility/CSRGraphFromCoordinates.hpp>
|
||||
|
||||
#include <cstddef>
|
||||
#include <optional>
|
||||
#include <utility>
|
||||
@@ -117,415 +119,62 @@ namespace Opm { namespace data {
|
||||
template <class MessageBufferType>
|
||||
void write(MessageBufferType& buffer) const
|
||||
{
|
||||
this->csr_.write(buffer);
|
||||
this->connections_.write(buffer);
|
||||
this->writeVector(this->rates_, buffer);
|
||||
}
|
||||
|
||||
// MessageBufferType API should be similar to Dune::MessageBufferIF
|
||||
template <class MessageBufferType>
|
||||
void read(MessageBufferType& buffer)
|
||||
{
|
||||
auto other = CSR{};
|
||||
other.read(buffer);
|
||||
this->connections_.read(buffer);
|
||||
|
||||
this->uncompressed_
|
||||
.add(other.maxRowIdx(),
|
||||
other.maxColIdx(),
|
||||
other.coordinateFormatRowIndices(),
|
||||
other.columnIndices(),
|
||||
other.values());
|
||||
auto rates = RateBuffer{};
|
||||
this->readVector(buffer, rates);
|
||||
this->appendRates(rates);
|
||||
}
|
||||
|
||||
/// Clear all internal buffers, but preserve allocated capacity.
|
||||
void clear();
|
||||
|
||||
private:
|
||||
/// Coordinate format representation of individual contributions to
|
||||
/// inter-region flows.
|
||||
class Connections
|
||||
// VertexID = int, TrackCompressedIdx = true.
|
||||
using Graph = utility::CSRGraphFromCoordinates<int, true>;
|
||||
|
||||
Graph connections_{};
|
||||
RateBuffer rates_{};
|
||||
|
||||
template <typename T, class A, class MessageBufferType>
|
||||
void writeVector(const std::vector<T,A>& vec,
|
||||
MessageBufferType& buffer) const
|
||||
{
|
||||
public:
|
||||
/// Add contributions from a single inter-region connection.
|
||||
///
|
||||
/// \param[in] r1 Source region. Zero-based region index/ID.
|
||||
///
|
||||
/// \param[in] r2 Destination region. Zero-based region index.
|
||||
///
|
||||
/// \param[in] rates Flow rates of single inter-region
|
||||
/// connection.
|
||||
void add(const int r1, const int r2, const FlowRates& rates);
|
||||
const auto n = vec.size();
|
||||
buffer.write(n);
|
||||
|
||||
/// Add contributions from multiple inter-region connections.
|
||||
///
|
||||
/// \param[in] maxRowIdx Maximum row (source region) index
|
||||
/// across all new inter-region connection contributions.
|
||||
///
|
||||
/// \param[in] maxColIdx Maximum column (destination region)
|
||||
/// index across all new inter-region contributions.
|
||||
///
|
||||
/// \param[in] rows Source region indices for all new
|
||||
/// inter-region connection contributions.
|
||||
///
|
||||
/// \param[in] cols Destination region indices for all new
|
||||
/// inter-region connection contributions.
|
||||
///
|
||||
/// \param[in] rates Flow rate values for all new inter-region
|
||||
/// connection contributions.
|
||||
void add(const int maxRowIdx,
|
||||
const int maxColIdx,
|
||||
const Neighbours& rows,
|
||||
const Neighbours& cols,
|
||||
const RateBuffer& rates);
|
||||
for (const auto& x : vec) {
|
||||
buffer.write(x);
|
||||
}
|
||||
}
|
||||
|
||||
/// Clear internal tables. Preserve allocated capacity.
|
||||
void clear();
|
||||
|
||||
/// Predicate.
|
||||
///
|
||||
/// \return Whether or not internal tables are empty.
|
||||
bool empty() const;
|
||||
|
||||
/// Whether or not internal tables meet size consistency
|
||||
/// requirements.
|
||||
bool isValid() const;
|
||||
|
||||
/// Maximum zero-based row (source region) index.
|
||||
int maxRow() const;
|
||||
|
||||
/// Maximum zero-based column (destination region) index.
|
||||
int maxCol() const;
|
||||
|
||||
/// Number of uncompressed contributions in internal tables.
|
||||
Neighbours::size_type numContributions() const;
|
||||
|
||||
/// Read-only access to uncompressed row indices.
|
||||
const Neighbours& rowIndices() const;
|
||||
|
||||
/// Read-only access to uncompressed column indices.
|
||||
const Neighbours& columnIndices() const;
|
||||
|
||||
/// Read-only access to uncompressed flow rate values.
|
||||
const RateBuffer& values() const;
|
||||
|
||||
private:
|
||||
/// Zero-based row/source region indices.
|
||||
Neighbours i_{};
|
||||
|
||||
/// Zero-based column/destination region indices.
|
||||
Neighbours j_{};
|
||||
|
||||
/// Uncompressed flow rate values. Window::bufferSize() entries
|
||||
/// per connection.
|
||||
RateBuffer v_{};
|
||||
|
||||
/// Maximum row index in \code this->i_ \endcode.
|
||||
int max_i_{ -1 };
|
||||
|
||||
/// Maximum column index in \code this->j_ \endcode.
|
||||
int max_j_{ -1 };
|
||||
};
|
||||
|
||||
/// Compressed sparse row representation of inter-region flow rates
|
||||
///
|
||||
/// Row and column indices are zero-based region IDs. Column
|
||||
/// indices ascendingly sorted per row. Value type is window,
|
||||
/// backed by a pair of iterators, of aggregate flow rates per
|
||||
/// region pair.
|
||||
class CSR
|
||||
template <typename T, class A, class MessageBufferType>
|
||||
void readVector(MessageBufferType& buffer,
|
||||
std::vector<T,A>& vec)
|
||||
{
|
||||
public:
|
||||
/// Merge coordinate format into existing CSR map.
|
||||
///
|
||||
/// \param[in] conns Coordinate representation of new
|
||||
/// contributions.
|
||||
///
|
||||
/// \param[in] numRegions Maximum number of regions in this
|
||||
/// region set. Common values/settings are
|
||||
///
|
||||
/// -# Maximum one-based region ID on local MPI rank
|
||||
/// -# Maximum one-based region ID across all MPI ranks
|
||||
/// -# Maximum *possible* one-based region ID in model
|
||||
/// ("NTFIP"), from TABDIMS(5) and/or REGDIMS(1).
|
||||
///
|
||||
/// If this value is smaller than the maximum one-based
|
||||
/// region ID on the local MPI rank, then it will be ignored
|
||||
/// and the local rank's maximum one-based region ID will be
|
||||
/// used instead.
|
||||
void merge(const Connections& conns,
|
||||
const Offset numRegions);
|
||||
auto n = 0 * vec.size();
|
||||
buffer.read(n);
|
||||
|
||||
/// Read-only access to flow rates of given region ID pair.
|
||||
///
|
||||
/// \param[in] i Source region. Zero-based region ID.
|
||||
///
|
||||
/// \param[in] j Destination region. Zero-based region ID.
|
||||
///
|
||||
/// \return Flow rates of region ID pair. Nullopt if no such
|
||||
/// pair exists.
|
||||
std::optional<ReadOnlyWindow> getWindow(const int i, const int j) const;
|
||||
vec.resize(n);
|
||||
|
||||
/// Total number of rows in compressed map structure.
|
||||
Offset numRows() const;
|
||||
|
||||
/// Maximum zero-based row index encountered mapped structure.
|
||||
int maxRowIdx() const;
|
||||
|
||||
/// Maximum zero-based column index encountered mapped structure.
|
||||
int maxColIdx() const;
|
||||
|
||||
/// Read-only access to compressed structure's start pointers.
|
||||
const Start& startPointers() const;
|
||||
|
||||
/// Read-only access to compressed structure's column indices,
|
||||
/// ascendingly sorted per rwo.
|
||||
const Neighbours& columnIndices() const;
|
||||
|
||||
/// Read-only access to compressed, unique, linearised flow rate
|
||||
/// values. \code Window::bufferSize() \endcode entries per
|
||||
/// non-zero element.
|
||||
const RateBuffer& values() const;
|
||||
|
||||
/// Coordinate format row index vector. Expanded from \code
|
||||
/// startPointers() \endcode.
|
||||
Neighbours coordinateFormatRowIndices() const;
|
||||
|
||||
// MessageBufferType API should be similar to Dune::MessageBufferIF
|
||||
template <class MessageBufferType>
|
||||
void write(MessageBufferType& buffer) const
|
||||
{
|
||||
this->writeVector(this->ia_, buffer);
|
||||
this->writeVector(this->ja_, buffer);
|
||||
this->writeVector(this->sa_, buffer);
|
||||
this->writeVector(this->compressedIdx_, buffer);
|
||||
|
||||
buffer.write(this->numRows_);
|
||||
buffer.write(this->numCols_);
|
||||
for (auto& x : vec) {
|
||||
buffer.read(x);
|
||||
}
|
||||
}
|
||||
|
||||
// MessageBufferType API should be similar to Dune::MessageBufferIF
|
||||
template <class MessageBufferType>
|
||||
void read(MessageBufferType& buffer)
|
||||
{
|
||||
this->readVector(buffer, this->ia_);
|
||||
this->readVector(buffer, this->ja_);
|
||||
this->readVector(buffer, this->sa_);
|
||||
this->readVector(buffer, this->compressedIdx_);
|
||||
|
||||
buffer.read(this->numRows_);
|
||||
buffer.read(this->numCols_);
|
||||
}
|
||||
|
||||
/// Clear internal tables. Preserve allocated capacity.
|
||||
void clear();
|
||||
|
||||
private:
|
||||
/// Start pointers.
|
||||
Start ia_{};
|
||||
|
||||
/// Column indices. Ascendingly sorted per row once structure
|
||||
/// is fully established.
|
||||
Neighbours ja_{};
|
||||
|
||||
/// Compressed, unique, linearised flow rate values. \code
|
||||
/// Window::bufferSize() \endcode entries per non-zero map
|
||||
/// element.
|
||||
RateBuffer sa_{};
|
||||
|
||||
/// Destination index in compressed representation. Size NNZ.
|
||||
Start compressedIdx_{};
|
||||
|
||||
/// Number of active rows in compressed map structure.
|
||||
int numRows_{ 0 };
|
||||
|
||||
/// Number of active columns in compressed map structure.
|
||||
/// Tracked as the maximum column index plus one.
|
||||
int numCols_{ 0 };
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// Implementation of read()/write()
|
||||
// ---------------------------------------------------------
|
||||
|
||||
template <typename T, class A, class MessageBufferType>
|
||||
void writeVector(const std::vector<T,A>& vec,
|
||||
MessageBufferType& buffer) const
|
||||
{
|
||||
const auto n = vec.size();
|
||||
buffer.write(n);
|
||||
|
||||
for (const auto& x : vec) {
|
||||
buffer.write(x);
|
||||
}
|
||||
}
|
||||
|
||||
template <class MessageBufferType, typename T, class A>
|
||||
void readVector(MessageBufferType& buffer,
|
||||
std::vector<T,A>& vec)
|
||||
{
|
||||
auto n = 0 * vec.size();
|
||||
buffer.read(n);
|
||||
|
||||
vec.resize(n);
|
||||
|
||||
for (auto& x : vec) {
|
||||
buffer.read(x);
|
||||
}
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// Implementation of merge()
|
||||
// ---------------------------------------------------------
|
||||
|
||||
/// Incorporate new, coordinate format contributions into
|
||||
/// existing, possibly empty, CSR mapping structure.
|
||||
///
|
||||
/// On exit the ia_ array holds the proper start pointers while
|
||||
/// ja_ holds the corresponding column indices albeit possibly
|
||||
/// repeated and unsorted.
|
||||
///
|
||||
/// \param[in] rows Row indices of all, possibly repeated,
|
||||
/// coordinate format input contributions. Start pointers \c
|
||||
/// ia_ updated to account for new entries.
|
||||
///
|
||||
/// \param[in] cols Column index of coordinate format intput
|
||||
/// structure. Inserted into \c ja_ according to its
|
||||
/// corresponding row index.
|
||||
///
|
||||
/// \param[in] maxRowIdx Maximum index in \p rows. Needed to
|
||||
/// ensure proper size of \c ia_.
|
||||
///
|
||||
/// \param[in] maxColIdx Maximum index in \p cols.
|
||||
void assemble(const Neighbours& rows,
|
||||
const Neighbours& cols,
|
||||
const int maxRowIdx,
|
||||
const int maxColIdx);
|
||||
|
||||
|
||||
/// Sort column indices per row and compress repeated column
|
||||
/// indices down to a single unique element per row. Sum
|
||||
/// repeated values
|
||||
///
|
||||
/// On exit the \c ia_, \c ja_, and \c sa_ arrays all have their
|
||||
/// expected, canonical structure.
|
||||
///
|
||||
/// \param[in] numRegions Maximum number of regions supported by
|
||||
/// final compressed mapping structure. Ignored if less than
|
||||
/// active number of rows.
|
||||
///
|
||||
/// \param[in] rates Uncompressed flow rate values from
|
||||
/// coordinate format contributions.
|
||||
void compress(const Offset numRegions,
|
||||
const RateBuffer& rates);
|
||||
|
||||
/// Sort column indices within each mapped row.
|
||||
///
|
||||
/// On exit \c ja_ has ascendingly sorted column indices, albeit
|
||||
/// possibly with repeated entries. This function also updates
|
||||
/// \c compressedIdx_ to account for the new locations of the
|
||||
/// non-zero elements in the grouped structure.
|
||||
void sortColumnIndicesPerRow();
|
||||
|
||||
/// Condense repeated column indices per row down to a single
|
||||
/// unique entry for each.
|
||||
///
|
||||
/// Assumes that each row has ascendingly sorted column indices
|
||||
/// in \c ja_ and must therefore be called after member function
|
||||
/// sortColumnIndicesPerRow(). On exit, \c ja_ has its final
|
||||
/// canonical structure and \c compressedIdx_ knows the final
|
||||
/// location of each non-zero contribution in the input
|
||||
/// coordinate format.
|
||||
void condenseDuplicates();
|
||||
|
||||
/// Sum coordinate format flow rates into compressed map
|
||||
/// structure.
|
||||
///
|
||||
/// Repeated (row,column) index pairs in the input coordinate
|
||||
/// format add to the same compressed map element. This
|
||||
/// function assumes that \c compressedIdx_ knows the final
|
||||
/// compressed location of each non-zero contribution in the
|
||||
/// input coordinate format and must therefore be called after
|
||||
/// member function condenseDuplicates(). On exit \c sa_ has
|
||||
/// incorporated all entries from the input coordinate
|
||||
/// structure.
|
||||
///
|
||||
/// \param[in] v Uncompressed flow rate values from coordinate
|
||||
/// format contributions.
|
||||
void accumulateFlowRates(const RateBuffer& v);
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// Implementation of assemble()
|
||||
// ---------------------------------------------------------
|
||||
|
||||
/// Position end pointers at start of row to prepare for column
|
||||
/// index grouping by corresponding row index.
|
||||
///
|
||||
/// Also counts total number of non-zero elements, possibly
|
||||
/// including repetitions, in \code this->ia_[0] \endcode.
|
||||
///
|
||||
/// \param[in] Number of rows in final compressed structure.
|
||||
/// Used to allocate \code this->ia_ \endcode.
|
||||
///
|
||||
/// \param[in] Row indices of all, possibly repeated, coordinate
|
||||
/// format input contributions. Needed to count the number
|
||||
/// of possibly repeated column index entries per row.
|
||||
void preparePushbackRowGrouping(const int numRows,
|
||||
const Neighbours& rowIdx);
|
||||
|
||||
/// Group column indices by corresponding row index and track
|
||||
/// grouped location of original coordinate format element
|
||||
///
|
||||
/// Appends grouped location to \c compressedIdx_.
|
||||
///
|
||||
/// \param[in] rowIdx Row index of coordinate format input
|
||||
/// structure. Used as grouping key.
|
||||
///
|
||||
/// \param[in] colIdx Column index of coordinate format intput
|
||||
/// structure. Inserted into \c ja_ according to its
|
||||
/// corresponding row index.
|
||||
void groupAndTrackColumnIndicesByRow(const Neighbours& rowIdx,
|
||||
const Neighbours& colIdx);
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// General utilities
|
||||
// ---------------------------------------------------------
|
||||
|
||||
/// Transpose connectivity structure.
|
||||
///
|
||||
/// Essentially swaps the roles of rows and columns. Also used
|
||||
/// as a basic building block for sortColumnIndicesPerRow().
|
||||
void transpose();
|
||||
|
||||
/// Condense sequences of repeated column indices in a single
|
||||
/// map row down to a single copy of each unique column index.
|
||||
///
|
||||
/// Appends new unique column indices to \code ja_ \endcode
|
||||
///
|
||||
/// Assumes that the map row has ascendingly sorted column
|
||||
/// indices and therefore has the same requirements as
|
||||
/// std::unique. Will also update the internal compressedIdx_
|
||||
/// mapping to record new compressed locations for the current,
|
||||
/// uncompressed, non-zero map elements.
|
||||
///
|
||||
/// \param[in] begin Start of map row that contains possibly
|
||||
/// repeated column indices.
|
||||
///
|
||||
/// \param[in] end One-past-end of map row that contains
|
||||
/// possibly repeated column indices.
|
||||
void condenseAndTrackUniqueColumnsForSingleRow(Neighbours::const_iterator begin,
|
||||
Neighbours::const_iterator end);
|
||||
|
||||
/// Update \c compressedIdx_ mapping to account for column index
|
||||
/// reshuffling.
|
||||
///
|
||||
/// \param[in] compressedIdx New compressed index locations of
|
||||
/// the non-zero map entries.
|
||||
void remapCompressedIndex(Start&& compressedIdx);
|
||||
};
|
||||
|
||||
/// Accumulated coordinate format contributions that have not yet
|
||||
/// been added to the final CSR structure.
|
||||
Connections uncompressed_;
|
||||
|
||||
/// Canonical representation of unique inter-region flow rates.
|
||||
CSR csr_;
|
||||
template <typename Rates>
|
||||
void appendRates(const Rates& rates)
|
||||
{
|
||||
this->rates_.insert(this->rates_.end(), rates.begin(), rates.end());
|
||||
}
|
||||
};
|
||||
|
||||
}} // namespace Opm::data
|
||||
|
||||
@@ -27,6 +27,8 @@
|
||||
|
||||
#include <opm/output/data/InterRegFlow.hpp>
|
||||
|
||||
#include <opm/common/utility/CSRGraphFromCoordinates.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <cassert>
|
||||
#include <exception>
|
||||
@@ -37,475 +39,6 @@
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
// ---------------------------------------------------------------------
|
||||
// Class Opm::data::InterRegFlowMap::Connections
|
||||
// ---------------------------------------------------------------------
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::
|
||||
Connections::add(const int r1, const int r2, const FlowRates& v)
|
||||
{
|
||||
using ElmT = Window::ElmT;
|
||||
|
||||
const auto one = ElmT{1};
|
||||
const auto sign = (r1 < r2) ? one : -one;
|
||||
|
||||
auto low = r1, high = r2;
|
||||
if (std::signbit(sign)) {
|
||||
std::swap(low, high);
|
||||
}
|
||||
|
||||
this->i_.push_back(low);
|
||||
this->j_.push_back(high);
|
||||
|
||||
this->max_i_ = std::max(this->max_i_, this->i_.back());
|
||||
this->max_j_ = std::max(this->max_j_, this->j_.back());
|
||||
|
||||
const auto start = this->v_.size();
|
||||
this->v_.insert(this->v_.end(), Window::bufferSize(), ElmT{0});
|
||||
|
||||
Window{ this->v_.begin() + start, this->v_.end() }.addFlow(sign, v);
|
||||
}
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::
|
||||
Connections::add(const int maxRowIdx,
|
||||
const int maxColIdx,
|
||||
const Neighbours& rows,
|
||||
const Neighbours& cols,
|
||||
const RateBuffer& rates)
|
||||
{
|
||||
if (cols.size() != rows.size()) {
|
||||
throw std::invalid_argument {
|
||||
"Coordinate format column index table size does not match "
|
||||
"row index table size"
|
||||
};
|
||||
}
|
||||
|
||||
if (rates.size() != Window::bufferSize() * rows.size()) {
|
||||
throw std::invalid_argument {
|
||||
"Coordinate format value table size does not match "
|
||||
"row index table size"
|
||||
};
|
||||
}
|
||||
|
||||
this->i_.insert(this->i_.end(), rows .begin(), rows .end());
|
||||
this->j_.insert(this->j_.end(), cols .begin(), cols .end());
|
||||
this->v_.insert(this->v_.end(), rates.begin(), rates.end());
|
||||
|
||||
this->max_i_ = std::max(this->max_i_, maxRowIdx);
|
||||
this->max_j_ = std::max(this->max_j_, maxColIdx);
|
||||
}
|
||||
|
||||
void Opm::data::InterRegFlowMap::Connections::clear()
|
||||
{
|
||||
this->v_.clear();
|
||||
this->j_.clear();
|
||||
this->i_.clear();
|
||||
|
||||
this->max_i_ = -1;
|
||||
this->max_j_ = -1;
|
||||
}
|
||||
|
||||
bool Opm::data::InterRegFlowMap::Connections::empty() const
|
||||
{
|
||||
return this->i_.empty();
|
||||
}
|
||||
|
||||
bool Opm::data::InterRegFlowMap::Connections::isValid() const
|
||||
{
|
||||
return (this->i_.size() == this->j_.size())
|
||||
&& (this->v_.size() == this->i_.size()*Window::bufferSize());
|
||||
}
|
||||
|
||||
int Opm::data::InterRegFlowMap::Connections::maxRow() const
|
||||
{
|
||||
return this->max_i_;
|
||||
}
|
||||
|
||||
int Opm::data::InterRegFlowMap::Connections::maxCol() const
|
||||
{
|
||||
return this->max_j_;
|
||||
}
|
||||
|
||||
Opm::data::InterRegFlowMap::Neighbours::size_type
|
||||
Opm::data::InterRegFlowMap::Connections::numContributions() const
|
||||
{
|
||||
return this->i_.size();
|
||||
}
|
||||
|
||||
const Opm::data::InterRegFlowMap::Neighbours&
|
||||
Opm::data::InterRegFlowMap::Connections::rowIndices() const
|
||||
{
|
||||
return this->i_;
|
||||
}
|
||||
|
||||
const Opm::data::InterRegFlowMap::Neighbours&
|
||||
Opm::data::InterRegFlowMap::Connections::columnIndices() const
|
||||
{
|
||||
return this->j_;
|
||||
}
|
||||
|
||||
const Opm::data::InterRegFlowMap::RateBuffer&
|
||||
Opm::data::InterRegFlowMap::Connections::values() const
|
||||
{
|
||||
return this->v_;
|
||||
}
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// ---------------------------------------------------------------------
|
||||
// Class Opm::data::InterRegFlowMap::CSR
|
||||
// ---------------------------------------------------------------------
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::
|
||||
CSR::merge(const Connections& conns, const Offset numRegions)
|
||||
{
|
||||
if (! conns.empty() &&
|
||||
(static_cast<Offset>(conns.maxRow()) >= numRegions))
|
||||
{
|
||||
throw std::invalid_argument {
|
||||
"Input graph contains more "
|
||||
"source regions than are "
|
||||
"implied by explicit size of "
|
||||
"adjacency matrix"
|
||||
};
|
||||
}
|
||||
|
||||
this->assemble(conns.rowIndices(), conns.columnIndices(),
|
||||
conns.maxRow(), conns.maxCol());
|
||||
this->compress(numRegions, conns.values());
|
||||
}
|
||||
|
||||
std::optional<Opm::data::InterRegFlowMap::ReadOnlyWindow>
|
||||
Opm::data::InterRegFlowMap::CSR::getWindow(const int i, const int j) const
|
||||
{
|
||||
if ((i < 0) || (i >= static_cast<int>(this->numRows())) ||
|
||||
(j < 0) || (j > this->maxColIdx()))
|
||||
{
|
||||
// Entity pair IDs out of range.
|
||||
return std::nullopt;
|
||||
}
|
||||
|
||||
auto begin = this->columnIndices().begin() + this->startPointers()[i + 0];
|
||||
auto end = this->columnIndices().begin() + this->startPointers()[i + 1];
|
||||
|
||||
auto pos = std::lower_bound(begin, end, j);
|
||||
if ((pos == end) || (*pos > j)) {
|
||||
// Entity 'j' does not connect to entity 'i'.
|
||||
return std::nullopt;
|
||||
}
|
||||
|
||||
// Entity 'j' connects to 'i'. Form read-only view into sub-range
|
||||
// pertaining to this entity pair.
|
||||
const auto sz = ReadOnlyWindow::bufferSize();
|
||||
const auto windowID = pos - this->columnIndices().begin();
|
||||
auto start = this->values().begin() + windowID*sz;
|
||||
|
||||
return { ReadOnlyWindow{ start, start + sz } };
|
||||
}
|
||||
|
||||
Opm::data::InterRegFlowMap::Offset
|
||||
Opm::data::InterRegFlowMap::CSR::numRows() const
|
||||
{
|
||||
return this->startPointers().empty()
|
||||
? 0 : this->startPointers().size() - 1;
|
||||
}
|
||||
|
||||
int Opm::data::InterRegFlowMap::CSR::maxRowIdx() const
|
||||
{
|
||||
return this->numRows_ - 1;
|
||||
}
|
||||
|
||||
int Opm::data::InterRegFlowMap::CSR::maxColIdx() const
|
||||
{
|
||||
return this->numCols_ - 1;
|
||||
}
|
||||
|
||||
const Opm::data::InterRegFlowMap::Start&
|
||||
Opm::data::InterRegFlowMap::CSR::startPointers() const
|
||||
{
|
||||
return this->ia_;
|
||||
}
|
||||
|
||||
const Opm::data::InterRegFlowMap::Neighbours&
|
||||
Opm::data::InterRegFlowMap::CSR::columnIndices() const
|
||||
{
|
||||
return this->ja_;
|
||||
}
|
||||
|
||||
const Opm::data::InterRegFlowMap::RateBuffer&
|
||||
Opm::data::InterRegFlowMap::CSR::values() const
|
||||
{
|
||||
return this->sa_;
|
||||
}
|
||||
|
||||
std::vector<int>
|
||||
Opm::data::InterRegFlowMap::CSR::coordinateFormatRowIndices() const
|
||||
{
|
||||
auto rowIdx = std::vector<int>{};
|
||||
|
||||
if (this->ia_.empty()) {
|
||||
return rowIdx;
|
||||
}
|
||||
|
||||
rowIdx.reserve(this->ia_.back());
|
||||
|
||||
auto row = 0;
|
||||
|
||||
const auto m = this->ia_.size() - 1;
|
||||
for (auto i = 0*m; i < m; ++i, ++row) {
|
||||
const auto n = this->ia_[i + 1] - this->ia_[i + 0];
|
||||
|
||||
rowIdx.insert(rowIdx.end(), n, row);
|
||||
}
|
||||
|
||||
return rowIdx;
|
||||
}
|
||||
|
||||
void Opm::data::InterRegFlowMap::CSR::clear()
|
||||
{
|
||||
this->ia_.clear();
|
||||
this->ja_.clear();
|
||||
this->sa_.clear();
|
||||
this->compressedIdx_.clear();
|
||||
|
||||
this->numRows_ = 0;
|
||||
this->numCols_ = 0;
|
||||
}
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::
|
||||
CSR::assemble(const Neighbours& rows,
|
||||
const Neighbours& cols,
|
||||
const int maxRowIdx,
|
||||
const int maxColIdx)
|
||||
{
|
||||
auto i = this->coordinateFormatRowIndices();
|
||||
i.insert(i.end(), rows.begin(), rows.end());
|
||||
|
||||
auto j = this->ja_;
|
||||
j.insert(j.end(), cols.begin(), cols.end());
|
||||
|
||||
const auto thisNumRows = std::max(this->numRows_, maxRowIdx + 1);
|
||||
const auto thisNumCols = std::max(this->numCols_, maxColIdx + 1);
|
||||
|
||||
this->preparePushbackRowGrouping(thisNumRows, i);
|
||||
|
||||
this->groupAndTrackColumnIndicesByRow(i, j);
|
||||
|
||||
this->numRows_ = thisNumRows;
|
||||
this->numCols_ = thisNumCols;
|
||||
}
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::
|
||||
CSR::compress(const Offset numRegions,
|
||||
const RateBuffer& rates)
|
||||
{
|
||||
this->sortColumnIndicesPerRow();
|
||||
|
||||
// Must be called *after* sortColumnIndicesPerRow().
|
||||
this->condenseDuplicates();
|
||||
|
||||
{
|
||||
auto v = this->values();
|
||||
v.insert(v.end(), rates.begin(), rates.end());
|
||||
|
||||
this->accumulateFlowRates(v);
|
||||
}
|
||||
|
||||
const auto nRows = this->startPointers().size() - 1;
|
||||
if (nRows < numRegions) {
|
||||
this->ia_.insert(this->ia_.end(),
|
||||
numRegions - nRows,
|
||||
this->startPointers().back());
|
||||
}
|
||||
}
|
||||
|
||||
void Opm::data::InterRegFlowMap::CSR::sortColumnIndicesPerRow()
|
||||
{
|
||||
// Transposition is, in this context, effectively a linear time (O(nnz))
|
||||
// bucket insertion procedure. In other words transposing the structure
|
||||
// twice creates a structure with column indices in (ascendingly) sorted
|
||||
// order.
|
||||
|
||||
this->transpose();
|
||||
this->transpose();
|
||||
}
|
||||
|
||||
void Opm::data::InterRegFlowMap::CSR::condenseDuplicates()
|
||||
{
|
||||
// Note: Must be called *after* sort().
|
||||
|
||||
const auto colIdx = this->ja_;
|
||||
auto compressedIdx = this->compressedIdx_;
|
||||
auto end = colIdx.begin();
|
||||
|
||||
this->ja_.clear();
|
||||
this->compressedIdx_.clear();
|
||||
|
||||
const auto numRows = this->ia_.size() - 1;
|
||||
for (auto row = 0*numRows; row < numRows; ++row) {
|
||||
auto begin = end;
|
||||
|
||||
std::advance(end, this->ia_[row + 1] - this->ia_[row + 0]);
|
||||
|
||||
const auto q = this->ja_.size();
|
||||
|
||||
this->condenseAndTrackUniqueColumnsForSingleRow(begin, end);
|
||||
|
||||
this->ia_[row + 0] = q;
|
||||
}
|
||||
|
||||
this->remapCompressedIndex(std::move(compressedIdx));
|
||||
|
||||
// Record final table sizes.
|
||||
this->ia_.back() = this->ja_.size();
|
||||
}
|
||||
|
||||
void Opm::data::InterRegFlowMap::CSR::accumulateFlowRates(const RateBuffer& v)
|
||||
{
|
||||
constexpr auto sz = Window::bufferSize();
|
||||
|
||||
if (v.size() != this->compressedIdx_.size()*sz) {
|
||||
throw std::logic_error {
|
||||
"Flow rates must be provided for each connection"
|
||||
};
|
||||
}
|
||||
|
||||
auto dst = [this](const Offset start) -> Window
|
||||
{
|
||||
auto begin = this->sa_.begin() + start*sz;
|
||||
|
||||
return Window { begin, begin + sz };
|
||||
};
|
||||
|
||||
auto src = [&v](const Offset start) -> ReadOnlyWindow
|
||||
{
|
||||
auto begin = v.begin() + start*sz;
|
||||
|
||||
return ReadOnlyWindow { begin, begin + sz };
|
||||
};
|
||||
|
||||
this->sa_.assign(this->ja_.size() * sz, Window::ElmT{0});
|
||||
|
||||
const auto numRates = this->compressedIdx_.size();
|
||||
for (auto rateID = 0*numRates; rateID < numRates; ++rateID) {
|
||||
dst(this->compressedIdx_[rateID]) += src(rateID);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::
|
||||
CSR::preparePushbackRowGrouping(const int numRows,
|
||||
const Neighbours& rowIdx)
|
||||
{
|
||||
assert (numRows >= 0);
|
||||
|
||||
this->ia_.assign(numRows + 1, 0);
|
||||
|
||||
for (const auto& row : rowIdx) {
|
||||
this->ia_[row + 1] += 1;
|
||||
}
|
||||
|
||||
// Note index range: 1..numRows inclusive.
|
||||
for (Start::size_type i = 1, n = numRows; i <= n; ++i) {
|
||||
this->ia_[0] += this->ia_[i];
|
||||
this->ia_[i] = this->ia_[0] - this->ia_[i];
|
||||
}
|
||||
|
||||
assert (this->ia_[0] == rowIdx.size());
|
||||
}
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::
|
||||
CSR::groupAndTrackColumnIndicesByRow(const Neighbours& rowIdx,
|
||||
const Neighbours& colIdx)
|
||||
{
|
||||
assert (this->ia_[0] == rowIdx.size());
|
||||
|
||||
const auto nnz = rowIdx.size();
|
||||
|
||||
this->ja_.resize(nnz);
|
||||
|
||||
this->compressedIdx_.clear();
|
||||
this->compressedIdx_.reserve(nnz);
|
||||
|
||||
for (auto nz = 0*nnz; nz < nnz; ++nz) {
|
||||
const auto k = this->ia_[rowIdx[nz] + 1] ++;
|
||||
|
||||
this->ja_[k] = colIdx[nz];
|
||||
this->compressedIdx_.push_back(k);
|
||||
}
|
||||
|
||||
this->ia_[0] = 0;
|
||||
}
|
||||
|
||||
void Opm::data::InterRegFlowMap::CSR::transpose()
|
||||
{
|
||||
auto compressedIdx = this->compressedIdx_;
|
||||
|
||||
{
|
||||
const auto rowIdx = this->coordinateFormatRowIndices();
|
||||
const auto colIdx = this->ja_;
|
||||
|
||||
this->preparePushbackRowGrouping(this->numCols_, colIdx);
|
||||
|
||||
// Note parameter order. Transposition switches role of rows and
|
||||
// columns.
|
||||
this->groupAndTrackColumnIndicesByRow(colIdx, rowIdx);
|
||||
}
|
||||
|
||||
this->remapCompressedIndex(std::move(compressedIdx));
|
||||
|
||||
std::swap(this->numRows_, this->numCols_);
|
||||
}
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::CSR::
|
||||
condenseAndTrackUniqueColumnsForSingleRow(Neighbours::const_iterator begin,
|
||||
Neighbours::const_iterator end)
|
||||
{
|
||||
// We assume that we're only called *after* sortColumnIndicesPerRow()
|
||||
// whence duplicate elements appear consecutively in [begin, end).
|
||||
//
|
||||
// Note: This is essentially the same as std::unique(begin, end) save
|
||||
// for the return value and the fact that we additionally record the
|
||||
// 'compressedIdx_' mapping. That mapping enables subsequent, decoupled
|
||||
// accumulation of the 'sa_' contributions.
|
||||
|
||||
while (begin != end) {
|
||||
// Note: Order of ja_ and compressedIdx_ matters here.
|
||||
this->compressedIdx_.push_back(this->ja_.size());
|
||||
this->ja_ .push_back(*begin);
|
||||
|
||||
while ((++begin != end) &&
|
||||
( *begin == this->ja_.back()))
|
||||
{
|
||||
this->compressedIdx_.push_back(this->compressedIdx_.back());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::CSR::
|
||||
remapCompressedIndex(Start&& compressedIdx)
|
||||
{
|
||||
for (auto& i : compressedIdx) {
|
||||
i = this->compressedIdx_[i];
|
||||
}
|
||||
|
||||
this->compressedIdx_.swap(compressedIdx);
|
||||
}
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// ---------------------------------------------------------------------
|
||||
// Class Opm::data::InterRegFlowMap
|
||||
// ---------------------------------------------------------------------
|
||||
|
||||
void
|
||||
Opm::data::InterRegFlowMap::
|
||||
addConnection(const int r1,
|
||||
@@ -525,26 +58,61 @@ addConnection(const int r1,
|
||||
return;
|
||||
}
|
||||
|
||||
this->uncompressed_.add(r1, r2, rates);
|
||||
const auto one = Window::ElmT{1};
|
||||
const auto sign = (r1 < r2) ? one : -one;
|
||||
const auto start = this->rates_.size();
|
||||
|
||||
auto low = r1, high = r2;
|
||||
if (std::signbit(sign)) {
|
||||
std::swap(low, high);
|
||||
}
|
||||
|
||||
this->connections_.addConnection(low, high);
|
||||
|
||||
this->rates_.insert(this->rates_.end(), Window::bufferSize(), Window::ElmT{0});
|
||||
Window { this->rates_.begin() + start, this->rates_.end() }.addFlow(sign, rates);
|
||||
}
|
||||
|
||||
void Opm::data::InterRegFlowMap::compress(const std::size_t numRegions)
|
||||
{
|
||||
if (! this->uncompressed_.isValid()) {
|
||||
this->connections_.compress(numRegions);
|
||||
|
||||
const auto v = this->rates_;
|
||||
constexpr auto sz = Window::bufferSize();
|
||||
const auto& dstIx = this->connections_.compressedIndexMap();
|
||||
|
||||
if (v.size() != dstIx.size()*sz) {
|
||||
throw std::logic_error {
|
||||
"Cannot compress invalid connection list"
|
||||
"Flow rates must be provided for each connection"
|
||||
};
|
||||
}
|
||||
|
||||
this->csr_.merge(this->uncompressed_, numRegions);
|
||||
auto dst = [this](const Offset start) -> Window
|
||||
{
|
||||
auto begin = this->rates_.begin() + start*sz;
|
||||
|
||||
this->uncompressed_.clear();
|
||||
return Window { begin, begin + sz };
|
||||
};
|
||||
|
||||
auto src = [&v](const Offset start) -> ReadOnlyWindow
|
||||
{
|
||||
auto begin = v.begin() + start*sz;
|
||||
|
||||
return ReadOnlyWindow { begin, begin + sz };
|
||||
};
|
||||
|
||||
this->rates_.assign(this->connections_.columnIndices().size() * sz, Window::ElmT{0});
|
||||
|
||||
const auto numRates = dstIx.size();
|
||||
for (auto rateID = 0*numRates; rateID < numRates; ++rateID) {
|
||||
dst(dstIx[rateID]) += src(rateID);
|
||||
}
|
||||
}
|
||||
|
||||
Opm::data::InterRegFlowMap::Offset
|
||||
Opm::data::InterRegFlowMap::numRegions() const
|
||||
{
|
||||
return this->csr_.numRows();
|
||||
return this->connections_.numVertices();
|
||||
}
|
||||
|
||||
std::optional<std::pair<
|
||||
@@ -580,17 +148,27 @@ Opm::data::InterRegFlowMap::getInterRegFlows(const int r1, const int r2) const
|
||||
std::swap(low, high);
|
||||
}
|
||||
|
||||
auto window = this->csr_.getWindow(low, high);
|
||||
if (! window.has_value()) {
|
||||
const auto& ia = this->connections_.startPointers();
|
||||
const auto& ja = this->connections_.columnIndices();
|
||||
|
||||
auto begin = ja.begin() + ia[low + 0];
|
||||
auto end = ja.begin() + ia[low + 1];
|
||||
auto pos = std::lower_bound(begin, end, high);
|
||||
if ((pos == end) || (*pos > high)) {
|
||||
// High is not connected to low.
|
||||
return std::nullopt;
|
||||
}
|
||||
|
||||
return std::make_pair(std::move(window.value()), sign);
|
||||
const auto sz = ReadOnlyWindow::bufferSize();
|
||||
const auto windowID = pos - ja.begin();
|
||||
|
||||
auto rateStart = this->rates_.begin() + windowID*sz;
|
||||
|
||||
return { std::pair { ReadOnlyWindow { rateStart, rateStart + sz }, sign } };
|
||||
}
|
||||
|
||||
void Opm::data::InterRegFlowMap::clear()
|
||||
{
|
||||
this->uncompressed_.clear();
|
||||
this->csr_.clear();
|
||||
this->connections_.clear();
|
||||
this->rates_.clear();
|
||||
}
|
||||
|
||||
1449
tests/test_CSRGraphFromCoordinates.cpp
Normal file
1449
tests/test_CSRGraphFromCoordinates.cpp
Normal file
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user