OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Fixed unsafe type conversion.

Browse Source
This commit is contained in:
Kjetil Olsen Lye
2023-03-29 14:21:41 +02:00
parent 3f1cbbe1b8
commit cc8db08e83
2 changed files with 68 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

62
opm/simulators/linalg/cuistl/CuSparseMatrix.cpp
View File

@@ -44,15 +44,15 @@ template <class T>
CuSparseMatrix<T>::CuSparseMatrix(const T* nonZeroElements,
const int* rowIndices,
const int* columnIndices,
int numberOfNonzeroBlocks,
int blockSize,
int numberOfRows)
size_t numberOfNonzeroBlocks,
size_t blockSize,
size_t numberOfRows)
: m_nonZeroElements(nonZeroElements, numberOfNonzeroBlocks * blockSize * blockSize)
, m_columnIndices(columnIndices, numberOfNonzeroBlocks)
, m_rowIndices(rowIndices, numberOfRows + 1)
, m_numberOfNonzeroBlocks(numberOfNonzeroBlocks)
, m_numberOfRows(numberOfRows)
, m_blockSize(blockSize)
, m_numberOfNonzeroBlocks(detail::to_int(numberOfNonzeroBlocks))
, m_numberOfRows(detail::to_int(numberOfRows))
, m_blockSize(detail::to_int(blockSize))
, m_matrixDescription(detail::createMatrixDescription())
, m_cusparseHandle(detail::CuSparseHandle::getInstance())
{
@@ -75,11 +75,11 @@ CuSparseMatrix<T>::fromMatrix(const MatrixType& matrix)
rowIndices.push_back(0);
const int blockSize = matrix[0][0].N();
const int numberOfRows = matrix.N();
const int numberOfNonzeroBlocks = matrix.nonzeroes();
const int numberOfNonzeroElements = blockSize * blockSize * numberOfNonzeroBlocks;
const size_t blockSize = matrix[0][0].N();
const size_t numberOfRows = matrix.N();
const size_t numberOfNonzeroBlocks = matrix.nonzeroes();
#ifndef CUSPARSE_ASSUME_UNSAFE_SPARSITY
const size_t numberOfNonzeroElements = blockSize * blockSize * numberOfNonzeroBlocks;
std::vector<T> nonZeroElementsData;
// TODO: [perf] Can we avoid building nonZeroElementsData?
nonZeroElementsData.reserve(numberOfNonzeroElements);
@@ -90,14 +90,14 @@ CuSparseMatrix<T>::fromMatrix(const MatrixType& matrix)
for (auto columnIterator = row.begin(); columnIterator != row.end(); ++columnIterator) {
columnIndices.push_back(columnIterator.index());
#ifndef CUSPARSE_ASSUME_UNSAFE_SPARSITY
for (int c = 0; c < blockSize; ++c) {
for (int d = 0; d < blockSize; ++d) {
for (size_t c = 0; c < blockSize; ++c) {
for (size_t d = 0; d < blockSize; ++d) {
nonZeroElementsData.push_back((*columnIterator)[c][d]);
}
}
#endif
}
rowIndices.push_back(columnIndices.size());
rowIndices.push_back(detail::to_int(columnIndices.size()));
}
#ifndef CUSPARSE_ASSUME_UNSAFE_SPARSITY
auto nonZeroElements = nonZeroElementsData.data();
@@ -106,7 +106,7 @@ CuSparseMatrix<T>::fromMatrix(const MatrixType& matrix)
#endif
// Sanity check
// h_rows and h_cols could be changed to 'unsigned int', but cusparse expects 'int'
OPM_ERROR_IF(size_t(rowIndices[matrix.N()]) != matrix.nonzeroes(),
OPM_ERROR_IF(rowIndices[matrix.N()] != detail::to_int(matrix.nonzeroes()),
"Error size of rows do not sum to number of nonzeroes in CuSparseMatrix.");
OPM_ERROR_IF(rowIndices.size() != numberOfRows + 1, "Row indices do not match for CuSparseMatrix.");
OPM_ERROR_IF(columnIndices.size() != numberOfNonzeroBlocks, "Column indices do not match for CuSparseMatrix.");
@@ -125,17 +125,18 @@ CuSparseMatrix<T>::updateNonzeroValues(const MatrixType& matrix)
OPM_ERROR_IF(matrix[0][0].N() != blockSize(), "Matrix does not have the same blocksize.");
OPM_ERROR_IF(matrix.N() != N(), "Matrix does not have the same number of rows.");
const int numberOfRows = N();
const int numberOfNonzeroBlocks = nonzeroes();
const int numberOfNonzeroElements = blockSize() * blockSize() * numberOfNonzeroBlocks;
#ifndef CUSPARSE_ASSUME_UNSAFE_SPARSITY
const size_t numberOfRows = N();
const size_t numberOfNonzeroBlocks = nonzeroes();
const size_t numberOfNonzeroElements = blockSize() * blockSize() * numberOfNonzeroBlocks;
std::vector<T> nonZeroElementsData;
// TODO: [perf] Can we avoid building nonZeroElementsData?
nonZeroElementsData.reserve(numberOfNonzeroElements);
for (auto& row : matrix) {
for (auto columnIterator = row.begin(); columnIterator != row.end(); ++columnIterator) {
for (int c = 0; c < blockSize(); ++c) {
for (int d = 0; d < blockSize(); ++d) {
for (size_t c = 0; c < blockSize(); ++c) {
for (size_t d = 0; d < blockSize(); ++d) {
nonZeroElementsData.push_back((*columnIterator)[c][d]);
}
}
@@ -182,13 +183,11 @@ CuSparseMatrix<T>::mv(const CuVector<T>& x, CuVector<T>& y) const
{
CHECK_SIZE(x)
CHECK_SIZE(y)
if (blockSize() < 2) {
if (blockSize() < 2u) {
OPM_THROW(
std::invalid_argument,
"CuSparseMatrix<T>::usmv and CuSparseMatrix<T>::mv are only implemented for block sizes greater than 1.");
}
const auto numberOfRows = N();
const auto numberOfNonzeroBlocks = nonzeroes();
const auto nonzeroValues = getNonZeroValues().data();
auto rowIndices = getRowIndices().data();
@@ -198,9 +197,9 @@ CuSparseMatrix<T>::mv(const CuVector<T>& x, CuVector<T>& y) const
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrmv(m_cusparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
CUSPARSE_OPERATION_NON_TRANSPOSE,
numberOfRows,
numberOfRows,
numberOfNonzeroBlocks,
m_numberOfRows,
m_numberOfRows,
m_numberOfNonzeroBlocks,
&alpha,
m_matrixDescription->get(),
nonzeroValues,
@@ -218,13 +217,12 @@ CuSparseMatrix<T>::umv(const CuVector<T>& x, CuVector<T>& y) const
{
CHECK_SIZE(x)
CHECK_SIZE(y)
if (blockSize() < 2) {
if (blockSize() < 2u) {
OPM_THROW(
std::invalid_argument,
"CuSparseMatrix<T>::usmv and CuSparseMatrix<T>::mv are only implemented for block sizes greater than 1.");
}
const auto numberOfRows = N();
const auto numberOfNonzeroBlocks = nonzeroes();
const auto nonzeroValues = getNonZeroValues().data();
auto rowIndices = getRowIndices().data();
@@ -234,15 +232,15 @@ CuSparseMatrix<T>::umv(const CuVector<T>& x, CuVector<T>& y) const
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrmv(m_cusparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
CUSPARSE_OPERATION_NON_TRANSPOSE,
numberOfRows,
numberOfRows,
numberOfNonzeroBlocks,
m_numberOfRows,
m_numberOfRows,
m_numberOfNonzeroBlocks,
&alpha,
m_matrixDescription->get(),
nonzeroValues,
rowIndices,
columnIndices,
blockSize(),
m_blockSize,
x.data(),
&beta,
y.data()));

47
opm/simulators/linalg/cuistl/CuSparseMatrix.hpp
View File

@@ -25,6 +25,7 @@
#include <opm/simulators/linalg/cuistl/CuVector.hpp>
#include <opm/simulators/linalg/cuistl/detail/CuMatrixDescription.hpp>
#include <opm/simulators/linalg/cuistl/detail/CuSparseHandle.hpp>
#include <opm/simulators/linalg/cuistl/detail/safe_conversion.hpp>
#include <vector>
namespace Opm::cuistl
@@ -56,12 +57,15 @@ public:
//! \param[in] numberOfNonzeroElements number of nonzero elements
//! \param[in] blockSize size of each block matrix (typically 3)
//! \param[in] numberOfRows the number of rows
//!
//! \note We assume numberOfNonzeroBlocks, blockSize and numberOfRows all are representable as int due to
//! restrictions in the current version of cusparse. This might change in future versions.
CuSparseMatrix(const T* nonZeroElements,
const int* rowIndices,
const int* columnIndices,
int numberOfNonzeroBlocks,
int blockSize,
int numberOfRows);
size_t numberOfNonzeroBlocks,
size_t blockSize,
size_t numberOfRows);
// TODO: Handle copy ctor and operator=
@@ -101,7 +105,12 @@ public:
*/
size_t N() const
{
return m_numberOfRows;
// Technically this safe conversion is not needed since we enforce these to be
// non-negative in the constructor, but keeping them for added sanity for now.
//
// We don't believe this will yield any performance penality (it's used too far away from the inner loop),
// but should that be false, they can be removed.
return detail::to_size_t(m_numberOfRows);
}
/**
@@ -110,7 +119,12 @@ public:
*/
size_t nonzeroes() const
{
return m_numberOfNonzeroBlocks;
// Technically this safe conversion is not needed since we enforce these to be
// non-negative in the constructor, but keeping them for added sanity for now.
//
// We don't believe this will yield any performance penality (it's used too far away from the inner loop),
// but should that be false, they can be removed.
return detail::to_size_t(m_numberOfNonzeroBlocks);
}
/**
@@ -179,17 +193,27 @@ public:
* This is equivalent to matrix.N() * matrix.blockSize()
* @return matrix.N() * matrix.blockSize()
*/
int dim() const
size_t dim() const
{
return m_blockSize * m_numberOfRows;
// Technically this safe conversion is not needed since we enforce these to be
// non-negative in the constructor, but keeping them for added sanity for now.
//
// We don't believe this will yield any performance penality (it's used too far away from the inner loop),
// but should that be false, they can be removed.
return detail::to_size_t(m_blockSize) * detail::to_size_t(m_numberOfRows);
}
/**
* @brief blockSize size of the blocks
*/
int blockSize() const
size_t blockSize() const
{
return m_blockSize;
// Technically this safe conversion is not needed since we enforce these to be
// non-negative in the constructor, but keeping them for added sanity for now.
//
// We don't believe this will yield any performance penality (it's used too far away from the inner loop),
// but should that be false, they can be removed.
return detail::to_size_t(m_blockSize);
}
/**
@@ -244,6 +268,11 @@ private:
CuVector<T> m_nonZeroElements;
CuVector<int> m_columnIndices;
CuVector<int> m_rowIndices;
// Notice that we store these three as int to make sure we are cusparse compatible.
//
// This gives the added benifit of checking the size constraints at construction of the matrix
// rather than in some call to cusparse.
const int m_numberOfNonzeroBlocks;
const int m_numberOfRows;
const int m_blockSize;