CuSparseMatrix.hpp

Go to the documentation of this file.

/*
  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_CUSPARSEMATRIX_HPP
#define OPM_CUSPARSEMATRIX_HPP
#include <cusparse.h>
#include <iostream>
#include <memory>
#include <opm/common/ErrorMacros.hpp>
#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
{
 
template <typename T>
class CuSparseMatrix
{
public:
    CuSparseMatrix(const T* nonZeroElements,
                   const int* rowIndices,
                   const int* columnIndices,
                   size_t numberOfNonzeroBlocks,
                   size_t blockSize,
                   size_t numberOfRows);
 
    CuSparseMatrix(const CuSparseMatrix&) = delete;
 
    CuSparseMatrix& operator=(const CuSparseMatrix&) = delete;
 
    virtual ~CuSparseMatrix();
 
    template <class MatrixType>
    static CuSparseMatrix<T> fromMatrix(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false);
 
    void setUpperTriangular();
 
    void setLowerTriangular();
 
    void setUnitDiagonal();
 
 
    void setNonUnitDiagonal();
 
    size_t N() const
    {
        // 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);
    }
 
    size_t nonzeroes() const
    {
        // 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);
    }
 
    CuVector<T>& getNonZeroValues()
    {
        return m_nonZeroElements;
    }
 
    const CuVector<T>& getNonZeroValues() const
    {
        return m_nonZeroElements;
    }
 
    CuVector<int>& getRowIndices()
    {
        return m_rowIndices;
    }
 
    const CuVector<int>& getRowIndices() const
    {
        return m_rowIndices;
    }
 
    CuVector<int>& getColumnIndices()
    {
        return m_columnIndices;
    }
 
    const CuVector<int>& getColumnIndices() const
    {
        return m_columnIndices;
    }
 
    size_t dim() const
    {
        // 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);
    }
 
    size_t blockSize() const
    {
        // 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::CuSparseMatrixDescription& getDescription()
    {
        return *m_matrixDescription;
    }
 
    virtual void mv(const CuVector<T>& x, CuVector<T>& y) const;
 
    virtual void umv(const CuVector<T>& x, CuVector<T>& y) const;
 
 
    virtual void usmv(T alpha, const CuVector<T>& x, CuVector<T>& y) const;
 
    template <class MatrixType>
    void updateNonzeroValues(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false);
 
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 benefit 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;
 
    detail::CuSparseMatrixDescriptionPtr m_matrixDescription;
    detail::CuSparseHandle& m_cusparseHandle;
 
    template <class VectorType>
    void assertSameSize(const VectorType& vector) const;
};
} // namespace Opm::cuistl
#endif