Loading...
Searching...
No Matches
GpuVector.hpp
/*
Copyright 2022-2023 SINTEF AS
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_GPUVECTOR_HEADER_HPP
#define OPM_GPUVECTOR_HEADER_HPP
#include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
#include <fmt/core.h>
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/gpuistl/detail/CuBlasHandle.hpp>
#include <opm/simulators/linalg/gpuistl/detail/safe_conversion.hpp>
#include <opm/simulators/linalg/gpuistl/detail/gpu_constants.hpp>
#include <vector>
#include <string>
namespace Opm::gpuistl
{
template <typename T>
class GpuVector
{
public:
using field_type = T;
using size_type = size_t;
GpuVector(const GpuVector<T>& other);
explicit GpuVector(const std::vector<T>& data);
GpuVector& operator=(const GpuVector<T>& other);
template<int BlockDimension>
explicit GpuVector(const Dune::BlockVector<Dune::FieldVector<T, BlockDimension>>& bvector)
: GpuVector(bvector.dim())
{
copyFromHost(bvector);
}
GpuVector& operator=(T scalar);
GpuVector() : m_dataOnDevice(nullptr), m_numberOfElements(0), m_cuBlasHandle(detail::CuBlasHandle::getInstance()) {}
explicit GpuVector(const size_t numberOfElements);
GpuVector(const T* dataOnHost, const size_t numberOfElements);
virtual ~GpuVector();
T* data();
const T* data() const;
template <int BlockDimension>
void copyFromHost(const Dune::BlockVector<Dune::FieldVector<T, BlockDimension>>& bvector)
{
// TODO: [perf] vector.dim() can be replaced by bvector.N() * BlockDimension
if (detail::to_size_t(m_numberOfElements) != bvector.dim()) {
OPM_THROW(std::runtime_error,
fmt::format("Given incompatible vector size. GpuVector has size {}, \n"
"however, BlockVector has N() = {}, and dim = {}.",
m_numberOfElements,
bvector.N(),
bvector.dim()));
}
const auto dataPointer = static_cast<const T*>(&(bvector[0][0]));
copyFromHost(dataPointer, m_numberOfElements);
}
template <int BlockDimension>
void copyFromHostAsync(const Dune::BlockVector<Dune::FieldVector<T, BlockDimension>>& bvector, cudaStream_t stream = detail::DEFAULT_STREAM)
{
// TODO: [perf] vector.dim() can be replaced by bvector.N() * BlockDimension
if (detail::to_size_t(m_numberOfElements) != bvector.dim()) {
OPM_THROW(std::runtime_error,
fmt::format("Given incompatible vector size. GpuVector has size {}, \n"
"however, BlockVector has N() = {}, and dim = {}.",
m_numberOfElements,
bvector.N(),
bvector.dim()));
}
const auto dataPointer = static_cast<const T*>(&(bvector[0][0]));
copyFromHostAsync(dataPointer, m_numberOfElements, stream);
}
template <int BlockDimension>
void copyToHost(Dune::BlockVector<Dune::FieldVector<T, BlockDimension>>& bvector) const
{
// TODO: [perf] vector.dim() can be replaced by bvector.N() * BlockDimension
if (detail::to_size_t(m_numberOfElements) != bvector.dim()) {
OPM_THROW(std::runtime_error,
fmt::format("Given incompatible vector size. GpuVector has size {},\n however, the BlockVector "
"has has N() = {}, and dim() = {}.",
m_numberOfElements,
bvector.N(),
bvector.dim()));
}
const auto dataPointer = static_cast<T*>(&(bvector[0][0]));
copyToHost(dataPointer, m_numberOfElements);
}
template <int BlockDimension>
void copyToHostAsync(Dune::BlockVector<Dune::FieldVector<T, BlockDimension>>& bvector, cudaStream_t stream = detail::DEFAULT_STREAM) const
{
// TODO: [perf] vector.dim() can be replaced by bvector.N() * BlockDimension
if (detail::to_size_t(m_numberOfElements) != bvector.dim()) {
OPM_THROW(std::runtime_error,
fmt::format("Given incompatible vector size. GpuVector has size {},\n however, the BlockVector "
"has has N() = {}, and dim() = {}.",
m_numberOfElements,
bvector.N(),
bvector.dim()));
}
const auto dataPointer = static_cast<T*>(&(bvector[0][0]));
copyToHostAsync(dataPointer, m_numberOfElements, stream);
}
void copyFromHost(const T* dataPointer, size_t numberOfElements);
void copyFromHostAsync(const T* dataPointer, size_t numberOfElements, cudaStream_t stream = detail::DEFAULT_STREAM);
void copyToHost(T* dataPointer, size_t numberOfElements) const;
void copyToHostAsync(T* dataPointer, size_t numberOfElements, cudaStream_t stream = detail::DEFAULT_STREAM) const;
void copyFromHost(const std::vector<T>& data);
void copyFromHostAsync(const std::vector<T>& data, cudaStream_t stream = detail::DEFAULT_STREAM);
void copyToHost(std::vector<T>& data) const;
void copyToHostAsync(std::vector<T>& data, cudaStream_t stream = detail::DEFAULT_STREAM) const;
void copyFromDeviceToDevice(const GpuVector<T>& other) const;
void prepareSendBuf(GpuVector<T>& buffer, const GpuVector<int>& indexSet) const;
void syncFromRecvBuf(GpuVector<T>& buffer, const GpuVector<int>& indexSet) const;
GpuVector<T>& operator*=(const T& scalar);
GpuVector<T>& axpy(T alpha, const GpuVector<T>& y);
GpuVector<T>& operator+=(const GpuVector<T>& other);
GpuVector<T>& operator-=(const GpuVector<T>& other);
T dot(const GpuVector<T>& other) const;
T two_norm() const;
T dot(const GpuVector<T>& other, const GpuVector<int>& indexSet, GpuVector<T>& buffer) const;
T two_norm(const GpuVector<int>& indexSet, GpuVector<T>& buffer) const;
T dot(const GpuVector<T>& other, const GpuVector<int>& indexSet) const;
T two_norm(const GpuVector<int>& indexSet) const;
size_type dim() const;
void resize(size_t new_size);
std::vector<T> asStdVector() const;
template <int blockSize>
Dune::BlockVector<Dune::FieldVector<T, blockSize>> asDuneBlockVector() const
{
OPM_ERROR_IF(dim() % blockSize != 0,
fmt::format("blockSize is not a multiple of dim(). Given blockSize = {}, and dim() = {}",
blockSize,
dim()));
Dune::BlockVector<Dune::FieldVector<T, blockSize>> returnValue(dim() / blockSize);
copyToHost(returnValue);
return returnValue;
}
void setZeroAtIndexSet(const GpuVector<int>& indexSet);
// Slow method that creates a string representation of a GpuVector for debug purposes
std::string toDebugString()
{
std::vector<T> v = asStdVector();
std::string res = "";
for (T element : v){
res += std::to_string(element) + " ";
}
res += std::to_string(v[v.size()-1]);
return res;
}
private:
T* m_dataOnDevice = nullptr;
// Note that we store this as int to make sure we are always cublas compatible.
// This gives the added benefit that a size_t to int conversion error occurs during construction.
int m_numberOfElements;
detail::CuBlasHandle& m_cuBlasHandle;
void assertSameSize(const GpuVector<T>& other) const;
void assertSameSize(int size) const;
void assertHasElements() const;
};
} // namespace Opm::gpuistl
// ADL bridge: convert GPU vector to Dune BlockVector and delegate to Dune's writer
namespace Opm::gpuistl
{
template <typename T>
inline void writeMatrixMarket(const GpuVector<T>& vectorOnDevice, std::ostream& ostr)
{
const auto hostBlockVector = vectorOnDevice.template asDuneBlockVector<1>();
writeMatrixMarket(hostBlockVector, ostr);
}
} // namespace Opm::gpuistl
#endif
Opm::gpuistl::GpuVector
Definition gpu_type_detection.hpp:30
Opm::gpuistl::detail::setZeroAtIndexSet
void setZeroAtIndexSet(T *deviceData, size_t numberOfElements, const int *indices)
setZeroAtIndexSet sets deviceData to zero in the indices of contained in indices
Opm::gpuistl::detail::DEFAULT_STREAM
constexpr cudaStream_t DEFAULT_STREAM
The default GPU stream (stream 0).
Definition gpu_constants.hpp:31
Opm::gpuistl::detail::to_size_t
__host__ __device__ std::size_t to_size_t(int i)
to_size_t converts a (on most relevant platforms) a 32 bit signed int to a 64 bits unsigned int
Definition safe_conversion.hpp:86