diff --git a/src/TNL/Algorithms/Segments/SegmentsPrinting.h b/src/TNL/Algorithms/Segments/SegmentsPrinting.h index f8fd7412edc134c3c5454d2fca819b53b16e0474..491f059fa6b66c300f68d37f06a25772ddfd2b78 100644 --- a/src/TNL/Algorithms/Segments/SegmentsPrinting.h +++ b/src/TNL/Algorithms/Segments/SegmentsPrinting.h @@ -49,6 +49,8 @@ std::ostream& printSegments( const Segments& segments, std::ostream& str ) return str; } +/// This is to prevent from appearing in Doxygen documentation. +/// \cond HIDDEN_CLASS template< typename Segments, typename Fetch > struct SegmentsPrinter @@ -130,6 +132,7 @@ std::ostream& printSegments( const Segments& segments, Fetch&& fetch, std::ostre } return str; } +/// \endcond } // namespace Segments } // namespace Algorithms diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h new file mode 100644 index 0000000000000000000000000000000000000000..6a2a6a5652e286fe5c1152d02e4c9f903880f046 --- /dev/null +++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h @@ -0,0 +1,1176 @@ +/*************************************************************************** + SparseSandboxMatrix.h - description + ------------------- + begin : Apr 19, 2021 + copyright : (C) 2021 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +#include <map> +#include <TNL/Matrices/Matrix.h> +#include <TNL/Matrices/MatrixType.h> +#include <TNL/Allocators/Default.h> +#include <TNL/Algorithms/Segments/CSR.h> +#include <TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h> +#include <TNL/Matrices/Sandbox/SparseSandboxMatrixView.h> +#include <TNL/Matrices/DenseMatrix.h> + +namespace TNL { + namespace Matrices { + namespace Sandbox { + +/** + * \brief Template of a sparse matrix that can be used for testing of new sparse-matrix formats. + * + * \tparam Real is a type of matrix elements. If \e Real equals \e bool the matrix is treated + * as binary and so the matrix elements values are not stored in the memory since we need + * to remember only coordinates of non-zero elements( which equal one). + * \tparam Device is a device where the matrix is allocated. + * \tparam Index is a type for indexing of the matrix elements. + * \tparam MatrixType specifies a symmetry of matrix. See \ref MatrixType. Symmetric + * matrices store only lower part of the matrix and its diagonal. The upper part is reconstructed on the fly. + * GeneralMatrix with no symmetry is used by default. + * \tparam RealAllocator is allocator for the matrix elements values. + * \tparam IndexAllocator is allocator for the matrix elements column indexes. + * + * This class can be used for rapid testing and development of new formats for sparse matrices. One may profit from + * several TNL tools compatible with interface of this templated class like: + * + * 1. Large set of existing unit tests. + * 3. Matrix reading from MTX files - to use \ref TNL::Matrices::MatrixReader, the following methods must be functional + * a. \ref TNL::Matrices::SandboxSparseMatrix::setRowCapacities + * b. \ref TNL::Matrices::SandboxSparseMatrix::setElement + * c. \ref TNL::Matrices::SandboxSparseMatrix::operator= between different devices + * 4. Matrix benchmarks - the following methods must be functional + * a. \ref TNL::Matrices::SandboxSparseMatrix::vectorProduct - for SpMV benchmark + * 5. Linear solvers + * 6. Simple comparison of performance with other matrix formats + * + * In the core of this class there is: + * + * 1. Vector 'values` (\ref TNL::Matrices::Matrix::values) which is inheritted from \ref TNL::Matrices::Matrix. This vector is used for storing + * of matrix elements values. + * 2. Vector `columnIndexes` (\ref TNL::Matrices::SendboxMatrix::columnIndexes). This vector is used for storing of matrix elements column indexes. + * + * This class contains fully functional implementation of CSR format and so the user have to replace just what he needs to. Once you have + * successfully implemented the sparse matrix format in this form, you may consider to extract it into a form of segments to make it accessible + * even for other algorithms then SpMV. + * + * Parts of the code, that need to be modified are marked by SANDBOX_TODO tag. The whole implementation consits of the following classes: + * + * 1. \ref TNL::Matrices::Sandbox::SparseSandboxMatrix - this class, it serves for matrix setup and performing of the main operations. + * 2. \ref TNL::Matrices::Sandbox::SparseSandboxMatrixView - view class which is necessary mainly for passing the matrix to GPU kernels. Most methods of `SparseSandboxMatrix` are common + * with `SparseSandboxMatrixView` and in this case they are implemented in the view class (and there is just redirection from this class). For this reason, `SparseSandboxMatrix` contains instance of the view class + * (\ref TNL::Matrices::Sandbox::SparseSandboxMatrix::view) which needs to be regularly updated each time when metadata are changed. This is usually done by the means of + * method \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::getView. + * 3. \ref TNL::Matrices::Sandbox::SparseSandboxMatrixRowView - is a class for accessing particular matrix rows. It will, likely, require some changes as well. + * + * We suggest the following way of implementation of the new sparse matrix format: + * + * 1. Add metadata required by your format next to \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::rowPointers but do not replace the row pointers. It will allow you + * to implement your new format next to the original CSR and to check/compare with the valid CSR implementation any time you get into troubles. The adventage is that all + * unit tests are working properly and you may just focus on modifying one method after another. The unit tests are called from + * `src/UnitTests/Matrices/SparseMatrixTests_SandboxMatrix.h` and `src/UnitTests/Matrices/SparseMatrixVectorProductTests_SandboxMatrix.h` + * 2. Modify first the method \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::setRowCapacities which is responsible for the setup of the format metadata. + * 3. Continue with modification of constructors, view class, \ref TNL::Matrices::Sandbox::SparseSandoxMatrix::getView and \ref TNL::Matrices::Sandbox::SparseSandoxMatrix::getConstView. + * 4. Next you need to modify \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::setElement and \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::getElement methods and assignment operator + * at least for copying the matrix across different devices (i.e. from CPU to GPU). It will allow you to use \ref TNL::Matrices::MatrixReader. We recommend to have the same data layout + * on both CPU and GPU so that the transfer of the matrix from CPU to GPU is trivial. + * 5. Finally proceed to \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::vectorProduct to implement SpMV operation. We recommend to implement first the CPU version which is easier to + * debug. Next proceed to GPU version. + * 6. When SpMV works it is time to delete the original CSR implementation, i.e. everything around `rowPointers`. + * 7. Optimize your implementation to the best performance and test with `tnl-benchmark-spmv` - you need to include your new matrix to `src/Benchmarks/SpMV/spmv.h` and modify this file + * accordingly. + * 8. If you want, you may now generalize SpMV to \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::reduceRows method. + * 9. If you have `reduceRows` implemented, you may use the original implementation of SpMV based just on the `reduceRows` method. + * 10. You may implement \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::forRows and \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::forElements. + * 11. Now you have complete implementation of new sparse matrix format. You may turn it into new type of segments (\ref TNL::Algorithms::Segments). + * + * During the implementation some unit tests may crash. If you do not need them at the moment, you may comment them in files + * `src/UnitTests/Matrices/SparseMatrixTests.h` and `src/UnitTests/Matrices/SparseMatrixVectorProductTests.h` + */ +template< typename Real = double, + typename Device = Devices::Host, + typename Index = int, + typename MatrixType = GeneralMatrix, + typename RealAllocator = typename Allocators::Default< Device >::template Allocator< Real >, + typename IndexAllocator = typename Allocators::Default< Device >::template Allocator< Index > > +class SparseSandboxMatrix : public Matrix< Real, Device, Index, RealAllocator > +{ + static_assert( + ! MatrixType::isSymmetric() || + ! std::is_same< Device, Devices::Cuda >::value || + ( std::is_same< Real, float >::value || std::is_same< Real, double >::value || std::is_same< Real, int >::value || std::is_same< Real, long long int >::value ), + "Given Real type is not supported by atomic operations on GPU which are necessary for symmetric operations." ); + + public: + + // Supporting types - they are not important for the user + using BaseType = Matrix< Real, Device, Index, RealAllocator >; + using ValuesVectorType = typename Matrix< Real, Device, Index, RealAllocator >::ValuesType; + using ValuesViewType = typename ValuesVectorType::ViewType; + using ConstValuesViewType = typename ValuesViewType::ConstViewType; + using ColumnsIndexesVectorType = Containers::Vector< typename TNL::copy_const< Index >::template from< Real >::type, Device, Index, IndexAllocator >; + using ColumnsIndexesViewType = typename ColumnsIndexesVectorType::ViewType; + using ConstColumnsIndexesViewType = typename ColumnsIndexesViewType::ConstViewType; + using RowsCapacitiesType = Containers::Vector< std::remove_const_t< Index >, Device, Index, IndexAllocator >; + using RowsCapacitiesView = Containers::VectorView< std::remove_const_t< Index >, Device, Index >; + using ConstRowsCapacitiesView = typename RowsCapacitiesView::ConstViewType; + + /** + * \brief Test of symmetric matrix type. + * + * \return \e true if the matrix is stored as symmetric and \e false otherwise. + */ + static constexpr bool isSymmetric() { return MatrixType::isSymmetric(); }; + + /** + * \brief Test of binary matrix type. + * + * \return \e true if the matrix is stored as binary and \e false otherwise. + */ + static constexpr bool isBinary() { return std::is_same< Real, bool >::value; }; + + /** + * \brief The type of matrix elements. + */ + using RealType = std::remove_const_t< Real >; + + /** + * \brief The device where the matrix is allocated. + */ + using DeviceType = Device; + + /** + * \brief The type used for matrix elements indexing. + */ + using IndexType = Index; + + /** + * \brief The allocator for matrix elements values. + */ + using RealAllocatorType = RealAllocator; + + /** + * \brief The allocator for matrix elements column indexes. + */ + using IndexAllocatorType = IndexAllocator; + + /** + * \brief Type of related matrix view. + * + * See \ref SparseSandboxMatrixView. + */ + using ViewType = SparseSandboxMatrixView< Real, Device, Index, MatrixType >; + + /** + * \brief Matrix view type for constant instances. + * + * See \ref SparseSandboxMatrixView. + */ + using ConstViewType = SparseSandboxMatrixView< std::add_const_t< Real >, Device, Index, MatrixType >; + + /** + * \brief Type for accessing matrix rows. + */ + using RowView = SparseSandboxMatrixRowView< ValuesViewType, ColumnsIndexesViewType, isBinary() >; + + /** + * \brief Type for accessing constant matrix rows. + */ + using ConstRowView = SparseSandboxMatrixRowView< ConstValuesViewType, ConstColumnsIndexesViewType, isBinary() >;; + + /** + * \brief Helper type for getting self type or its modifications. + */ + template< typename _Real = Real, + typename _Device = Device, + typename _Index = Index, + typename _MatrixType = MatrixType, + typename _RealAllocator = typename Allocators::Default< _Device >::template Allocator< _Real >, + typename _IndexAllocator = typename Allocators::Default< _Device >::template Allocator< _Index > > + using Self = SparseSandboxMatrix< _Real, _Device, _Index, _MatrixType, _RealAllocator, _IndexAllocator >; + + /** + * \brief Type of container for CSR row pointers. + * + * SANDBOX_TODO: You may replace it with containers for metadata of your format. + */ + using RowPointers = TNL::Containers::Vector< IndexType, DeviceType, IndexType >; + + /** + * \brief Constructor only with values and column indexes allocators. + * + * \param realAllocator is used for allocation of matrix elements values. + * \param indexAllocator is used for allocation of matrix elements column indexes. + */ + SparseSandboxMatrix( const RealAllocatorType& realAllocator = RealAllocatorType(), + const IndexAllocatorType& indexAllocator = IndexAllocatorType() ); + + /** + * \brief Copy constructor. + * + * \param matrix is the source matrix + */ + SparseSandboxMatrix( const SparseSandboxMatrix& matrix1 ) = default; + + /** + * \brief Move constructor. + * + * \param matrix is the source matrix + */ + SparseSandboxMatrix( SparseSandboxMatrix&& matrix ) = default; + + /** + * \brief Constructor with matrix dimensions. + * + * \param rows is number of matrix rows. + * \param columns is number of matrix columns. + * \param realAllocator is used for allocation of matrix elements values. + * \param indexAllocator is used for allocation of matrix elements column indexes. + */ + template< typename Index_t, std::enable_if_t< std::is_integral< Index_t >::value, int > = 0 > + SparseSandboxMatrix( const Index_t rows, + const Index_t columns, + const RealAllocatorType& realAllocator = RealAllocatorType(), + const IndexAllocatorType& indexAllocator = IndexAllocatorType() ); + + /** + * \brief Constructor with matrix rows capacities and number of columns. + * + * The number of matrix rows is given by the size of \e rowCapacities list. + * + * \tparam ListIndex is the initializer list values type. + * \param rowCapacities is a list telling how many matrix elements must be + * allocated in each row. + * \param columns is the number of matrix columns. + * \param realAllocator is used for allocation of matrix elements values. + * \param indexAllocator is used for allocation of matrix elements column indexes. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_Constructor_init_list_1.cpp + * \par Output + * \include SparseMatrixExample_Constructor_init_list_1.out + */ + template< typename ListIndex > + explicit SparseSandboxMatrix( const std::initializer_list< ListIndex >& rowCapacities, + const IndexType columns, + const RealAllocatorType& realAllocator = RealAllocatorType(), + const IndexAllocatorType& indexAllocator = IndexAllocatorType() ); + + /** + * \brief Constructor with matrix rows capacities given as a vector and number of columns. + * + * The number of matrix rows is given by the size of \e rowCapacities vector. + * + * \tparam RowCapacitiesVector is the row capacities vector type. Usually it is some of + * \ref TNL::Containers::Array, \ref TNL::Containers::ArrayView, \ref TNL::Containers::Vector or + * \ref TNL::Containers::VectorView. + * \param rowCapacities is a vector telling how many matrix elements must be + * allocated in each row. + * \param columns is the number of matrix columns. + * \param realAllocator is used for allocation of matrix elements values. + * \param indexAllocator is used for allocation of matrix elements column indexes. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_Constructor_rowCapacities_vector.cpp + * \par Output + * \include SparseMatrixExample_Constructor_rowCapacities_vector.out + */ + template< typename RowCapacitiesVector, std::enable_if_t< TNL::IsArrayType< RowCapacitiesVector >::value, int > = 0 > + explicit SparseSandboxMatrix( const RowCapacitiesVector& rowCapacities, + const IndexType columns, + const RealAllocatorType& realAllocator = RealAllocatorType(), + const IndexAllocatorType& indexAllocator = IndexAllocatorType() ); + + /** + * \brief Constructor with matrix dimensions and data in initializer list. + * + * The matrix elements values are given as a list \e data of triples: + * { { row1, column1, value1 }, + * { row2, column2, value2 }, + * ... }. + * + * \param rows is number of matrix rows. + * \param columns is number of matrix columns. + * \param data is a list of matrix elements values. + * \param realAllocator is used for allocation of matrix elements values. + * \param indexAllocator is used for allocation of matrix elements column indexes. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_Constructor_init_list_2.cpp + * \par Output + * \include SparseMatrixExample_Constructor_init_list_2.out + */ + explicit SparseSandboxMatrix( const IndexType rows, + const IndexType columns, + const std::initializer_list< std::tuple< IndexType, IndexType, RealType > >& data, + const RealAllocatorType& realAllocator = RealAllocatorType(), + const IndexAllocatorType& indexAllocator = IndexAllocatorType() ); + + /** + * \brief Constructor with matrix dimensions and data in std::map. + * + * The matrix elements values are given as a map \e data where keys are + * std::pair of matrix coordinates ( {row, column} ) and value is the + * matrix element value. + * + * \tparam MapIndex is a type for indexing rows and columns. + * \tparam MapValue is a type for matrix elements values in the map. + * + * \param rows is number of matrix rows. + * \param columns is number of matrix columns. + * \param map is std::map containing matrix elements. + * \param realAllocator is used for allocation of matrix elements values. + * \param indexAllocator is used for allocation of matrix elements column indexes. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_Constructor_std_map.cpp + * \par Output + * \include SparseMatrixExample_Constructor_std_map.out + */ + template< typename MapIndex, + typename MapValue > + explicit SparseSandboxMatrix( const IndexType rows, + const IndexType columns, + const std::map< std::pair< MapIndex, MapIndex >, MapValue >& map, + const RealAllocatorType& realAllocator = RealAllocatorType(), + const IndexAllocatorType& indexAllocator = IndexAllocatorType() ); + + /** + * \brief Returns a modifiable view of the sparse matrix. + * + * See \ref SparseSandboxMatrixView. + * + * \return sparse matrix view. + */ + ViewType getView() const; // TODO: remove const + + /** + * \brief Returns a non-modifiable view of the sparse matrix. + * + * See \ref SparseSandboxMatrixView. + * + * \return sparse matrix view. + */ + ConstViewType getConstView() const; + + /** + * \brief Returns string with serialization type. + * + * The string has a form `Matrices::SparseSandboxMatrix< RealType, [any_device], IndexType, General/Symmetric, Format, [any_allocator] >`. + * + * \return \ref String with the serialization type. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_getSerializationType.cpp + * \par Output + * \include SparseMatrixExample_getSerializationType.out + */ + static String getSerializationType(); + + /** + * \brief Returns string with serialization type. + * + * See \ref SparseSandboxMatrix::getSerializationType. + * + * \return \e String with the serialization type. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_getSerializationType.cpp + * \par Output + * \include SparseMatrixExample_getSerializationType.out + */ + virtual String getSerializationTypeVirtual() const override; + + /** + * \brief Set number of rows and columns of this matrix. + * + * \param rows is the number of matrix rows. + * \param columns is the number of matrix columns. + */ + virtual void setDimensions( const IndexType rows, + const IndexType columns ) override; + + /** + * \brief Set the number of matrix rows and columns by the given matrix. + * + * \tparam Matrix is matrix type. This can be any matrix having methods + * \ref getRows and \ref getColumns. + * + * \param matrix in the input matrix dimensions of which are to be adopted. + */ + template< typename Matrix > + void setLike( const Matrix& matrix ); + + /** + * \brief Allocates memory for non-zero matrix elements. + * + * The size of the input vector must be equal to the number of matrix rows. + * The number of allocated matrix elements for each matrix row depends on + * the sparse matrix format. Some formats may allocate more elements than + * required. + * + * \tparam RowsCapacitiesVector is a type of vector/array used for row + * capacities setting. + * + * \param rowCapacities is a vector telling the number of required non-zero + * matrix elements in each row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_setRowCapacities.cpp + * \par Output + * \include SparseMatrixExample_setRowCapacities.out + */ + template< typename RowsCapacitiesVector > + void setRowCapacities( const RowsCapacitiesVector& rowCapacities ); + + /** + * \brief Compute capacities of all rows. + * + * The row capacities are not stored explicitly and must be computed. + * + * \param rowCapacities is a vector where the row capacities will be stored. + */ + template< typename Vector > + void getRowCapacities( Vector& rowCapacities ) const; + + /** + * \brief This method sets the sparse matrix elements from initializer list. + * + * The number of matrix rows and columns must be set already. + * The matrix elements values are given as a list \e data of triples: + * { { row1, column1, value1 }, + * { row2, column2, value2 }, + * ... }. + * + * \param data is a initializer list of initializer lists representing + * list of matrix rows. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_setElements.cpp + * \par Output + * \include SparseMatrixExample_setElements.out + */ + void setElements( const std::initializer_list< std::tuple< IndexType, IndexType, RealType > >& data ); + + /** + * \brief This method sets the sparse matrix elements from std::map. + * + * The matrix elements values are given as a map \e data where keys are + * std::pair of matrix coordinates ( {row, column} ) and value is the + * matrix element value. + * + * \tparam MapIndex is a type for indexing rows and columns. + * \tparam MapValue is a type for matrix elements values in the map. + * + * \param map is std::map containing matrix elements. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_setElements_map.cpp + * \par Output + * \include SparseMatrixExample_setElements_map.out + */ + template< typename MapIndex, + typename MapValue > + void setElements( const std::map< std::pair< MapIndex, MapIndex > , MapValue >& map ); + + /** + * \brief Computes number of non-zeros in each row. + * + * \param rowLengths is a vector into which the number of non-zeros in each row + * will be stored. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_getCompressedRowLengths.cpp + * \par Output + * \include SparseMatrixExample_getCompressedRowLengths.out + */ + template< typename Vector > + void getCompressedRowLengths( Vector& rowLengths ) const; + + + /** + * \brief Returns capacity of given matrix row. + * + * \param row index of matrix row. + * \return number of matrix elements allocated for the row. + */ + __cuda_callable__ + IndexType getRowCapacity( const IndexType row ) const; + + /** + * \brief Returns number of non-zero matrix elements. + * + * This method really counts the non-zero matrix elements and so + * it returns zero for matrix having all allocated elements set to zero. + * + * \return number of non-zero matrix elements. + */ + IndexType getNonzeroElementsCount() const; + + /** + * \brief Resets the matrix to zero dimensions. + */ + void reset(); + + /** + * \brief Constant getter of simple structure for accessing given matrix row. + * + * \param rowIdx is matrix row index. + * + * \return RowView for accessing given matrix row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_getConstRow.cpp + * \par Output + * \include SparseMatrixExample_getConstRow.out + * + * See \ref SparseMatrixRowView. + */ + __cuda_callable__ + const ConstRowView getRow( const IndexType& rowIdx ) const; + + /** + * \brief Non-constant getter of simple structure for accessing given matrix row. + * + * \param rowIdx is matrix row index. + * + * \return RowView for accessing given matrix row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_getRow.cpp + * \par Output + * \include SparseMatrixExample_getRow.out + * + * See \ref SparseMatrixRowView. + */ + __cuda_callable__ + RowView getRow( const IndexType& rowIdx ); + + /** + * \brief Sets element at given \e row and \e column to given \e value. + * + * This method can be called from the host system (CPU) no matter + * where the matrix is allocated. If the matrix is allocated on GPU this method + * can be called even from device kernels. If the matrix is allocated in GPU device + * this method is called from CPU, it transfers values of each matrix element separately and so the + * performance is very low. For higher performance see. \ref SparseMatrix::getRow + * or \ref SparseMatrix::forElements and \ref SparseMatrix::forAllElements. + * The call may fail if the matrix row capacity is exhausted. + * + * \param row is row index of the element. + * \param column is columns index of the element. + * \param value is the value the element will be set to. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_setElement.cpp + * \par Output + * \include SparseMatrixExample_setElement.out + */ + __cuda_callable__ + void setElement( const IndexType row, + const IndexType column, + const RealType& value ); + + /** + * \brief Add element at given \e row and \e column to given \e value. + * + * This method can be called from the host system (CPU) no matter + * where the matrix is allocated. If the matrix is allocated on GPU this method + * can be called even from device kernels. If the matrix is allocated in GPU device + * this method is called from CPU, it transfers values of each matrix element separately and so the + * performance is very low. For higher performance see. \ref SparseMatrix::getRow + * or \ref SparseMatrix::forElements and \ref SparseMatrix::forAllElements. + * The call may fail if the matrix row capacity is exhausted. + * + * \param row is row index of the element. + * \param column is columns index of the element. + * \param value is the value the element will be set to. + * \param thisElementMultiplicator is multiplicator the original matrix element + * value is multiplied by before addition of given \e value. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_addElement.cpp + * \par Output + * \include SparseMatrixExample_addElement.out + * + */ + __cuda_callable__ + void addElement( const IndexType row, + const IndexType column, + const RealType& value, + const RealType& thisElementMultiplicator ); + + /** + * \brief Returns value of matrix element at position given by its row and column index. + * + * This method can be called from the host system (CPU) no matter + * where the matrix is allocated. If the matrix is allocated on GPU this method + * can be called even from device kernels. If the matrix is allocated in GPU device + * this method is called from CPU, it transfers values of each matrix element separately and so the + * performance is very low. For higher performance see. \ref SparseMatrix::getRow + * or \ref SparseMatrix::forElements and \ref SparseMatrix::forAllElements. + * + * \param row is a row index of the matrix element. + * \param column i a column index of the matrix element. + * + * \return value of given matrix element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_getElement.cpp + * \par Output + * \include SparseMatrixExample_getElement.out + * + */ + __cuda_callable__ + RealType getElement( const IndexType row, + const IndexType column ) const; + + /** + * \brief Method for performing general reduction on matrix rows. + * + * \tparam Fetch is a type of lambda function for data fetch declared as + * `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`. + * The return type of this lambda can be any non void. + * \tparam Reduce is a type of lambda function for reduction declared as + * `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`. + * \tparam Keep is a type of lambda function for storing results of reduction in each row. + * It is declared as `keep( const IndexType rowIdx, const double& value )`. + * \tparam FetchValue is type returned by the Fetch lambda function. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param fetch is an instance of lambda function for data fetch. + * \param reduce is an instance of lambda function for reduction. + * \param keep in an instance of lambda function for storing results. + * \param zero is zero of given reduction operation also known as idempotent element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_reduceRows.cpp + * \par Output + * \include SparseMatrixExample_reduceRows.out + */ + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > + void reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ); + + /** + * \brief Method for performing general reduction on matrix rows for constant instances. + * + * \tparam Fetch is a type of lambda function for data fetch declared as + * `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`. + * The return type of this lambda can be any non void. + * \tparam Reduce is a type of lambda function for reduction declared as + * `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`. + * \tparam Keep is a type of lambda function for storing results of reduction in each row. + * It is declared as `keep( const IndexType rowIdx, const double& value )`. + * \tparam FetchValue is type returned by the Fetch lambda function. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param fetch is an instance of lambda function for data fetch. + * \param reduce is an instance of lambda function for reduction. + * \param keep in an instance of lambda function for storing results. + * \param zero is zero of given reduction operation also known as idempotent element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_reduceRows.cpp + * \par Output + * \include SparseMatrixExample_reduceRows.out + */ + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > + void reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const; + + /** + * \brief Method for performing general reduction on all matrix rows. + * + * \tparam Fetch is a type of lambda function for data fetch declared as + * `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`. + * The return type of this lambda can be any non void. + * \tparam Reduce is a type of lambda function for reduction declared as + * `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`. + * \tparam Keep is a type of lambda function for storing results of reduction in each row. + * It is declared as `keep( const IndexType rowIdx, const double& value )`. + * \tparam FetchValue is type returned by the Fetch lambda function. + * + * \param fetch is an instance of lambda function for data fetch. + * \param reduce is an instance of lambda function for reduction. + * \param keep in an instance of lambda function for storing results. + * \param zero is zero of given reduction operation also known as idempotent element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_reduceAllRows.cpp + * \par Output + * \include SparseMatrixExample_reduceAllRows.out + */ + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > + void reduceAllRows( Fetch&& fetch, const Reduce&& reduce, Keep&& keep, const FetchReal& zero ); + + /** + * \brief Method for performing general reduction on all matrix rows for constant instances. + * + * \tparam Fetch is a type of lambda function for data fetch declared as + * `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`. + * The return type of this lambda can be any non void. + * \tparam Reduce is a type of lambda function for reduction declared as + * `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`. + * \tparam Keep is a type of lambda function for storing results of reduction in each row. + * It is declared as `keep( const IndexType rowIdx, const double& value )`. + * \tparam FetchValue is type returned by the Fetch lambda function. + * + * \param fetch is an instance of lambda function for data fetch. + * \param reduce is an instance of lambda function for reduction. + * \param keep in an instance of lambda function for storing results. + * \param zero is zero of given reduction operation also known as idempotent element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_reduceAllRows.cpp + * \par Output + * \include SparseMatrixExample_reduceAllRows.out + */ + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > + void reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const; + + /** + * \brief Method for parallel iteration over matrix elements of given rows for constant instances. + * + * \tparam Function is type of lambda function that will operate on matrix elements. + * + * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed. + * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed. + * \param function is an instance of the lambda function to be called for element of given rows. + * + * The lambda function `function` should be declared like follows: + * + * ``` + * auto function = [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value ) { ... }; + * ``` + * + * The \e localIdx parameter is a rank of the non-zero element in given row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp + * \par Output + * \include SparseMatrixExample_forElements.out + */ + template< typename Function > + void forElements( IndexType begin, IndexType end, Function&& function ) const; + + /** + * \brief Method for parallel iteration over all matrix elements of given rows for non-constant instances. + * + * \tparam Function is type of lambda function that will operate on matrix elements. + * + * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed. + * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed. + * \param function is an instance of the lambda function to be called for each element of given rows. + * + * The lambda function `function` should be declared like follows: + * + * ``` + * auto function = [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value ) mutable { ... } + * ``` + * + * The \e localIdx parameter is a rank of the non-zero element in given row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp + * \par Output + * \include SparseMatrixExample_forElements.out + */ + template< typename Function > + void forElements( IndexType begin, IndexType end, Function&& function ); + + /** + * \brief Method for parallel iteration over all matrix elements for constant instances. + * + * See \ref SparseMatrix::forElements. + * + * \tparam Function is a type of lambda function that will operate on matrix elements. + * \param function is an instance of the lambda function to be called for each matrix element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp + * \par Output + * \include SparseMatrixExample_forElements.out + */ + template< typename Function > + void forAllElements( Function&& function ) const; + + /** + * \brief Method for parallel iteration over all matrix elements for non-constant instances. + * + * See \ref SparseMatrix::forElements. + * + * \tparam Function is a type of lambda function that will operate on matrix elements. + * \param function is an instance of the lambda function to be called for each matrix element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp + * \par Output + * \include SparseMatrixExample_forElements.out + */ + template< typename Function > + void forAllElements( Function&& function ); + + /** + * \brief Method for parallel iteration over matrix rows from interval [ \e begin, \e end). + * + * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method + * \ref SparseMatrix::forElements where more than one thread can be mapped to each row. + * + * \tparam Function is type of the lambda function. + * + * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed. + * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed. + * \param function is an instance of the lambda function to be called for each row. + * + * ``` + * auto function = [] __cuda_callable__ ( RowView& row ) mutable { ... }; + * ``` + * + * \e RowView represents matrix row - see \ref TNL::Matrices::SparseMatrix::RowView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_forRows.cpp + * \par Output + * \include SparseMatrixExample_forRows.out + */ + template< typename Function > + void forRows( IndexType begin, IndexType end, Function&& function ); + + /** + * \brief Method for parallel iteration over matrix rows from interval [ \e begin, \e end) for constant instances. + * + * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method + * \ref SparseMatrix::forElements where more than one thread can be mapped to each row. + * + * \tparam Function is type of the lambda function. + * + * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed. + * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed. + * \param function is an instance of the lambda function to be called for each row. + * + * ``` + * auto function = [] __cuda_callable__ ( RowView& row ) { ... }; + * ``` + * + * \e RowView represents matrix row - see \ref TNL::Matrices::SparseMatrix::RowView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_forRows.cpp + * \par Output + * \include SparseMatrixExample_forRows.out + */ + template< typename Function > + void forRows( IndexType begin, IndexType end, Function&& function ) const; + + /** + * \brief Method for parallel iteration over all matrix rows. + * + * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method + * \ref SparseMatrix::forAllElements where more than one thread can be mapped to each row. + * + * \tparam Function is type of the lambda function. + * + * \param function is an instance of the lambda function to be called for each row. + * + * ``` + * auto function = [] __cuda_callable__ ( RowView& row ) mutable { ... }; + * ``` + * + * \e RowView represents matrix row - see \ref TNL::Matrices::SparseMatrix::RowView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_forRows.cpp + * \par Output + * \include SparseMatrixExample_forRows.out + */ + template< typename Function > + void forAllRows( Function&& function ); + + /** + * \brief Method for parallel iteration over all matrix rows for constant instances. + * + * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method + * \ref SparseMatrix::forAllElements where more than one thread can be mapped to each row. + * + * \tparam Function is type of the lambda function. + * + * \param function is an instance of the lambda function to be called for each row. + * + * ``` + * auto function = [] __cuda_callable__ ( RowView& row ) { ... }; + * ``` + * + * \e RowView represents matrix row - see \ref TNL::Matrices::SparseMatrix::RowView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_forRows.cpp + * \par Output + * \include SparseMatrixExample_forRows.out + */ + template< typename Function > + void forAllRows( Function&& function ) const; + + /** + * \brief Method for sequential iteration over all matrix rows for constant instances. + * + * \tparam Function is type of lambda function that will operate on matrix elements. + * It is should have form like + * `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`. + * The column index repeats twice only for compatibility with sparse matrices. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param function is an instance of the lambda function to be called in each row. + */ + template< typename Function > + void sequentialForRows( IndexType begin, IndexType end, Function& function ) const; + + /** + * \brief Method for sequential iteration over all matrix rows for non-constant instances. + * + * \tparam Function is type of lambda function that will operate on matrix elements. + * It is should have form like + * `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`. + * The column index repeats twice only for compatibility with sparse matrices. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param function is an instance of the lambda function to be called in each row. + */ + template< typename Function > + void sequentialForRows( IndexType begin, IndexType end, Function& function ); + + /** + * \brief This method calls \e sequentialForRows for all matrix rows (for constant instances). + * + * See \ref SparseMatrix::sequentialForRows. + * + * \tparam Function is a type of lambda function that will operate on matrix elements. + * \param function is an instance of the lambda function to be called in each row. + */ + template< typename Function > + void sequentialForAllRows( Function& function ) const; + + /** + * \brief This method calls \e sequentialForRows for all matrix rows. + * + * See \ref SparseMatrix::sequentialForAllRows. + * + * \tparam Function is a type of lambda function that will operate on matrix elements. + * \param function is an instance of the lambda function to be called in each row. + */ + template< typename Function > + void sequentialForAllRows( Function& function ); + + /** + * \brief Computes product of matrix and vector. + * + * More precisely, it computes: + * + * `outVector = matrixMultiplicator * ( * this ) * inVector + outVectorMultiplicator * outVector` + * + * \tparam InVector is type of input vector. It can be \ref Vector, + * \ref VectorView, \ref Array, \ref ArraView or similar container. + * \tparam OutVector is type of output vector. It can be \ref Vector, + * \ref VectorView, \ref Array, \ref ArraView or similar container. + * + * \param inVector is input vector. + * \param outVector is output vector. + * \param matrixMultiplicator is a factor by which the matrix is multiplied. It is one by default. + * \param outVectorMultiplicator is a factor by which the outVector is multiplied before added + * to the result of matrix-vector product. It is zero by default. + * \param begin is the beginning of the rows range for which the vector product + * is computed. It is zero by default. + * \param end is the end of the rows range for which the vector product + * is computed. It is number if the matrix rows by default. + */ + template< typename InVector, + typename OutVector > + void vectorProduct( const InVector& inVector, + OutVector& outVector, + const RealType& matrixMultiplicator = 1.0, + const RealType& outVectorMultiplicator = 0.0, + const IndexType firstRow = 0, + const IndexType lastRow = 0 ) const; + + /*template< typename Real2, typename Index2 > + void addMatrix( const SparseMatrix< Real2, Segments, Device, Index2 >& matrix, + const RealType& matrixMultiplicator = 1.0, + const RealType& thisMatrixMultiplicator = 1.0 ); + + template< typename Real2, typename Index2 > + void getTransposition( const SparseMatrix< Real2, Segments, Device, Index2 >& matrix, + const RealType& matrixMultiplicator = 1.0 ); + */ + + template< typename Vector1, typename Vector2 > + bool performSORIteration( const Vector1& b, + const IndexType row, + Vector2& x, + const RealType& omega = 1.0 ) const; + + /** + * \brief Assignment of exactly the same matrix type. + * + * \param matrix is input matrix for the assignment. + * \return reference to this matrix. + */ + SparseSandboxMatrix& operator=( const SparseSandboxMatrix& matrix ); + + /** + * \brief Assignment of exactly the same matrix type but different device. + * + * \param matrix is input matrix for the assignment. + * \return reference to this matrix. + */ + template< typename Device_ > + SparseSandboxMatrix& operator=( const SparseSandboxMatrix< RealType, Device_, IndexType, MatrixType, RealAllocator, IndexAllocator >& matrix ); + + /** + * \brief Assignment of dense matrix + * + * \param matrix is input matrix for the assignment. + * \return reference to this matrix. + */ + template< typename Real_, typename Device_, typename Index_, ElementsOrganization Organization, typename RealAllocator_ > + SparseSandboxMatrix& operator=( const DenseMatrix< Real_, Device_, Index_, Organization, RealAllocator_ >& matrix ); + + + /** + * \brief Assignment of any matrix type other then this and dense. + * + * **Warning: Assignment of symmetric sparse matrix to general sparse matrix does not give correct result, currently. Only the diagonal and the lower part of the matrix is assigned.** + * + * \param matrix is input matrix for the assignment. + * \return reference to this matrix. + */ + template< typename RHSMatrix > + SparseSandboxMatrix& operator=( const RHSMatrix& matrix ); + + /** + * \brief Comparison operator with another arbitrary matrix type. + * + * \param matrix is the right-hand side matrix. + * \return \e true if the RHS matrix is equal, \e false otherwise. + */ + template< typename Matrix > + bool operator==( const Matrix& m ) const; + + /** + * \brief Comparison operator with another arbitrary matrix type. + * + * \param matrix is the right-hand side matrix. + * \return \e true if the RHS matrix is equal, \e false otherwise. + */ + template< typename Matrix > + bool operator!=( const Matrix& m ) const; + + /** + * \brief Method for saving the matrix to the file with given filename. + * + * \param fileName is name of the file. + */ + void save( const String& fileName ) const; + + /** + * \brief Method for loading the matrix from the file with given filename. + * + * \param fileName is name of the file. + */ + void load( const String& fileName ); + + /** + * \brief Method for saving the matrix to a file. + * + * \param file is the output file. + */ + virtual void save( File& file ) const override; + + /** + * \brief Method for loading the matrix from a file. + * + * \param file is the input file. + */ + virtual void load( File& file ) override; + + /** + * \brief Method for printing the matrix to output stream. + * + * \param str is the output stream. + */ + virtual void print( std::ostream& str ) const override; + + /** + * \brief Returns a padding index value. + * + * Padding index is used for column indexes of padding zeros. Padding zeros + * are used in some sparse matrix formats for better data alignment in memory. + * + * \return value of the padding index. + */ + __cuda_callable__ + IndexType getPaddingIndex() const; + + /** + * \brief Getter of segments for non-constant instances. + * + * \e Segments are a structure for addressing the matrix elements columns and values. + * In fact, \e Segments represent the sparse matrix format. + * + * \return Non-constant reference to segments. + */ + //SegmentsType& getSegments(); + + /** + * \brief Getter of segments for constant instances. + * + * \e Segments are a structure for addressing the matrix elements columns and values. + * In fact, \e Segments represent the sparse matrix format. + * + * \return Constant reference to segments. + */ + //const SegmentsType& getSegments() const; + + /** + * \brief Getter of column indexes for constant instances. + * + * \return Constant reference to a vector with matrix elements column indexes. + */ + const ColumnsIndexesVectorType& getColumnIndexes() const; + + /** + * \brief Getter of column indexes for nonconstant instances. + * + * \return Reference to a vector with matrix elements column indexes. + */ + ColumnsIndexesVectorType& getColumnIndexes(); + + protected: + + ColumnsIndexesVectorType columnIndexes; + + IndexAllocator indexAllocator; + + ViewType view; + + /** + * \brief Container for CSR row pointers. + * + * SANDBOX_TODO: You may replace it with containers and metadata required by you format. + */ + + RowPointers rowPointers; +}; + + } // namespace Sandbox + } // namespace Matrices +} // namespace TNL + +#include <TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp> diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp new file mode 100644 index 0000000000000000000000000000000000000000..63f49e6c84746add2b368c5381e307d21b06eb0b --- /dev/null +++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp @@ -0,0 +1,1197 @@ +/*************************************************************************** + SparseSandboxMatrix.hpp - description + ------------------- + begin : Apr 19, 2021 + copyright : (C) 2021 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +#include <functional> +#include <sstream> +#include <TNL/Algorithms/Reduction.h> +#include <TNL/Matrices/Sandbox/SparseSandboxMatrix.h> + +namespace TNL { + namespace Matrices { + namespace Sandbox { + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +SparseSandboxMatrix( const RealAllocatorType& realAllocator, + const IndexAllocatorType& indexAllocator ) +: BaseType( realAllocator ), columnIndexes( indexAllocator ), rowPointers( ( IndexType ) 1, ( IndexType ) 0, indexAllocator ) +{ + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Index_t, std::enable_if_t< std::is_integral< Index_t >::value, int > > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +SparseSandboxMatrix( const Index_t rows, + const Index_t columns, + const RealAllocatorType& realAllocator, + const IndexAllocatorType& indexAllocator ) +: BaseType( rows, columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rows + 1, ( IndexType ) 0, indexAllocator ) +{ + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename ListIndex > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +SparseSandboxMatrix( const std::initializer_list< ListIndex >& rowCapacities, + const IndexType columns, + const RealAllocatorType& realAllocator, + const IndexAllocatorType& indexAllocator ) +: BaseType( rowCapacities.size(), columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rowCapacities.size() + 1, ( IndexType ) 0, indexAllocator ) +{ + this->setRowCapacities( RowsCapacitiesType( rowCapacities ) ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename RowCapacitiesVector, std::enable_if_t< TNL::IsArrayType< RowCapacitiesVector >::value, int > > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +SparseSandboxMatrix( const RowCapacitiesVector& rowCapacities, + const IndexType columns, + const RealAllocatorType& realAllocator, + const IndexAllocatorType& indexAllocator ) +: BaseType( rowCapacities.getSize(), columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rowCapacities.getSize() + 1, ( IndexType ) 0, indexAllocator ) +{ + this->setRowCapacities( rowCapacities ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +SparseSandboxMatrix( const IndexType rows, + const IndexType columns, + const std::initializer_list< std::tuple< IndexType, IndexType, RealType > >& data, + const RealAllocatorType& realAllocator, + const IndexAllocatorType& indexAllocator ) +: BaseType( rows, columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rows + 1, ( IndexType ) 0, indexAllocator ) +{ + this->setElements( data ); + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename MapIndex, + typename MapValue > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +SparseSandboxMatrix( const IndexType rows, + const IndexType columns, + const std::map< std::pair< MapIndex, MapIndex > , MapValue >& map, + const RealAllocatorType& realAllocator, + const IndexAllocatorType& indexAllocator ) +: BaseType( rows, columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rows + 1, ( IndexType ) 0, indexAllocator ) +{ + this->setDimensions( rows, columns ); + this->setElements( map ); + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +auto +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getView() const -> ViewType +{ + return ViewType( this->getRows(), + this->getColumns(), + const_cast< SparseSandboxMatrix* >( this )->getValues().getView(), // TODO: remove const_cast + const_cast< SparseSandboxMatrix* >( this )->columnIndexes.getView(), + const_cast< SparseSandboxMatrix* >( this )->rowPointers.getView() ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +auto +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getConstView() const -> ConstViewType +{ + return ConstViewType( this->getRows(), + this->getColumns(), + this->getValues().getConstView(), + this->columnIndexes.getConstView(), + this->segments.getConstView() ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +String +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getSerializationType() +{ + return ViewType::getSerializationType(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +String +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getSerializationTypeVirtual() const +{ + return this->getSerializationType(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +setDimensions( const IndexType rows, + const IndexType columns ) +{ + BaseType::setDimensions( rows, columns ); + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Matrix_ > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +setLike( const Matrix_& matrix ) +{ + BaseType::setLike( matrix ); + // SANDBOX_TODO: Replace the following line with assignment of metadata required by your format. + // Do not assign matrix elements here. + this->rowPointers = matrix.rowPointers; + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename RowsCapacitiesVector > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +setRowCapacities( const RowsCapacitiesVector& rowsCapacities ) +{ + TNL_ASSERT_EQ( rowsCapacities.getSize(), this->getRows(), "Number of matrix rows does not fit with rowCapacities vector size." ); + using RowsCapacitiesVectorDevice = typename RowsCapacitiesVector::DeviceType; + + // SANDBOX_TODO: Replace the following lines with the setup of your sparse matrix format based on + // `rowsCapacities`. This container has the same number of elements as is the number of + // rows of this matrix. Each element says how many nonzero elements the user needs to have + // in each row. This number can be increased if the sparse matrix format uses padding zeros. + this->rowPointers.setSize( this->getRows() + 1 ); + if( std::is_same< DeviceType, RowsCapacitiesVectorDevice >::value ) + { + // GOTCHA: when this->getRows() == 0, getView returns a full view with size == 1 + if( this->getRows() > 0 ) { + auto view = this->rowPointers.getView( 0, this->getRows() ); + view = rowsCapacities; + } + } + else + { + RowsCapacitiesType thisRowsCapacities; + thisRowsCapacities = rowsCapacities; + if( this->getRows() > 0 ) { + auto view = this->rowPointers.getView( 0, this->getRows() ); + view = thisRowsCapacities; + } + } + this->rowPointers.setElement( this->getRows(), 0 ); + this->rowPointers.template scan< Algorithms::ScanType::Exclusive >(); + // End of sparse matrix format initiation. + + // SANDBOX_TODO: Compute number of all elements that need to be allocated by your format. + const auto storageSize = rowPointers.getElement( this->getRows() ); + + // The rest of this methods needs no changes. + if( ! isBinary() ) + { + this->values.setSize( storageSize ); + this->values = ( RealType ) 0; + } + this->columnIndexes.setSize( storageSize ); + this->columnIndexes = this->getPaddingIndex(); + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Vector > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getRowCapacities( Vector& rowCapacities ) const +{ + this->view.getRowCapacities( rowCapacities ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +setElements( const std::initializer_list< std::tuple< IndexType, IndexType, RealType > >& data ) +{ + const auto& rows = this->getRows(); + const auto& columns = this->getColumns(); + Containers::Vector< IndexType, Devices::Host, IndexType > rowCapacities( rows, 0 ); + for( const auto& i : data ) + { + if( std::get< 0 >( i ) >= rows ) + { + std::stringstream s; + s << "Wrong row index " << std::get< 0 >( i ) << " in an initializer list"; + throw std::logic_error( s.str() ); + } + rowCapacities[ std::get< 0 >( i ) ]++; + } + SparseSandboxMatrix< Real, Devices::Host, Index, MatrixType > hostMatrix( rows, columns ); + hostMatrix.setRowCapacities( rowCapacities ); + for( const auto& i : data ) + { + if( std::get< 1 >( i ) >= columns ) + { + std::stringstream s; + s << "Wrong column index " << std::get< 1 >( i ) << " in an initializer list"; + throw std::logic_error( s.str() ); + } + hostMatrix.setElement( std::get< 0 >( i ), std::get< 1 >( i ), std::get< 2 >( i ) ); + } + ( *this ) = hostMatrix; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename MapIndex, + typename MapValue > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +setElements( const std::map< std::pair< MapIndex, MapIndex > , MapValue >& map ) +{ + Containers::Vector< IndexType, Devices::Host, IndexType > rowsCapacities( this->getRows(), 0 ); + for( auto element : map ) + rowsCapacities[ element.first.first ]++; + if( !std::is_same< DeviceType, Devices::Host >::value ) + { + SparseSandboxMatrix< Real, Devices::Host, Index, MatrixType > hostMatrix( this->getRows(), this->getColumns() ); + hostMatrix.setRowCapacities( rowsCapacities ); + for( auto element : map ) + hostMatrix.setElement( element.first.first, element.first.second, element.second ); + *this = hostMatrix; + } + else + { + this->setRowCapacities( rowsCapacities ); + for( auto element : map ) + this->setElement( element.first.first, element.first.second, element.second ); + } +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Vector > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getCompressedRowLengths( Vector& rowLengths ) const +{ + this->view.getCompressedRowLengths( rowLengths ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +__cuda_callable__ +Index +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getRowCapacity( const IndexType row ) const +{ + return this->view.getRowCapacity( row ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +Index +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getNonzeroElementsCount() const +{ + return this->view.getNonzeroElementsCount(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +reset() +{ + BaseType::reset(); + this->columnIndexes.reset(); + // SANDBOX_TODO: Reset the metadata required by your format here. + this->rowPointers.reset(); + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +__cuda_callable__ auto +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getRow( const IndexType& rowIdx ) const -> const ConstRowView +{ + return this->view.getRow( rowIdx ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +__cuda_callable__ auto +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getRow( const IndexType& rowIdx ) -> RowView +{ + return this->view.getRow( rowIdx ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +__cuda_callable__ void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +setElement( const IndexType row, + const IndexType column, + const RealType& value ) +{ + this->view.setElement( row, column, value ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +__cuda_callable__ void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +addElement( const IndexType row, + const IndexType column, + const RealType& value, + const RealType& thisElementMultiplicator ) +{ + this->view.addElement( row, column, value, thisElementMultiplicator ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +__cuda_callable__ +auto +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getElement( const IndexType row, + const IndexType column ) const -> RealType +{ + return this->view.getElement( row, column ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +template< typename InVector, + typename OutVector > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +vectorProduct( const InVector& inVector, + OutVector& outVector, + const RealType& matrixMultiplicator, + const RealType& outVectorMultiplicator, + const IndexType firstRow, + const IndexType lastRow ) const +{ + this->view.vectorProduct( inVector, outVector, matrixMultiplicator, outVectorMultiplicator, firstRow, lastRow ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Fetch, typename Reduce, typename Keep, typename FetchValue > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchValue& zero ) +{ + this->view.reduceRows( begin, end, fetch, reduce, keep, zero ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Fetch, typename Reduce, typename Keep, typename FetchValue > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchValue& zero ) const +{ + this->view.reduceRows( begin, end, fetch, reduce, keep, zero ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +reduceAllRows( Fetch&& fetch, const Reduce&& reduce, Keep&& keep, const FetchReal& zero ) +{ + this->reduceRows( 0, this->getRows(), fetch, reduce, keep, zero ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const +{ + this->reduceRows( 0, this->getRows(), fetch, reduce, keep, zero ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +forElements( IndexType begin, IndexType end, Function&& function ) const +{ + this->view.forElements( begin, end, function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +forElements( IndexType begin, IndexType end, Function&& function ) +{ + this->view.forElements( begin, end, function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +forAllElements( Function&& function ) const +{ + this->forElements( 0, this->getRows(), function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +forAllElements( Function&& function ) +{ + this->forElements( 0, this->getRows(), function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +forRows( IndexType begin, IndexType end, Function&& function ) +{ + this->getView().forRows( begin, end, function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +forRows( IndexType begin, IndexType end, Function&& function ) const +{ + this->getConstView().forRows( begin, end, function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +forAllRows( Function&& function ) +{ + this->getView().forAllRows( function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +forAllRows( Function&& function ) const +{ + this->getConsView().forAllRows( function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +sequentialForRows( IndexType begin, IndexType end, Function& function ) const +{ + this->view.sequentialForRows( begin, end, function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +sequentialForRows( IndexType first, IndexType last, Function& function ) +{ + this->view.sequentialForRows( first, last, function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +sequentialForAllRows( Function& function ) const +{ + this->sequentialForRows( 0, this->getRows(), function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Function > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +sequentialForAllRows( Function& function ) +{ + this->sequentialForRows( 0, this->getRows(), function ); +} + + +/*template< typename Real, + template< typename, typename, typename > class Segments, + typename Device, + typename Index, + typename RealAllocator, + typename IndexAllocator > +template< typename Real2, template< typename, typename > class Segments2, typename Index2, typename RealAllocator2, typename IndexAllocator2 > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +addMatrix( const SparseSandboxMatrix< Real2, Segments2, Device, Index2, RealAllocator2, IndexAllocator2 >& matrix, + const RealType& matrixMultiplicator, + const RealType& thisMatrixMultiplicator ) +{ + +} + +template< typename Real, + template< typename, typename, typename > class Segments, + typename Device, + typename Index, + typename RealAllocator, + typename IndexAllocator > +template< typename Real2, typename Index2 > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getTransposition( const SparseSandboxMatrix< Real2, Device, Index2 >& matrix, + const RealType& matrixMultiplicator ) +{ + +}*/ + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +template< typename Vector1, typename Vector2 > +bool +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +performSORIteration( const Vector1& b, + const IndexType row, + Vector2& x, + const RealType& omega ) const +{ + return false; +} + +// copy assignment +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >& +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +operator=( const SparseSandboxMatrix& matrix ) +{ + Matrix< Real, Device, Index >::operator=( matrix ); + this->columnIndexes = matrix.columnIndexes; + // SANDBOX_TODO: Replace the following line with an assignment of metadata required by you sparse matrix format. + this->rowPointers = matrix.rowPointers; + this->view = this->getView(); + return *this; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Device_ > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >& +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +operator=( const SparseSandboxMatrix< RealType, Device_, IndexType, MatrixType, RealAllocator, IndexAllocator >& matrix ) +{ + Matrix< Real, Device, Index >::operator=( matrix ); + this->columnIndexes = matrix.columnIndexes; + // SANDBOX_TODO: Replace the following line with an assignment of metadata required by you sparse matrix format. + this->rowPointers = matrix.rowPointers; + this->view = this->getView(); + return *this; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Real_, typename Device_, typename Index_, ElementsOrganization Organization, typename RealAllocator_ > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >& +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +operator=( const DenseMatrix< Real_, Device_, Index_, Organization, RealAllocator_ >& matrix ) +{ + using RHSMatrix = DenseMatrix< Real_, Device_, Index_, Organization, RealAllocator_ >; + using RHSIndexType = typename RHSMatrix::IndexType; + using RHSRealType = typename RHSMatrix::RealType; + using RHSDeviceType = typename RHSMatrix::DeviceType; + using RHSRealAllocatorType = typename RHSMatrix::RealAllocatorType; + + Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rowLengths; + matrix.getCompressedRowLengths( rowLengths ); + this->setLike( matrix ); + this->setRowCapacities( rowLengths ); + Containers::Vector< IndexType, DeviceType, IndexType > rowLocalIndexes( matrix.getRows() ); + rowLocalIndexes = 0; + + // TODO: use getConstView when it works + const auto matrixView = const_cast< RHSMatrix& >( matrix ).getView(); + const IndexType paddingIndex = this->getPaddingIndex(); + auto columns_view = this->columnIndexes.getView(); + auto values_view = this->values.getView(); + auto rowLocalIndexes_view = rowLocalIndexes.getView(); + columns_view = paddingIndex; + + if( std::is_same< DeviceType, RHSDeviceType >::value ) + { + const auto segments_view = this->segments.getView(); + auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIdx, const RHSRealType& value ) mutable { + if( value != 0.0 ) + { + IndexType thisGlobalIdx = segments_view.getGlobalIndex( rowIdx, rowLocalIndexes_view[ rowIdx ]++ ); + columns_view[ thisGlobalIdx ] = columnIdx; + if( ! isBinary() ) + values_view[ thisGlobalIdx ] = value; + } + }; + matrix.forAllElements( f ); + } + else + { + const IndexType maxRowLength = matrix.getColumns(); + const IndexType bufferRowsCount( 128 ); + const size_t bufferSize = bufferRowsCount * maxRowLength; + Containers::Vector< RHSRealType, RHSDeviceType, RHSIndexType, RHSRealAllocatorType > matrixValuesBuffer( bufferSize ); + Containers::Vector< RealType, DeviceType, IndexType, RealAllocatorType > thisValuesBuffer( bufferSize ); + Containers::Vector< IndexType, DeviceType, IndexType, IndexAllocatorType > thisColumnsBuffer( bufferSize ); + auto matrixValuesBuffer_view = matrixValuesBuffer.getView(); + auto thisValuesBuffer_view = thisValuesBuffer.getView(); + + IndexType baseRow( 0 ); + const IndexType rowsCount = this->getRows(); + while( baseRow < rowsCount ) + { + const IndexType lastRow = min( baseRow + bufferRowsCount, rowsCount ); + thisColumnsBuffer = paddingIndex; + + //// + // Copy matrix elements into buffer + auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value ) mutable { + const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx; + matrixValuesBuffer_view[ bufferIdx ] = value; + }; + matrix.forElements( baseRow, lastRow, f1 ); + + //// + // Copy the source matrix buffer to this matrix buffer + thisValuesBuffer_view = matrixValuesBuffer_view; + + //// + // Copy matrix elements from the buffer to the matrix and ignoring + // zero matrix elements. + const IndexType matrix_columns = this->getColumns(); + auto f2 = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIndex, RealType& value ) mutable { + RealType inValue( 0.0 ); + IndexType bufferIdx, column( rowLocalIndexes_view[ rowIdx ] ); + while( inValue == 0.0 && column < matrix_columns ) + { + bufferIdx = ( rowIdx - baseRow ) * maxRowLength + column++; + inValue = thisValuesBuffer_view[ bufferIdx ]; + } + rowLocalIndexes_view[ rowIdx ] = column; + if( inValue == 0.0 ) + { + columnIndex = paddingIndex; + value = 0.0; + } + else + { + columnIndex = column - 1; + value = inValue; + } + }; + this->forElements( baseRow, lastRow, f2 ); + baseRow += bufferRowsCount; + } + //std::cerr << "This matrix = " << std::endl << *this << std::endl; + } + this->view = this->getView(); + return *this; + +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename RHSMatrix > +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >& +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +operator=( const RHSMatrix& matrix ) +{ + using RHSIndexType = typename RHSMatrix::IndexType; + using RHSRealType = typename RHSMatrix::RealType; + using RHSDeviceType = typename RHSMatrix::DeviceType; + using RHSRealAllocatorType = typename RHSMatrix::RealAllocatorType; + + Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rowCapacities; + matrix.getRowCapacities( rowCapacities ); + this->setDimensions( matrix.getRows(), matrix.getColumns() ); + this->setRowCapacities( rowCapacities ); + Containers::Vector< IndexType, DeviceType, IndexType > rowLocalIndexes( matrix.getRows() ); + rowLocalIndexes = 0; + + // TODO: use getConstView when it works + const auto matrixView = const_cast< RHSMatrix& >( matrix ).getView(); + const IndexType paddingIndex = this->getPaddingIndex(); + auto columns_view = this->columnIndexes.getView(); + auto values_view = this->values.getView(); + auto rowLocalIndexes_view = rowLocalIndexes.getView(); + columns_view = paddingIndex; + + // SANDBOX_TODO: Modify the follwoing accoring to your format + auto row_pointers_view = this->rowPointers.getView(); + if( std::is_same< DeviceType, RHSDeviceType >::value ) + { + auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx_, RHSIndexType columnIndex, const RHSRealType& value ) mutable { + IndexType localIdx( rowLocalIndexes_view[ rowIdx ] ); + IndexType thisRowBegin = row_pointers_view[ rowIdx ]; + if( value != 0.0 && columnIndex != paddingIndex ) + { + IndexType thisGlobalIdx = thisRowBegin + localIdx++; + columns_view[ thisGlobalIdx ] = columnIndex; + if( ! isBinary() ) + values_view[ thisGlobalIdx ] = value; + rowLocalIndexes_view[ rowIdx ] = localIdx; + } + }; + matrix.forAllElements( f ); + } + else + { + const IndexType maxRowLength = max( rowCapacities ); + const IndexType bufferRowsCount( 128 ); + const size_t bufferSize = bufferRowsCount * maxRowLength; + Containers::Vector< RHSRealType, RHSDeviceType, RHSIndexType, RHSRealAllocatorType > matrixValuesBuffer( bufferSize ); + Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > matrixColumnsBuffer( bufferSize ); + Containers::Vector< RealType, DeviceType, IndexType, RealAllocatorType > thisValuesBuffer( bufferSize ); + Containers::Vector< IndexType, DeviceType, IndexType > thisColumnsBuffer( bufferSize ); + Containers::Vector< IndexType, DeviceType, IndexType > thisRowLengths; + Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rhsRowLengths; + matrix.getCompressedRowLengths( rhsRowLengths ); + thisRowLengths= rhsRowLengths; + auto matrixValuesBuffer_view = matrixValuesBuffer.getView(); + auto matrixColumnsBuffer_view = matrixColumnsBuffer.getView(); + auto thisValuesBuffer_view = thisValuesBuffer.getView(); + auto thisColumnsBuffer_view = thisColumnsBuffer.getView(); + matrixValuesBuffer_view = 0.0; + + IndexType baseRow( 0 ); + const IndexType rowsCount = this->getRows(); + while( baseRow < rowsCount ) + { + const IndexType lastRow = min( baseRow + bufferRowsCount, rowsCount ); + thisColumnsBuffer = paddingIndex; + matrixColumnsBuffer_view = paddingIndex; + + //// + // Copy matrix elements into buffer + auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value ) mutable { + if( columnIndex != paddingIndex ) + { + TNL_ASSERT_LT( rowIdx - baseRow, bufferRowsCount, "" ); + TNL_ASSERT_LT( localIdx, maxRowLength, "" ); + const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx; + TNL_ASSERT_LT( bufferIdx, ( IndexType ) bufferSize, "" ); + matrixColumnsBuffer_view[ bufferIdx ] = columnIndex; + matrixValuesBuffer_view[ bufferIdx ] = value; + } + }; + matrix.forElements( baseRow, lastRow, f1 ); + + //// + // Copy the source matrix buffer to this matrix buffer + thisValuesBuffer_view = matrixValuesBuffer_view; + thisColumnsBuffer_view = matrixColumnsBuffer_view; + + //// + // Copy matrix elements from the buffer to the matrix and ignoring + // zero matrix elements + //const IndexType matrix_columns = this->getColumns(); + const auto thisRowLengths_view = thisRowLengths.getConstView(); + auto f2 = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIndex, RealType& value ) mutable { + RealType inValue( 0.0 ); + size_t bufferIdx; + IndexType bufferLocalIdx( rowLocalIndexes_view[ rowIdx ] ); + while( inValue == 0.0 && localIdx < thisRowLengths_view[ rowIdx ] ) + { + bufferIdx = ( rowIdx - baseRow ) * maxRowLength + bufferLocalIdx++; + TNL_ASSERT_LT( bufferIdx, bufferSize, "" ); + inValue = thisValuesBuffer_view[ bufferIdx ]; + } + rowLocalIndexes_view[ rowIdx ] = bufferLocalIdx; + if( inValue == 0.0 ) + { + columnIndex = paddingIndex; + value = 0.0; + } + else + { + columnIndex = thisColumnsBuffer_view[ bufferIdx ];//column - 1; + value = inValue; + } + }; + this->forElements( baseRow, lastRow, f2 ); + baseRow += bufferRowsCount; + } + } + this->view = this->getView(); + return *this; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Matrix > +bool +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +operator==( const Matrix& m ) const +{ + return view == m; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > + template< typename Matrix > +bool +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +operator!=( const Matrix& m ) const +{ + return view != m; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +save( File& file ) const +{ + this->view.save( file ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +load( File& file ) +{ + Matrix< RealType, DeviceType, IndexType >::load( file ); + file >> this->columnIndexes; + // SANDBOX_TODO: Replace the following line with loading of metadata required by your sparse matrix format. + file >> rowPointers; + this->view = this->getView(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +save( const String& fileName ) const +{ + Object::save( fileName ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +load( const String& fileName ) +{ + Object::load( fileName ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +void +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +print( std::ostream& str ) const +{ + this->view.print( str ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +__cuda_callable__ +Index +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getPaddingIndex() const +{ + return -1; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +auto +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getColumnIndexes() const -> const ColumnsIndexesVectorType& +{ + return this->columnIndexes; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType, + typename RealAllocator, + typename IndexAllocator > +auto +SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >:: +getColumnIndexes() -> ColumnsIndexesVectorType& +{ + return this->columnIndexes; +} + + } // namespace Sandbox + } // namespace Matrices +} // namespace TNL diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h new file mode 100644 index 0000000000000000000000000000000000000000..cabf7b7fd524c19d6c1564bf2c3a5a1c4795d708 --- /dev/null +++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h @@ -0,0 +1,282 @@ + /*************************************************************************** + SparseSandboxMatrixRowView.h - description + ------------------- + begin : Apr 20, 2021 + copyright : (C) 2021 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +#include <ostream> + +#include <TNL/Cuda/CudaCallable.h> +#include <TNL/Matrices/MatrixRowViewIterator.h> + +namespace TNL { + namespace Matrices { + namespace Sandbox { + +/** + * \brief RowView is a simple structure for accessing rows of sparse matrix. + * + * \tparam ValuesView is a vector view storing the matrix elements values. + * \tparam ColumnsIndexesView is a vector view storing the column indexes of the matrix element. + * \tparam isBinary tells if the the parent matrix is a binary matrix. + * + * See \ref SparseSandboxMatrix and \ref SparseSandboxMatrixView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_getRow.cpp + * \par Output + * \include SparseMatrixExample_getRow.out + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_getRow.cpp + * \par Output + * \include SparseMatrixViewExample_getRow.out + */ +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +class SparseSandboxMatrixRowView +{ + public: + + /** + * \brief The type of matrix elements. + */ + using RealType = typename ValuesView::RealType; + + /** + * \brief The type used for matrix elements indexing. + */ + using IndexType = typename ColumnsIndexesView::IndexType; + + /** + * \brief Type of container view used for storing the matrix elements values. + */ + using ValuesViewType = ValuesView; + + /** + * \brief Type of container view used for storing the column indexes of the matrix elements. + */ + using ColumnsIndexesViewType = ColumnsIndexesView; + + /** + * \brief Type of constant container view used for storing the matrix elements values. + */ + using ConstValuesViewType = typename ValuesViewType::ConstViewType; + + /** + * \brief Type of constant container view used for storing the column indexes of the matrix elements. + */ + using ConstColumnsIndexesViewType = typename ColumnsIndexesViewType::ConstViewType; + + /** + * \brief Type of sparse matrix row view. + */ + using RowView = SparseSandboxMatrixRowView< ValuesViewType, ColumnsIndexesViewType, isBinary_ >; + + /** + * \brief Type of constant sparse matrix row view. + */ + using ConstView = SparseSandboxMatrixRowView< ConstValuesViewType, ConstColumnsIndexesViewType, isBinary_ >; + + /** + * \brief The type of related matrix element. + */ + using MatrixElementType = SparseMatrixElement< RealType, IndexType >; + + /** + * \brief Type of iterator for the matrix row. + */ + using IteratorType = MatrixRowViewIterator< RowView >; + + /** + * \brief Tells whether the parent matrix is a binary matrix. + * @return `true` if the matrix is binary. + */ + static constexpr bool isBinary() { return isBinary_; }; + + /** + * \brief Constructor with \e segmentView, \e values and \e columnIndexes. + * + * \param rowIdx is row index. + * \param offset is the begining of the matrix row in arrays with values and column indexes of matrix elements. + * \param size is row size, i.e. number of nonzero matrix elements in the row. + * \param values is a container view for storing the matrix elements values. + * \param columnIndexes is a container view for storing the column indexes of the matrix elements. + */ + __cuda_callable__ + SparseSandboxMatrixRowView( IndexType rowIdx, + IndexType offset, + IndexType size, + const ValuesViewType& values, + const ColumnsIndexesViewType& columnIndexes ); + + /** + * \brief Returns size of the matrix row, i.e. number of matrix elements in this row. + * + * \return Size of the matrix row. + */ + __cuda_callable__ + IndexType getSize() const; + + /** + * \brief Returns the matrix row index. + * + * \return matrix row index. + */ + __cuda_callable__ + const IndexType& getRowIndex() const; + + /** + * \brief Returns constants reference to a column index of an element with given rank in the row. + * + * \param localIdx is the rank of the non-zero element in given row. + * + * \return constant reference to the matrix element column index. + */ + __cuda_callable__ + const IndexType& getColumnIndex( const IndexType localIdx ) const; + + /** + * \brief Returns non-constants reference to a column index of an element with given rank in the row. + * + * \param localIdx is the rank of the non-zero element in given row. + * + * \return non-constant reference to the matrix element column index. + */ + __cuda_callable__ + IndexType& getColumnIndex( const IndexType localIdx ); + + /** + * \brief Returns constants reference to value of an element with given rank in the row. + * + * \param localIdx is the rank of the non-zero element in given row. + * + * \return constant reference to the matrix element value. + */ + __cuda_callable__ + const RealType& getValue( const IndexType localIdx ) const; + + /** + * \brief Returns non-constants reference to value of an element with given rank in the row. + * + * \param localIdx is the rank of the non-zero element in given row. + * + * \return non-constant reference to the matrix element value. + */ + __cuda_callable__ + RealType& getValue( const IndexType localIdx ); + + /** + * \brief Sets a value of matrix element with given rank in the matrix row. + * + * \param localIdx is the rank of the matrix element in the row. + * \param value is the new value of the matrix element. + */ + __cuda_callable__ + void setValue( const IndexType localIdx, + const RealType& value ); + + /** + * \brief Sets a column index of matrix element with given rank in the matrix row. + * + * \param localIdx is the rank of the matrix element in the row. + * \param columnIndex is the new column index of the matrix element. + */ + __cuda_callable__ + void setColumnIndex( const IndexType localIdx, + const IndexType& columnIndex ); + + /** + * \brief Sets both a value and a column index of matrix element with given rank in the matrix row. + * + * \param localIdx is the rank of the matrix element in the row. + * \param columnIndex is the new column index of the matrix element. + * \param value is the new value of the matrix element. + */ + __cuda_callable__ + void setElement( const IndexType localIdx, + const IndexType columnIndex, + const RealType& value ); + + /** + * \brief Comparison of two matrix rows. + * + * The other matrix row can be from any other matrix. + * + * \param other is another matrix row. + * \return \e true if both rows are the same, \e false otherwise. + */ + template< typename _ValuesView, + typename _ColumnsIndexesView, + bool _isBinary > + __cuda_callable__ + bool operator==( const SparseSandboxMatrixRowView< _ValuesView, _ColumnsIndexesView, _isBinary >& other ) const; + + /** + * \brief Returns iterator pointing at the beginning of the matrix row. + * + * \return iterator pointing at the beginning. + */ + __cuda_callable__ + IteratorType begin(); + + /** + * \brief Returns iterator pointing at the end of the matrix row. + * + * \return iterator pointing at the end. + */ + __cuda_callable__ + IteratorType end(); + + /** + * \brief Returns constant iterator pointing at the beginning of the matrix row. + * + * \return iterator pointing at the beginning. + */ + __cuda_callable__ + const IteratorType cbegin() const; + + /** + * \brief Returns constant iterator pointing at the end of the matrix row. + * + * \return iterator pointing at the end. + */ + __cuda_callable__ + const IteratorType cend() const; + + protected: + + IndexType rowIdx, size; + + // SANDBOX_TODO: Replace the following line with data required by your format. + IndexType offset; + + ValuesViewType values; + + ColumnsIndexesViewType columnIndexes; +}; + +/** + * \brief Insertion operator for a sparse matrix row. + * + * \param str is an output stream. + * \param row is an input sparse matrix row. + * \return reference to the output stream. + */ +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +std::ostream& operator<<( std::ostream& str, const SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >& row ); + + } // namespace Sandbox + } // namespace Matrices +} // namespace TNL + +#include <TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp> diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp new file mode 100644 index 0000000000000000000000000000000000000000..09598edd04390d9098902f2b6b0949d7bccd7478 --- /dev/null +++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp @@ -0,0 +1,240 @@ +/*************************************************************************** + SparseSandboxMatrixRowView.hpp - description + ------------------- + begin : Apr 20, 2021 + copyright : (C) 2021 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +#include <TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h> +#include <TNL/Assert.h> + +namespace TNL { + namespace Matrices { + namespace Sandbox { + +// SANDBOX_TODO: Modify the follwing constructor by your needs +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +SparseSandboxMatrixRowView( IndexType rowIdx, + IndexType offset, + IndexType size, + const ValuesViewType& values, + const ColumnsIndexesViewType& columnIndexes ) + : rowIdx( rowIdx ), offset( offset ), size( size ), values( values ), columnIndexes( columnIndexes ) +{ +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +getSize() const -> IndexType +{ + return this->size; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ +auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +getRowIndex() const -> const IndexType& +{ + return this->rowIdx; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +getColumnIndex( const IndexType localIdx ) const -> const IndexType& +{ + TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." ); + // SANDBOX_TODO: Modify the following line to match with your sparse format. + return columnIndexes[ offset + localIdx ]; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +getColumnIndex( const IndexType localIdx ) -> IndexType& +{ + TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." ); + // SANDBOX_TODO: Modify the following line to match with your sparse format. + return columnIndexes[ offset + localIdx ]; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +getValue( const IndexType localIdx ) const -> const RealType& +{ + TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." ); + TNL_ASSERT_FALSE( isBinary(), "Cannot call this method for binary matrix row." ); + // SANDBOX_TODO: Modify the following line to match with your sparse format. + return values[ offset + localIdx ]; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +getValue( const IndexType localIdx ) -> RealType& +{ + TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." ); + TNL_ASSERT_FALSE( isBinary(), "Cannot call this method for binary matrix row." ); + // SANDBOX_TODO: Modify the following line to match with your sparse format. + return values[ offset + localIdx ]; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ void +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +setValue( const IndexType localIdx, + const RealType& value ) +{ + TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." ); + if( ! isBinary() ) { + // SANDBOX_TODO: Modify the following line to match with your sparse format. + const IndexType globalIdx = offset + localIdx; + values[ globalIdx ] = value; + } +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ void +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +setColumnIndex( const IndexType localIdx, + const IndexType& columnIndex ) +{ + TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." ); + // SANDBOX_TODO: Modify the following line to match with your sparse format. + const IndexType globalIdx = offset + localIdx; + this->columnIndexes[ globalIdx ] = columnIndex; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ void +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +setElement( const IndexType localIdx, + const IndexType column, + const RealType& value ) +{ + TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." ); + // SANDBOX_TODO: Modify the following line to match with your sparse format. + const IndexType globalIdx = offset + localIdx; + columnIndexes[ globalIdx ] = column; + if( ! isBinary() ) + values[ globalIdx ] = value; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > + template< typename _ValuesView, + typename _ColumnsIndexesView, + bool _isBinary > +__cuda_callable__ +bool +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +operator==( const SparseSandboxMatrixRowView< _ValuesView, _ColumnsIndexesView, _isBinary >& other ) const +{ + IndexType i = 0; + while( i < getSize() && i < other.getSize() ) { + if( getColumnIndex( i ) != other.getColumnIndex( i ) ) + return false; + if( ! _isBinary && getValue( i ) != other.getValue( i ) ) + return false; + ++i; + } + for( IndexType j = i; j < getSize(); j++ ) + // TODO: use ... != getPaddingIndex() + if( getColumnIndex( j ) >= 0 ) + return false; + for( IndexType j = i; j < other.getSize(); j++ ) + // TODO: use ... != getPaddingIndex() + if( other.getColumnIndex( j ) >= 0 ) + return false; + return true; +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +begin() -> IteratorType +{ + return IteratorType( *this, 0 ); +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +end() -> IteratorType +{ + return IteratorType( *this, this->getSize() ); +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +cbegin() const -> const IteratorType +{ + return IteratorType( *this, 0 ); +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +__cuda_callable__ auto +SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >:: +cend() const -> const IteratorType +{ + return IteratorType( *this, this->getSize() ); +} + +template< typename ValuesView, + typename ColumnsIndexesView, + bool isBinary_ > +std::ostream& operator<<( std::ostream& str, const SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >& row ) +{ + using NonConstIndex = std::remove_const_t< typename SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::IndexType >; + for( NonConstIndex i = 0; i < row.getSize(); i++ ) + if( isBinary_ ) + // TODO: check getPaddingIndex(), print only the column indices of non-zeros but not the values + str << " [ " << row.getColumnIndex( i ) << " ] = " << (row.getColumnIndex( i ) >= 0) << ", "; + else + str << " [ " << row.getColumnIndex( i ) << " ] = " << row.getValue( i ) << ", "; + return str; +} + + } // namespace Sandbox + } // namespace Matrices +} // namespace TNL diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.h b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.h new file mode 100644 index 0000000000000000000000000000000000000000..66247a349dfe9b9d63e5084abe0c4b41b6c922bb --- /dev/null +++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.h @@ -0,0 +1,871 @@ +/*************************************************************************** + SparseSandboxMatrixView.h - description + ------------------- + begin : Apr 20, 2021 + copyright : (C) 2021 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +#include <TNL/Matrices/Matrix.h> +#include <TNL/Matrices/MatrixType.h> +#include <TNL/Allocators/Default.h> +#include <TNL/Algorithms/Segments/CSR.h> +#include <TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h> +#include <TNL/TypeTraits.h> + +namespace TNL { + namespace Matrices { + namespace Sandbox { + +/** + * \brief Implementation of sparse sandbox matrix view. + * + * It serves as an accessor to \ref SparseSandboxMatrix for example when passing the + * matrix to lambda functions. SparseSandboxMatrix view can be also created in CUDA kernels. + * + * \tparam Real is a type of matrix elements. If \e Real equals \e bool the matrix is treated + * as binary and so the matrix elements values are not stored in the memory since we need + * to remember only coordinates of non-zero elements( which equal one). + * \tparam Device is a device where the matrix is allocated. + * \tparam Index is a type for indexing of the matrix elements. + * \tparam MatrixType specifies a symmetry of matrix. See \ref MatrixType. Symmetric + * matrices store only lower part of the matrix and its diagonal. The upper part is reconstructed on the fly. + * GeneralMatrix with no symmetry is used by default. + * \tparam Segments is a structure representing the sparse matrix format. Depending on the pattern of the non-zero elements + * different matrix formats can perform differently especially on GPUs. By default \ref CSR format is used. See also + * \ref Ellpack, \ref SlicedEllpack, \ref ChunkedEllpack or \ref BiEllpack. + * \tparam ComputeReal is the same as \e Real mostly but for binary matrices it is set to \e Index type. This can be changed + * bu the user, of course. + * + */ +template< typename Real, + typename Device = Devices::Host, + typename Index = int, + typename MatrixType = GeneralMatrix > +class SparseSandboxMatrixView : public MatrixView< Real, Device, Index > +{ + static_assert( + ! MatrixType::isSymmetric() || + ! std::is_same< Device, Devices::Cuda >::value || + ( std::is_same< Real, float >::value || std::is_same< Real, double >::value || std::is_same< Real, int >::value || std::is_same< Real, long long int >::value ), + "Given Real type is not supported by atomic operations on GPU which are necessary for symmetric operations." ); + + public: + + // Supporting types - they are not important for the user + using BaseType = MatrixView< Real, Device, Index >; + using ValuesViewType = typename BaseType::ValuesView; + using ConstValuesViewType = typename ValuesViewType::ConstViewType; + using ColumnsIndexesViewType = Containers::VectorView< typename TNL::copy_const< Index >::template from< Real >::type, Device, Index >; + using ConstColumnsIndexesViewType = typename ColumnsIndexesViewType::ConstViewType; + using RowsCapacitiesView = Containers::VectorView< Index, Device, Index >; + using ConstRowsCapacitiesView = typename RowsCapacitiesView::ConstViewType; + + /** + * \brief Test of symmetric matrix type. + * + * \return \e true if the matrix is stored as symmetric and \e false otherwise. + */ + static constexpr bool isSymmetric() { return MatrixType::isSymmetric(); }; + + /** + * \brief Test of binary matrix type. + * + * \return \e true if the matrix is stored as binary and \e false otherwise. + */ + static constexpr bool isBinary() { return std::is_same< Real, bool >::value; }; + + /** + * \brief The type of matrix elements. + */ + using RealType = Real; + + //using ComputeRealType = ComputeReal; + + /** + * \brief The device where the matrix is allocated. + */ + using DeviceType = Device; + + /** + * \brief The type used for matrix elements indexing. + */ + using IndexType = Index; + + /** + * \brief Templated type of segments view, i.e. sparse matrix format. + */ + //template< typename Device_, typename Index_ > + //using SegmentsViewTemplate = SegmentsView< Device_, Index_ >; + + /** + * \brief Type of segments view used by this matrix. It represents the sparse matrix format. + */ + //using SegmentsViewType = SegmentsView< Device, Index >; + + /** + * \brief Type of related matrix view. + */ + using ViewType = SparseSandboxMatrixView< Real, Device, Index, MatrixType >; + + /** + * \brief Matrix view type for constant instances. + */ + using ConstViewType = SparseSandboxMatrixView< std::add_const_t< Real >, Device, Index, MatrixType >; + + /** + * \brief Type for accessing matrix rows. + */ + using RowView = SparseSandboxMatrixRowView< ValuesViewType, ColumnsIndexesViewType, isBinary() >; + + /** + * \brief Type for accessing constant matrix rows. + */ + using ConstRowView = SparseSandboxMatrixRowView< ConstValuesViewType, ConstColumnsIndexesViewType, isBinary() >;; + + /** + * \brief Helper type for getting self type or its modifications. + */ + template< typename _Real = Real, + typename _Device = Device, + typename _Index = Index, + typename _MatrixType = MatrixType > + using Self = SparseSandboxMatrixView< _Real, _Device, _Index, _MatrixType >; + + /** + * \brief Type of container view for CSR row pointers. + * + * SANDBOX_TODO: You may replace it with containers views for metadata of your format. + */ + using RowPointersView = TNL::Containers::VectorView< IndexType, DeviceType, IndexType >; + + /** + * \brief Constructor with no parameters. + */ + __cuda_callable__ + SparseSandboxMatrixView(); + + /** + * \brief Constructor with all necessary data and views. + * + * \param rows is a number of matrix rows. + * \param columns is a number of matrix columns. + * \param values is a vector view with matrix elements values. + * \param columnIndexes is a vector view with matrix elements column indexes. + * \param rowPointers is a container view with row pointers. + * + * SANDBOX_TODO: Replace `rowPointers` with metadata by your needs. + */ + __cuda_callable__ + SparseSandboxMatrixView( const IndexType rows, + const IndexType columns, + const ValuesViewType& values, + const ColumnsIndexesViewType& columnIndexes, + const RowPointersView& rowPointers ); + + /** + * \brief Copy constructor. + * + * \param matrix is an input sparse matrix view. + */ + __cuda_callable__ + SparseSandboxMatrixView( const SparseSandboxMatrixView& matrix ) = default; + + /** + * \brief Move constructor. + * + * \param matrix is an input sparse matrix view. + */ + __cuda_callable__ + SparseSandboxMatrixView( SparseSandboxMatrixView&& matrix ) = default; + + /** + * \brief Returns a modifiable view of the sparse matrix. + * + * \return sparse matrix view. + */ + __cuda_callable__ + ViewType getView(); + + /** + * \brief Returns a non-modifiable view of the sparse matrix. + * + * \return sparse matrix view. + */ + __cuda_callable__ + ConstViewType getConstView() const; + + /** + * \brief Returns string with serialization type. + * + * The string has a form `Matrices::SparseSandboxMatrix< RealType, [any_device], IndexType, General/Symmetric, Format, [any_allocator] >`. + * + * \return \ref String with the serialization type. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_getSerializationType.cpp + * \par Output + * \include SparseMatrixViewExample_getSerializationType.out + */ + static String getSerializationType(); + + /** + * \brief Returns string with serialization type. + * + * See \ref SparseSandboxMatrix::getSerializationType. + * + * \return \e String with the serialization type. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixExample_getSerializationType.cpp + * \par Output + * \include SparseMatrixExample_getSerializationType.out + */ + virtual String getSerializationTypeVirtual() const; + + /** + * \brief Computes number of non-zeros in each row. + * + * \param rowLengths is a vector into which the number of non-zeros in each row + * will be stored. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_getCompressedRowLengths.cpp + * \par Output + * \include SparseMatrixViewExample_getCompressedRowLengths.out + */ + template< typename Vector > + void getCompressedRowLengths( Vector& rowLengths ) const; + + /** + * \brief Compute capacities of all rows. + * + * The row capacities are not stored explicitly and must be computed. + * + * \param rowCapacities is a vector where the row capacities will be stored. + */ + template< typename Vector > + void getRowCapacities( Vector& rowCapacities ) const; + + /** + * \brief Returns capacity of given matrix row. + * + * \param row index of matrix row. + * \return number of matrix elements allocated for the row. + */ + __cuda_callable__ + IndexType getRowCapacity( const IndexType row ) const; + + /** + * \brief Returns number of non-zero matrix elements. + * + * This method really counts the non-zero matrix elements and so + * it returns zero for matrix having all allocated elements set to zero. + * + * \return number of non-zero matrix elements. + */ + IndexType getNonzeroElementsCount() const; + + /** + * \brief Constant getter of simple structure for accessing given matrix row. + * + * \param rowIdx is matrix row index. + * + * \return RowView for accessing given matrix row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_getConstRow.cpp + * \par Output + * \include SparseMatrixViewExample_getConstRow.out + * + * See \ref SparseSandboxMatrixRowView. + */ + __cuda_callable__ + ConstRowView getRow( const IndexType& rowIdx ) const; + + /** + * \brief Non-constant getter of simple structure for accessing given matrix row. + * + * \param rowIdx is matrix row index. + * + * \return RowView for accessing given matrix row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_getRow.cpp + * \par Output + * \include SparseMatrixViewExample_getRow.out + * + * See \ref SparseSandboxMatrixRowView. + */ + __cuda_callable__ + RowView getRow( const IndexType& rowIdx ); + + /** + * \brief Sets element at given \e row and \e column to given \e value. + * + * This method can be called from the host system (CPU) no matter + * where the matrix is allocated. If the matrix is allocated on GPU this method + * can be called even from device kernels. If the matrix is allocated in GPU device + * this method is called from CPU, it transfers values of each matrix element separately and so the + * performance is very low. For higher performance see. \ref SparseSandboxMatrix::getRow + * or \ref SparseSandboxMatrix::forElements and \ref SparseSandboxMatrix::forAllElements. + * The call may fail if the matrix row capacity is exhausted. + * + * \param row is row index of the element. + * \param column is columns index of the element. + * \param value is the value the element will be set to. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_setElement.cpp + * \par Output + * \include SparseMatrixViewExample_setElement.out + */ + __cuda_callable__ + void setElement( const IndexType row, + const IndexType column, + const RealType& value ); + + /** + * \brief Add element at given \e row and \e column to given \e value. + * + * This method can be called from the host system (CPU) no matter + * where the matrix is allocated. If the matrix is allocated on GPU this method + * can be called even from device kernels. If the matrix is allocated in GPU device + * this method is called from CPU, it transfers values of each matrix element separately and so the + * performance is very low. For higher performance see. \ref SparseSandboxMatrix::getRow + * or \ref SparseSandboxMatrix::forElements and \ref SparseSandboxMatrix::forAllElements. + * The call may fail if the matrix row capacity is exhausted. + * + * \param row is row index of the element. + * \param column is columns index of the element. + * \param value is the value the element will be set to. + * \param thisElementMultiplicator is multiplicator the original matrix element + * value is multiplied by before addition of given \e value. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_addElement.cpp + * \par Output + * \include SparseMatrixViewExample_addElement.out + */ + __cuda_callable__ + void addElement( IndexType row, + IndexType column, + const RealType& value, + const RealType& thisElementMultiplicator = 1.0 ); + + /** + * \brief Returns value of matrix element at position given by its row and column index. + * + * This method can be called from the host system (CPU) no matter + * where the matrix is allocated. If the matrix is allocated on GPU this method + * can be called even from device kernels. If the matrix is allocated in GPU device + * this method is called from CPU, it transfers values of each matrix element separately and so the + * performance is very low. For higher performance see. \ref SparseSandboxMatrix::getRow + * or \ref SparseSandboxMatrix::forElements and \ref SparseSandboxMatrix::forAllElements. + * + * \param row is a row index of the matrix element. + * \param column i a column index of the matrix element. + * + * \return value of given matrix element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_getElement.cpp + * \par Output + * \include SparseMatrixViewExample_getElement.out + * + */ + __cuda_callable__ + RealType getElement( IndexType row, + IndexType column ) const; + + /** + * \brief Method for performing general reduction on matrix rows. + * + * \tparam Fetch is a type of lambda function for data fetch declared as + * `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`. + * The return type of this lambda can be any non void. + * \tparam Reduce is a type of lambda function for reduction declared as + * `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`. + * \tparam Keep is a type of lambda function for storing results of reduction in each row. + * It is declared as `keep( const IndexType rowIdx, const double& value )`. + * \tparam FetchValue is type returned by the Fetch lambda function. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param fetch is an instance of lambda function for data fetch. + * \param reduce is an instance of lambda function for reduction. + * \param keep in an instance of lambda function for storing results. + * \param zero is zero of given reduction operation also known as idempotent element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_reduceRows.cpp + * \par Output + * \include SparseMatrixViewExample_reduceRows.out + */ + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > + void reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ); + + /** + * \brief Method for performing general reduction on matrix rows for constant instances. + * + * \tparam Fetch is a type of lambda function for data fetch declared as + * `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`. + * The return type of this lambda can be any non void. + * \tparam Reduce is a type of lambda function for reduction declared as + * `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`. + * \tparam Keep is a type of lambda function for storing results of reduction in each row. + * It is declared as `keep( const IndexType rowIdx, const double& value )`. + * \tparam FetchValue is type returned by the Fetch lambda function. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param fetch is an instance of lambda function for data fetch. + * \param reduce is an instance of lambda function for reduction. + * \param keep in an instance of lambda function for storing results. + * \param zero is zero of given reduction operation also known as idempotent element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_reduceRows.cpp + * \par Output + * \include SparseMatrixViewExample_reduceRows.out + */ + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > + void reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const; + + /** + * \brief Method for performing general reduction on all matrix rows. + * + * \tparam Fetch is a type of lambda function for data fetch declared as + * `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`. + * The return type of this lambda can be any non void. + * \tparam Reduce is a type of lambda function for reduction declared as + * `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`. + * \tparam Keep is a type of lambda function for storing results of reduction in each row. + * It is declared as `keep( const IndexType rowIdx, const double& value )`. + * \tparam FetchValue is type returned by the Fetch lambda function. + * + * \param fetch is an instance of lambda function for data fetch. + * \param reduce is an instance of lambda function for reduction. + * \param keep in an instance of lambda function for storing results. + * \param zero is zero of given reduction operation also known as idempotent element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_reduceAllRows.cpp + * \par Output + * \include SparseMatrixViewExample_reduceAllRows.out + */ + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > + void reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ); + + /** + * \brief Method for performing general reduction on all matrix rows for constant instances. + * + * \tparam Fetch is a type of lambda function for data fetch declared as + * `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`. + * The return type of this lambda can be any non void. + * \tparam Reduce is a type of lambda function for reduction declared as + * `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`. + * \tparam Keep is a type of lambda function for storing results of reduction in each row. + * It is declared as `keep( const IndexType rowIdx, const double& value )`. + * \tparam FetchValue is type returned by the Fetch lambda function. + * + * \param fetch is an instance of lambda function for data fetch. + * \param reduce is an instance of lambda function for reduction. + * \param keep in an instance of lambda function for storing results. + * \param zero is zero of given reduction operation also known as idempotent element. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_reduceAllRows.cpp + * \par Output + * \include SparseMatrixViewExample_reduceAllRows.out + */ + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > + void reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const; + + /** + * \brief Method for iteration over all matrix rows for constant instances. + * + * \tparam Function is type of lambda function that will operate on matrix elements. + * It is should have form like + * `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`. + * The \e localIdx parameter is a rank of the non-zero element in given row. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param function is an instance of the lambda function to be called in each row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp + * \par Output + * \include SparseMatrixViewExample_forRows.out + */ + template< typename Function > + void forElements( IndexType begin, IndexType end, Function& function ) const; + + /** + * \brief Method for iteration over all matrix rows for non-constant instances. + * + * \tparam Function is type of lambda function that will operate on matrix elements. + * It is should have form like + * `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`. + * The \e localIdx parameter is a rank of the non-zero element in given row. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param function is an instance of the lambda function to be called in each row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp + * \par Output + * \include SparseMatrixViewExample_forRows.out + */ + template< typename Function > + void forElements( IndexType begin, IndexType end, Function& function ); + + /** + * \brief This method calls \e forElements for all matrix rows (for constant instances). + * + * See \ref SparseSandboxMatrix::forElements. + * + * \tparam Function is a type of lambda function that will operate on matrix elements. + * \param function is an instance of the lambda function to be called in each row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_forAllRows.cpp + * \par Output + * \include SparseMatrixViewExample_forAllRows.out + */ + template< typename Function > + void forAllElements( Function& function ) const; + + /** + * \brief This method calls \e forElements for all matrix rows. + * + * See \ref SparseSandboxMatrix::forElements. + * + * \tparam Function is a type of lambda function that will operate on matrix elements. + * \param function is an instance of the lambda function to be called in each row. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_forAllRows.cpp + * \par Output + * \include SparseMatrixViewExample_forAllRows.out + */ + template< typename Function > + void forAllElements( Function& function ); + + /** + * \brief Method for parallel iteration over matrix rows from interval [ \e begin, \e end). + * + * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method + * \ref SparseSandboxMatrixView::forElements where more than one thread can be mapped to each row. + + * + * \tparam Function is type of the lambda function. + * + * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed. + * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed. + * \param function is an instance of the lambda function to be called for each row. + * + * ``` + * auto function = [] __cuda_callable__ ( RowView& row ) mutable { ... }; + * ``` + * + * \e RowView represents matrix row - see \ref TNL::Matrices::SparseSandboxMatrixView::RowView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp + * \par Output + * \include SparseMatrixViewExample_forRows.out + */ + template< typename Function > + void forRows( IndexType begin, IndexType end, Function&& function ); + + /** + * \brief Method for parallel iteration over matrix rows from interval [ \e begin, \e end) for constant instances. + * + * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method + * \ref SparseSandboxMatrixView::forElements where more than one thread can be mapped to each row. + * + * \tparam Function is type of the lambda function. + * + * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed. + * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed. + * \param function is an instance of the lambda function to be called for each row. + * + * ``` + * auto function = [] __cuda_callable__ ( RowView& row ) { ... }; + * ``` + * + * \e RowView represents matrix row - see \ref TNL::Matrices::SparseSandboxMatrixView::RowView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp + * \par Output + * \include SparseMatrixViewExample_forRows.out + */ + template< typename Function > + void forRows( IndexType begin, IndexType end, Function&& function ) const; + + /** + * \brief Method for parallel iteration over all matrix rows. + * + * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method + * \ref SparseSandboxMatrixView::forAllElements where more than one thread can be mapped to each row. + * + * \tparam Function is type of the lambda function. + * + * \param function is an instance of the lambda function to be called for each row. + * + * ``` + * auto function = [] __cuda_callable__ ( RowView& row ) mutable { ... }; + * ``` + * + * \e RowView represents matrix row - see \ref TNL::Matrices::SparseSandboxMatrixView::RowView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp + * \par Output + * \include SparseMatrixViewExample_forRows.out + */ + template< typename Function > + void forAllRows( Function&& function ); + + /** + * \brief Method for parallel iteration over all matrix rows for constant instances. + * + * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method + * \ref SparseSandboxMatrixView::forAllElements where more than one thread can be mapped to each row. + * + * \tparam Function is type of the lambda function. + * + * \param function is an instance of the lambda function to be called for each row. + * + * ``` + * auto function = [] __cuda_callable__ ( RowView& row ) { ... }; + * ``` + * + * \e RowView represents matrix row - see \ref TNL::Matrices::SparseSandboxMatrixView::RowView. + * + * \par Example + * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp + * \par Output + * \include SparseMatrixViewExample_forRows.out + */ + template< typename Function > + void forAllRows( Function&& function ) const; + + /** + * \brief Method for sequential iteration over all matrix rows for constant instances. + * + * \tparam Function is type of lambda function that will operate on matrix elements. + * It is should have form like + * `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`. + * The column index repeats twice only for compatibility with sparse matrices. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param function is an instance of the lambda function to be called in each row. + */ + template< typename Function > + void sequentialForRows( IndexType begin, IndexType end, Function& function ) const; + + /** + * \brief Method for sequential iteration over all matrix rows for non-constant instances. + * + * \tparam Function is type of lambda function that will operate on matrix elements. + * It is should have form like + * `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`. + * The column index repeats twice only for compatibility with sparse matrices. + * + * \param begin defines beginning of the range [begin,end) of rows to be processed. + * \param end defines ending of the range [begin,end) of rows to be processed. + * \param function is an instance of the lambda function to be called in each row. + */ + template< typename Function > + void sequentialForRows( IndexType begin, IndexType end, Function& function ); + + /** + * \brief This method calls \e sequentialForRows for all matrix rows (for constant instances). + * + * See \ref SparseSandboxMatrixView::sequentialForRows. + * + * \tparam Function is a type of lambda function that will operate on matrix elements. + * \param function is an instance of the lambda function to be called in each row. + */ + template< typename Function > + void sequentialForAllRows( Function& function ) const; + + /** + * \brief This method calls \e sequentialForRows for all matrix rows. + * + * See \ref SparseSandboxMatrixView::sequentialForAllRows. + * + * \tparam Function is a type of lambda function that will operate on matrix elements. + * \param function is an instance of the lambda function to be called in each row. + */ + template< typename Function > + void sequentialForAllRows( Function& function ); + + /** + * \brief Computes product of matrix and vector. + * + * More precisely, it computes: + * + * `outVector = matrixMultiplicator * ( * this ) * inVector + outVectorMultiplicator * outVector` + * + * \tparam InVector is type of input vector. It can be \ref Vector, + * \ref VectorView, \ref Array, \ref ArraView or similar container. + * \tparam OutVector is type of output vector. It can be \ref Vector, + * \ref VectorView, \ref Array, \ref ArraView or similar container. + * + * \param inVector is input vector. + * \param outVector is output vector. + * \param matrixMultiplicator is a factor by which the matrix is multiplied. It is one by default. + * \param outVectorMultiplicator is a factor by which the outVector is multiplied before added + * to the result of matrix-vector product. It is zero by default. + * \param begin is the beginning of the rows range for which the vector product + * is computed. It is zero by default. + * \param end is the end of the rows range for which the vector product + * is computed. It is number if the matrix rows by default. + */ + template< typename InVector, + typename OutVector > + void vectorProduct( const InVector& inVector, + OutVector& outVector, + const RealType matrixMultiplicator = 1.0, + const RealType outVectorMultiplicator = 0.0, + const IndexType begin = 0, + IndexType end = 0 ) const; + + template< typename Vector1, typename Vector2 > + bool performSORIteration( const Vector1& b, + const IndexType row, + Vector2& x, + const RealType& omega = 1.0 ) const; + + /** + * \brief Assignment of any matrix type. + * . + * \param matrix is input matrix for the assignment. + * \return reference to this matrix. + */ + SparseSandboxMatrixView& operator=( const SparseSandboxMatrixView& matrix ); + + /** + * \brief Comparison operator with another arbitrary matrix type. + * + * \param matrix is the right-hand side matrix. + * \return \e true if the RHS matrix is equal, \e false otherwise. + */ + template< typename Matrix > + bool operator==( const Matrix& m ) const; + + /** + * \brief Comparison operator with another arbitrary matrix type. + * + * \param matrix is the right-hand side matrix. + * \return \e true if the RHS matrix is equal, \e false otherwise. + */ + template< typename Matrix > + bool operator!=( const Matrix& m ) const; + + /** + * \brief Method for saving the matrix to the file with given filename. + * + * \param fileName is name of the file. + */ + void save( const String& fileName ) const; + + /** + * \brief Method for saving the matrix to a file. + * + * \param file is the output file. + */ + void save( File& file ) const; + + /** + * \brief Method for printing the matrix to output stream. + * + * \param str is the output stream. + */ + void print( std::ostream& str ) const; + + /** + * \brief Getter of segments for non-constant instances. + * + * \e Segments are a structure for addressing the matrix elements columns and values. + * In fact, \e Segments represent the sparse matrix format. + * + * \return Non-constant reference to segments. + */ + //SegmentsViewType& getSegments(); + + /** + * \brief Getter of segments for constant instances. + * + * \e Segments are a structure for addressing the matrix elements columns and values. + * In fact, \e Segments represent the sparse matrix format. + * + * \return Constant reference to segments. + */ + //const SegmentsViewType& getSegments() const; + + /** + * \brief Getter of column indexes for constant instances. + * + * \return Constant reference to a vector with matrix elements column indexes. + */ + const ColumnsIndexesViewType& getColumnIndexes() const; + + /** + * \brief Getter of column indexes for nonconstant instances. + * + * \return Reference to a vector with matrix elements column indexes. + */ + ColumnsIndexesViewType& getColumnIndexes(); + + /** + * \brief Returns a padding index value. + * + * Padding index is used for column indexes of padding zeros. Padding zeros + * are used in some sparse matrix formats for better data alignment in memory. + * + * \return value of the padding index. + */ + __cuda_callable__ + IndexType getPaddingIndex() const; + + protected: + + ColumnsIndexesViewType columnIndexes; + + RowPointersView rowPointers; + //SegmentsViewType segments; + + private: + // TODO: this should be probably moved into a detail namespace + template< typename VectorOrView, + std::enable_if_t< HasSetSizeMethod< VectorOrView >::value, bool > = true > + static void set_size_if_resizable( VectorOrView& v, IndexType size ) + { + v.setSize( size ); + } + + template< typename VectorOrView, + std::enable_if_t< ! HasSetSizeMethod< VectorOrView >::value, bool > = true > + static void set_size_if_resizable( VectorOrView& v, IndexType size ) + { + TNL_ASSERT_EQ( v.getSize(), size, "view has wrong size" ); + } +}; + + } // namespace Sandbox + } // namespace Matrices +} // namespace TNL + +#include <TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp> diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp new file mode 100644 index 0000000000000000000000000000000000000000..fb50a7acf1db76bc40c0de0ade441206d51ba1fb --- /dev/null +++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp @@ -0,0 +1,1026 @@ +/*************************************************************************** + SparseSandboxMatrixView.hpp - description + ------------------- + begin : Apr 20, 2021 + copyright : (C) 2021 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +#include <functional> +#include <TNL/Matrices/Sandbox/SparseSandboxMatrixView.h> +#include <TNL/Algorithms/Reduction.h> +#include <TNL/Algorithms/AtomicOperations.h> +#include <TNL/Matrices/details/SparseMatrix.h> + +namespace TNL { + namespace Matrices { + namespace Sandbox { + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +SparseSandboxMatrixView() +{ +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +SparseSandboxMatrixView( const IndexType rows, + const IndexType columns, + const ValuesViewType& values, + const ColumnsIndexesViewType& columnIndexes, + const RowPointersView& rowPointers ) +: MatrixView< Real, Device, Index >( rows, columns, values ), columnIndexes( columnIndexes ), rowPointers( rowPointers ) +{ +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ +auto +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getView() -> ViewType +{ + return ViewType( this->getRows(), + this->getColumns(), + this->getValues().getView(), + this->columnIndexes.getView(), + this->segments.getView() ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ +auto +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getConstView() const -> ConstViewType +{ + return ConstViewType( this->getRows(), + this->getColumns(), + this->getValues().getConstView(), + this->getColumnsIndexes().getConstView(), + this->segments.getConstView() ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +String +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getSerializationType() +{ + return String( "Matrices::Sandbox::SparseMatrix< " ) + + TNL::getSerializationType< RealType >() + ", " + + TNL::getSerializationType< IndexType >() + ", " + + MatrixType::getSerializationType() + ", [any_allocator], [any_allocator] >"; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +String +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getSerializationTypeVirtual() const +{ + return this->getSerializationType(); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Vector > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getCompressedRowLengths( Vector& rowLengths ) const +{ + details::set_size_if_resizable( rowLengths, this->getRows() ); + rowLengths = 0; + auto rowLengths_view = rowLengths.getView(); + auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType { + return ( value != 0.0 ); + }; + auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable { + rowLengths_view[ rowIdx ] = value; + }; + this->reduceAllRows( fetch, std::plus<>{}, keep, 0 ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Vector > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getRowCapacities( Vector& rowLengths ) const +{ + details::set_size_if_resizable( rowLengths, this->getRows() ); + rowLengths = 0; + auto rowLengths_view = rowLengths.getView(); + auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType { + return 1; + }; + auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable { + rowLengths_view[ rowIdx ] = value; + }; + this->reduceAllRows( fetch, std::plus<>{}, keep, 0 ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ +Index +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getRowCapacity( const IndexType row ) const +{ + return this->segments.getSegmentSize( row ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +Index +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getNonzeroElementsCount() const +{ + const auto columns_view = this->columnIndexes.getConstView(); + const IndexType paddingIndex = this->getPaddingIndex(); + if( ! isSymmetric() ) + { + auto fetch = [=] __cuda_callable__ ( const IndexType i ) -> IndexType { + return ( columns_view[ i ] != paddingIndex ); + }; + return Algorithms::Reduction< DeviceType >::reduce( ( IndexType ) 0, this->columnIndexes.getSize(), fetch, std::plus<>{}, 0 ); + } + else + { + const auto rows = this->getRows(); + const auto columns = this->getColumns(); + Containers::Vector< IndexType, DeviceType, IndexType > row_sums( this->getRows(), 0 ); + auto row_sums_view = row_sums.getView(); + auto row_pointers_view = this->rowPointers.getConstView(); + const auto columnIndexesView = this->columnIndexes.getConstView(); + // SANDBOX_TODO: Replace the following lambda function (or more) with code compute number of nonzero matrix elements + // of symmetric matrix. Note, that this is required only by symmetric matrices and that the highest performance + // is not a priority here. + auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable { + auto begin = row_pointers_view[ rowIdx ]; + auto end = row_pointers_view[ rowIdx + 1 ]; + IndexType sum( 0 ); + for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ ) + { + const IndexType column = columnIndexesView[ globalIdx ]; + if( column != paddingIndex ) + sum += 1 + ( column != rowIdx && column < rows && rowIdx < columns ); + } + row_sums_view[ rowIdx ] = sum; + }; + TNL::Algorithms::ParallelFor< DeviceType >::exec( ( IndexType ) 0, this->getRows(), f ); + return sum( row_sums ); + } + return 0; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ auto +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getRow( const IndexType& rowIdx ) const -> ConstRowView +{ + TNL_ASSERT_LT( rowIdx, this->getRows(), "Row index is larger than number of matrix rows." ); + // SANDBOX_TODO: Replace the following with creation of RowView corresponding with your sparse matrix format. + return ConstRowView( rowIdx, // row index + this->rowPointers[ rowIdx ], // row begining + this->rowPointers[ rowIdx + 1 ] - this->rowPointers[ rowIdx ], // number of elemnts allocated for given row + this->values, this->columnIndexes ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ auto +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getRow( const IndexType& rowIdx ) -> RowView +{ + TNL_ASSERT_LT( rowIdx, this->getRows(), "Row index is larger than number of matrix rows." ); + // SANDBOX_TODO: Replace this with RowView constructor by your needs. + return RowView( rowIdx, // row index + rowPointers[ rowIdx ], // index of the first nonzero element in the row + rowPointers[ rowIdx + 1 ] - rowPointers[ rowIdx ], // number of nonzero elements in the row + this->values, + this->columnIndexes ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +setElement( const IndexType row, + const IndexType column, + const RealType& value ) +{ + this->addElement( row, column, value, 0.0 ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +addElement( IndexType row, + IndexType column, + const RealType& value, + const RealType& thisElementMultiplicator ) +{ + TNL_ASSERT_GE( row, 0, "Sparse matrix row index cannot be negative." ); + TNL_ASSERT_LT( row, this->getRows(), "Sparse matrix row index is larger than number of matrix rows." ); + TNL_ASSERT_GE( column, 0, "Sparse matrix column index cannot be negative." ); + TNL_ASSERT_LT( column, this->getColumns(), "Sparse matrix column index is larger than number of matrix columns." ); + + if( isSymmetric() && row < column ) + { + swap( row, column ); + TNL_ASSERT_LT( row, this->getRows(), "Column index is out of the symmetric part of the matrix after transposition." ); + TNL_ASSERT_LT( column,this->getColumns(), "Row index is out of the symmetric part of the matrix after transposition." ); + } + + // SANDBOX_TODO: Replace the following line with a code that computes number of matrix elements allocated for + // matrix row with indedx `row`. Note that the code works on both host and GPU kernel. To achieve + // the same effect, you may use macro __CUDA_ARCH__ as can be seen bellow in this method. + const IndexType rowSize = this->rowPointers.getElement( row + 1 ) - this->rowPointers.getElement( row ); + IndexType col( this->getPaddingIndex() ); + IndexType i; + IndexType globalIdx; + for( i = 0; i < rowSize; i++ ) + { + // SANDBOX_TODO: Replace the following line with a code that computes a global index of `i`-th nonzero matrix element + // in the `row`-th matrix row. The global index is a pointer to arrays `values` and `columnIndexes` storing + // the matrix elements values and column indexes respectively. + globalIdx = this->rowPointers.getElement( row ) + i; + TNL_ASSERT_LT( globalIdx, this->columnIndexes.getSize(), "" ); + col = this->columnIndexes.getElement( globalIdx ); + if( col == column ) + { + if( ! isBinary() ) + this->values.setElement( globalIdx, thisElementMultiplicator * this->values.getElement( globalIdx ) + value ); + return; + } + if( col == this->getPaddingIndex() || col > column ) + break; + } + if( i == rowSize ) + { +#ifndef __CUDA_ARCH__ + std::stringstream msg; + msg << "The capacity of the sparse matrix row number " << row << " was exceeded."; + throw std::logic_error( msg.str() ); +#else + TNL_ASSERT_TRUE( false, ""); + return; +#endif + } + if( col == this->getPaddingIndex() ) + { + this->columnIndexes.setElement( globalIdx, column ); + if( ! isBinary() ) + this->values.setElement( globalIdx, value ); + return; + } + else + { + IndexType j = rowSize - 1; + while( j > i ) + { + // SANDBOX_TODO: Replace the following two lines with a code that computes a global indexes of `j`-th and `j-1`-th nonzero matrix elements + // in the `row`-th matrix row. The global index is a pointer to arrays `values` and `columnIndexes` storing + // the matrix elements values and column indexes respectively. + const IndexType globalIdx1 = this->rowPointers.getElement( row ) + j; + const IndexType globalIdx2 = globalIdx1 - 1; + // End of code replacement. + TNL_ASSERT_LT( globalIdx1, this->columnIndexes.getSize(), "" ); + TNL_ASSERT_LT( globalIdx2, this->columnIndexes.getSize(), "" ); + this->columnIndexes.setElement( globalIdx1, this->columnIndexes.getElement( globalIdx2 ) ); + if( ! isBinary() ) + this->values.setElement( globalIdx1, this->values.getElement( globalIdx2 ) ); + j--; + } + + this->columnIndexes.setElement( globalIdx, column ); + if( ! isBinary() ) + this->values.setElement( globalIdx, value ); + return; + } +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ +auto +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getElement( IndexType row, + IndexType column ) const -> RealType +{ + TNL_ASSERT_GE( row, 0, "Sparse matrix row index cannot be negative." ); + TNL_ASSERT_LT( row, this->getRows(), "Sparse matrix row index is larger than number of matrix rows." ); + TNL_ASSERT_GE( column, 0, "Sparse matrix column index cannot be negative." ); + TNL_ASSERT_LT( column, this->getColumns(), "Sparse matrix column index is larger than number of matrix columns." ); + + if( isSymmetric() && row < column ) + { + swap( row, column ); + if( row >= this->getRows() || column >= this->getColumns() ) + return 0.0; + } + + // SANDBOX_TODO: Replace the following lines with a code for getting number of elements allocated for given row. + const IndexType rowSize = this->rowPointers.getElement( row + 1 ) - this->rowPointers.getElement( row ); + for( IndexType i = 0; i < rowSize; i++ ) + { + // SANDBOX_TODO: Replace the following line with a code for getting index of the matrix element in arrays `values` and `columnIdexes`. + const IndexType globalIdx = this->rowPointers.getElement( row ) + i; + TNL_ASSERT_LT( globalIdx, this->columnIndexes.getSize(), "" ); + const IndexType col = this->columnIndexes.getElement( globalIdx ); + if( col == column ) + { + if( isBinary() ) + return 1; + else + return this->values.getElement( globalIdx ); + } + } + return 0.0; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +template< typename InVector, + typename OutVector > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +vectorProduct( const InVector& inVector, + OutVector& outVector, + const RealType matrixMultiplicator, + const RealType outVectorMultiplicator, + const IndexType firstRow, + IndexType lastRow ) const +{ + TNL_ASSERT_EQ( this->getColumns(), inVector.getSize(), "Matrix columns do not fit with input vector." ); + TNL_ASSERT_EQ( this->getRows(), outVector.getSize(), "Matrix rows do not fit with output vector." ); + + using OutVectorReal = typename OutVector::RealType; + static_assert( + ! MatrixType::isSymmetric() || + ! std::is_same< Device, Devices::Cuda >::value || + ( std::is_same< OutVectorReal, float >::value || + std::is_same< OutVectorReal, double >::value || + std::is_same< OutVectorReal, int >::value || + std::is_same< OutVectorReal, long long int >::value ), + "Given Real type is not supported by atomic operations on GPU which are necessary for symmetric operations." ); + + const auto inVectorView = inVector.getConstView(); + auto outVectorView = outVector.getView(); + const auto valuesView = this->values.getConstView(); + const auto columnIndexesView = this->columnIndexes.getConstView(); + const auto rowPointersView = this->rowPointers.getConstView(); + const IndexType paddingIndex = this->getPaddingIndex(); +#define HAVE_SANDBOX_SIMPLE_SPMV + // SANDBOX_TODO: The following is simple direct implementation of SpMV operation with CSR format. We recommend to start by + // replacing this part with SpMV based on your sparse format. + if( std::is_same< DeviceType, TNL::Devices::Host >::value ) // this way you may easily specialize for different device types + { + // SANDBOX_TODO: This simple and naive implementation for CPU. + for( IndexType rowIdx = firstRow; rowIdx < lastRow; rowIdx++ ) + { + const auto begin = rowPointers[ rowIdx ]; + const auto end = rowPointers[ rowIdx + 1 ]; + RealType sum( 0.0 ); + for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ ) + sum += this->values[ globalIdx ] * inVector[ this->columnIndexes[ globalIdx ] ]; + // SANDBOX_TODO:The following is quite inefficient, its better to specialized the code for cases when + // `outVectorMultiplicator` is zero or `matrixMultiplicator` is one - see. the full implementation bellow. + outVector[ rowIdx ] = outVector[ rowIdx ] * outVectorMultiplicator + matrixMultiplicator * sum; + } + } + else + { + //SANDBOX_TODO: The following is general implementation based on ParallelFor and lambda function. It would work even on CPU. + auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable { + const auto begin = rowPointersView[ rowIdx ]; + const auto end = rowPointersView[ rowIdx + 1 ]; + RealType sum( 0.0 ); + for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ ) + sum + valuesView[ globalIdx ] * inVectorView[ columnIndexesView[ globalIdx ] ]; + outVectorView[ rowIdx ] = outVectorView[ rowIdx ] * outVectorMultiplicator + matrixMultiplicator * sum; + }; + TNL::Algorithms::ParallelFor< DeviceType >::exec( firstRow, lastRow, f ); + } +#ifdef HAVE_SANDBOX_SIMPLE_SPMV +#else + // SANDBOX_TODO: The following is fully functional implementation based on method `reduceRows`. + if( isSymmetric() ) + outVector *= outVectorMultiplicator; + auto symmetricFetch = [=] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType globalIdx, bool& compute ) mutable -> RealType { + const IndexType column = columnIndexesView[ globalIdx ]; + compute = ( column != paddingIndex ); + if( ! compute ) + return 0.0; + if( isSymmetric() && column < row ) + { + if( isBinary() ) + Algorithms::AtomicOperations< DeviceType >::add( outVectorView[ column ], ( OutVectorReal ) matrixMultiplicator * inVectorView[ row ] ); + else + Algorithms::AtomicOperations< DeviceType >::add( outVectorView[ column ], ( OutVectorReal ) matrixMultiplicator * valuesView[ globalIdx ] * inVectorView[ row ] ); + } + if( isBinary() ) + return inVectorView[ column ]; + return valuesView[ globalIdx ] * inVectorView[ column ]; + }; + auto fetch = [=] __cuda_callable__ ( IndexType globalIdx, bool& compute ) mutable -> RealType { + const IndexType column = columnIndexesView[ globalIdx ]; + if( isBinary() ) + return inVectorView[ column ]; + return valuesView[ globalIdx ] * inVectorView[ column ]; + }; + + auto keeperGeneral = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable { + if( isSymmetric() ) + { + typename OutVector::RealType aux = matrixMultiplicator * value; + Algorithms::AtomicOperations< DeviceType >::add( outVectorView[ row ], aux ); + } + else + { + if( outVectorMultiplicator == 0.0 ) + outVectorView[ row ] = matrixMultiplicator * value; + else + outVectorView[ row ] = outVectorMultiplicator * outVectorView[ row ] + matrixMultiplicator * value; + } + }; + auto keeperDirect = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable { + outVectorView[ row ] = value; + }; + auto keeperMatrixMult = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable { + outVectorView[ row ] = matrixMultiplicator * value; + }; + auto keeperVectorMult = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable { + outVectorView[ row ] = outVectorMultiplicator * outVectorView[ row ] + value; + }; + + if( lastRow == 0 ) + lastRow = this->getRows(); + if( isSymmetric() ) + this->reduceRows( firstRow, lastRow, symmetricFetch, std::plus<>{}, keeperGeneral, ( RealType ) 0.0 ); + else + { + if( outVectorMultiplicator == 0.0 ) + { + if( matrixMultiplicator == 1.0 ) + this->reduceRows( firstRow, lastRow, fetch, std::plus<>{}, keeperDirect, ( RealType ) 0.0 ); + else + this->reduceRows( firstRow, lastRow, fetch, std::plus<>{}, keeperMatrixMult, ( RealType ) 0.0 ); + } + else + { + if( matrixMultiplicator == 1.0 ) + this->reduceRows( firstRow, lastRow, fetch, std::plus<>{}, keeperVectorMult, ( RealType ) 0.0 ); + else + this->reduceRows( firstRow, lastRow, fetch, std::plus<>{}, keeperGeneral, ( RealType ) 0.0 ); + } + } +#endif +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Fetch, typename Reduce, typename Keep, typename FetchValue > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchValue& zero ) +{ + auto columns_view = this->columnIndexes.getView(); + auto values_view = this->values.getView(); + auto row_pointers_view = this->rowPointers.getConstView(); + const IndexType paddingIndex_ = this->getPaddingIndex(); + // SANDBOX_TODO: Replace the following code with the one for computing reduction in rows by your format. + // Note, that this method can be used for implementation of SpMV. + auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable { + const auto begin = row_pointers_view[ rowIdx ]; + const auto end = row_pointers_view[ rowIdx + 1 ]; + FetchValue sum = zero; + for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ ) + { + IndexType& columnIdx = columns_view[ globalIdx ]; + if( columnIdx != paddingIndex_ ) + { + if( isBinary() ) + sum = reduce( sum, fetch( rowIdx, columnIdx, 1 ) ); + else + sum = reduce( sum, fetch( rowIdx, columnIdx, values_view[ globalIdx ] ) ); + } + } + keep( rowIdx, sum ); + }; + TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Fetch, typename Reduce, typename Keep, typename FetchValue > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchValue& zero ) const +{ + auto columns_view = this->columnIndexes.getConstView(); + auto values_view = this->values.getConstView(); + const IndexType paddingIndex_ = this->getPaddingIndex(); + // SANDBOX_TODO: Replace the following code with the one for computing reduction in rows by your format. + // Note, that this method can be used for implementation of SpMV. + auto row_pointers_view = this->rowPointers.getConstView(); + auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable { + const auto begin = row_pointers_view[ rowIdx ]; + const auto end = row_pointers_view[ rowIdx + 1 ]; + FetchValue sum = zero; + for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ ) + { + const IndexType& columnIdx = columns_view[ globalIdx ]; + if( columnIdx != paddingIndex_ ) + { + if( isBinary() ) + sum = reduce( sum, fetch( rowIdx, columnIdx, 1 ) ); + else + sum = reduce( sum, fetch( rowIdx, columnIdx, values_view[ globalIdx ] ) ); + } + } + keep( rowIdx, sum ); + }; + TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) +{ + this->reduceRows( 0, this->getRows(), fetch, reduce, keep, zero ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Fetch, typename Reduce, typename Keep, typename FetchReal > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const +{ + this->reduceRows( 0, this->getRows(), fetch, reduce, keep, zero ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +forElements( IndexType begin, IndexType end, Function& function ) const +{ + const auto columns_view = this->columnIndexes.getConstView(); + const auto values_view = this->values.getConstView(); + // SANDBOX_TODO: Replace the following code with the one for iterating over all allocated matrix elements. + auto row_pointers_view = this->rowPointers.getConstView(); + auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable { + const auto begin = row_pointers_view[ rowIdx ]; + const auto end = row_pointers_view[ rowIdx + 1 ]; + IndexType localIdx( 0 ); + for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ ) + { + if( isBinary() ) + function( rowIdx, localIdx, columns_view[ globalIdx ], ( RealType ) 1 ); + else + function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ] ); + localIdx++; + } + }; + TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +forElements( IndexType begin, IndexType end, Function& function ) +{ + auto columns_view = this->columnIndexes.getView(); + auto values_view = this->values.getView(); + // SANDBOX_TODO: Replace the following code with the one for iterating over all allocated matrix elements. + auto row_pointers_view = this->rowPointers.getConstView(); + auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable { + const auto begin = row_pointers_view[ rowIdx ]; + const auto end = row_pointers_view[ rowIdx + 1 ]; + IndexType localIdx( 0 ); + RealType one( 1.0 ); + for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ ) + { + if( isBinary() ) + function( rowIdx, localIdx, columns_view[ globalIdx ], one ); // TODO: Fix this without using `one`. + else + function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ] ); + localIdx++; + } + }; + TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +forAllElements( Function& function ) const +{ + this->forElements( 0, this->getRows(), function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +forAllElements( Function& function ) +{ + this->forElements( 0, this->getRows(), function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +forRows( IndexType begin, IndexType end, Function&& function ) +{ + auto columns_view = this->columnIndexes.getView(); + auto values_view = this->values.getView(); + // SANDBOX_TODO: Replace the following code with the one for iteration over matrix rows. + auto row_pointers_view = this->rowPointers.getConstView(); + auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable { + auto rowView = RowView( rowIdx, // row index + row_pointers_view[ rowIdx ], // row begining + row_pointers_view[ rowIdx + 1 ] -row_pointers_view[ rowIdx ], // number of elemnts allocated for given matrix row + values_view, columns_view ); + function( rowView ); + }; + TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +forRows( IndexType begin, IndexType end, Function&& function ) const +{ + const auto columns_view = this->columnIndexes.getConstView(); + const auto values_view = this->values.getConstView(); + // SANDBOX_TODO: Replace the following code with the one for iteration over matrix rows. + auto row_pointers_view = this->rowPointers.getConstView(); + auto f = [=] __cuda_callable__ ( IndexType rowIdx ) { + auto rowView = ConstRowView( rowIdx, // row index + row_pointers_view[ rowIdx ], // row begining + row_pointers_view[ rowIdx + 1 ] -row_pointers_view[ rowIdx ], // number of elemnts allocated for given matrix row + values_view, columns_view ); + function( rowView ); + }; + TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +forAllRows( Function&& function ) +{ + this->forRows( 0, this->getRows(), function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +forAllRows( Function&& function ) const +{ + this->forRows( 0, this->getRows(), function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +sequentialForRows( IndexType begin, IndexType end, Function& function ) const +{ + for( IndexType row = begin; row < end; row ++ ) + this->forRows( row, row + 1, function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +sequentialForRows( IndexType begin, IndexType end, Function& function ) +{ + for( IndexType row = begin; row < end; row ++ ) + this->forRows( row, row + 1, function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +sequentialForAllRows( Function& function ) const +{ + this->sequentialForRows( 0, this->getRows(), function ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Function > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +sequentialForAllRows( Function& function ) +{ + this->sequentialForRows( 0, this->getRows(), function ); +} + +/*template< typename Real, + template< typename, typename > class SegmentsView, + typename Device, + typename Index, + typename RealAllocator, + typename IndexAllocator > +template< typename Real2, template< typename, typename > class Segments2, typename Index2, typename RealAllocator2, typename IndexAllocator2 > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +addMatrix( const SparseSandboxMatrixView< Real2, Segments2, Device, Index2, RealAllocator2, IndexAllocator2 >& matrix, + const RealType& matrixMultiplicator, + const RealType& thisMatrixMultiplicator ) +{ + +} + +template< typename Real, + template< typename, typename > class SegmentsView, + typename Device, + typename Index, + typename RealAllocator, + typename IndexAllocator > +template< typename Real2, typename Index2 > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getTransposition( const SparseSandboxMatrixView< Real2, Device, Index2 >& matrix, + const RealType& matrixMultiplicator ) +{ + +}*/ + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +template< typename Vector1, typename Vector2 > +bool +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +performSORIteration( const Vector1& b, + const IndexType row, + Vector2& x, + const RealType& omega ) const +{ + return false; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +SparseSandboxMatrixView< Real, Device, Index, MatrixType >& +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +operator=( const SparseSandboxMatrixView< Real, Device, Index, MatrixType >& matrix ) +{ + MatrixView< Real, Device, Index >::operator=( matrix ); + this->columnIndexes.bind( matrix.columnIndexes ); + // SANDBOX_TODO: Replace the following line with assignment of metadata required by your + // sparse format. + this->rowPointers.bind( matrix.rowPointers ); + return *this; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Matrix > +bool +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +operator==( const Matrix& m ) const +{ + const auto& view1 = *this; + // FIXME: getConstView does not work + //const auto view2 = m.getConstView(); + const auto view2 = m.getView(); + auto fetch = [=] __cuda_callable__ ( const IndexType i ) -> bool + { + return view1.getRow( i ) == view2.getRow( i ); + }; + return Algorithms::Reduction< DeviceType >::reduce( 0, this->getRows(), fetch, std::logical_and<>{}, true ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > + template< typename Matrix > +bool +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +operator!=( const Matrix& m ) const +{ + return ! operator==( m ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +save( File& file ) const +{ + MatrixView< RealType, DeviceType, IndexType >::save( file ); + file << this->columnIndexes + << this->rowPointers; // SANDBOX_TODO: Replace this with medata required by your format +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +save( const String& fileName ) const +{ + Object::save( fileName ); +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +void +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +print( std::ostream& str ) const +{ + if( isSymmetric() ) + { + for( IndexType row = 0; row < this->getRows(); row++ ) + { + str <<"Row: " << row << " -> "; + for( IndexType column = 0; column < this->getColumns(); column++ ) + { + auto value = this->getElement( row, column ); + if( value != ( RealType ) 0 ) + str << " Col:" << column << "->" << value << "\t"; + } + str << std::endl; + } + } + else + for( IndexType row = 0; row < this->getRows(); row++ ) + { + str <<"Row: " << row << " -> "; + // SANDBOX_TODO: Replace the followinf line with a code for computing number of elements allocated for given matrix row. + const auto rowLength = this->rowPointers.getElement( row + 1 ) - this->rowPointers.getElement( row ); + for( IndexType i = 0; i < rowLength; i++ ) + { + // SANDBOX_TODO: Replace the following line with a code for getting index of the matrix element in arrays `values` and `columnIdexes`. + const IndexType globalIdx = this->rowPointers.getElement( row ) + i; + const IndexType column = this->columnIndexes.getElement( globalIdx ); + if( column == this->getPaddingIndex() ) + break; + RealType value; + if( isBinary() ) + value = ( RealType ) 1.0; + else + value = this->values.getElement( globalIdx ); + if( value ) + { + std::stringstream str_; + str_ << std::setw( 4 ) << std::right << column << ":" << std::setw( 4 ) << std::left << value; + str << std::setw( 10 ) << str_.str(); + } + } + str << std::endl; + } +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +__cuda_callable__ +Index +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getPaddingIndex() const +{ + return -1; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +auto +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getColumnIndexes() const -> const ColumnsIndexesViewType& +{ + return this->columnIndexes; +} + +template< typename Real, + typename Device, + typename Index, + typename MatrixType > +auto +SparseSandboxMatrixView< Real, Device, Index, MatrixType >:: +getColumnIndexes() -> ColumnsIndexesViewType& +{ + return this->columnIndexes; +} + + } // namespace Sandbox + } //namespace Matrices +} // namespace TNL diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp index 40c568b4dfe1ae211a7ee9961718f16d037c2756..b58620c2d0d3c603e92760bc5e094280ede8bef3 100644 --- a/src/TNL/Matrices/SparseMatrix.hpp +++ b/src/TNL/Matrices/SparseMatrix.hpp @@ -417,6 +417,7 @@ SparseMatrix< Real, Device, Index, MatrixType, Segments, ComputeReal, RealAlloca reset() { BaseType::reset(); + this->columnIndexes.reset(); this->segments.reset(); this->view = this->getView(); TNL_ASSERT_EQ( this->getRows(), segments.getSegmentsCount(), "mismatched segments count" ); diff --git a/src/UnitTests/Matrices/CMakeLists.txt b/src/UnitTests/Matrices/CMakeLists.txt index a4b06708e6db2af107bb3b211918eeff2dfb93ad..3210920f58089773cfc798c2f4592037dc4937a6 100644 --- a/src/UnitTests/Matrices/CMakeLists.txt +++ b/src/UnitTests/Matrices/CMakeLists.txt @@ -28,6 +28,8 @@ set( COMMON_TESTS BinarySparseMatrixCopyTest SymmetricSparseMatrixTest_CSR LambdaMatrixTest + SparseMatrixTest_SandboxMatrix + SparseMatrixVectorProductTest_SandboxMatrix ) set( CPP_TESTS diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cpp b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cpp new file mode 100644 index 0000000000000000000000000000000000000000..dc856310e58f2fabf336efbe5c6d950be2998f9f --- /dev/null +++ b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cpp @@ -0,0 +1,11 @@ +/*************************************************************************** + SparseMatrixTest_SandboxMatrix.cpp - description + ------------------- + begin : Apr 19, 2021 + copyright : (C) 2021 by Tomas Oberhuber et al. + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#include "SparseMatrixTest_SandboxMatrix.h" diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cu b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cu new file mode 120000 index 0000000000000000000000000000000000000000..27787fdf2474a56b79948644678b201a81259688 --- /dev/null +++ b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cu @@ -0,0 +1 @@ +SparseMatrixTest_SandboxMatrix.cpp \ No newline at end of file diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.h b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.h new file mode 100644 index 0000000000000000000000000000000000000000..ad1a0c74d6d6f709a63aad6d500a5a103fb04292 --- /dev/null +++ b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.h @@ -0,0 +1,45 @@ +/*************************************************************************** + SandboxMatrixTest_SandboxMatrix.h - description + ------------------- + begin : Apr 19, 2021 + copyright : (C) 2021 by Tomas Oberhuber et al. + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#include <iostream> +#include <TNL/Matrices/Sandbox/SparseSandboxMatrix.h> + +#ifdef HAVE_GTEST +#include <gtest/gtest.h> + +const char* saveAndLoadFileName = "test_SandboxMatrixTest_SandboxMatrix"; + +// types for which MatrixTest is instantiated +using MatrixTypes = ::testing::Types +< + TNL::Matrices::Sandbox::SparseSandboxMatrix< int, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< long, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< float, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< double, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< int, TNL::Devices::Host, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< long, TNL::Devices::Host, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< float, TNL::Devices::Host, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< double, TNL::Devices::Host, long, TNL::Matrices::GeneralMatrix > +#ifdef HAVE_CUDA + ,TNL::Matrices::Sandbox::SparseSandboxMatrix< int, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< long, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< float, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< double, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< int, TNL::Devices::Cuda, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< long, TNL::Devices::Cuda, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< float, TNL::Devices::Cuda, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< double, TNL::Devices::Cuda, long, TNL::Matrices::GeneralMatrix > +#endif +>; + +#endif + +#include "SparseMatrixTest.h" +#include "../main.h" diff --git a/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cpp b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cpp new file mode 100644 index 0000000000000000000000000000000000000000..bfa16c02b91c7cddb03a63763e994ade31bcb0f8 --- /dev/null +++ b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cpp @@ -0,0 +1,11 @@ +/*************************************************************************** + SparseMatrixVectorProductTest_CSRScalar.cpp - description + ------------------- + begin : Mar 30, 2021 + copyright : (C) 2021 by Tomas Oberhuber et al. + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#include "SparseMatrixVectorProductTest_CSRScalar.h" diff --git a/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cu b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cu new file mode 120000 index 0000000000000000000000000000000000000000..bd87e1ad0c244dd2117317355391a5d450f8de98 --- /dev/null +++ b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cu @@ -0,0 +1 @@ +SparseMatrixVectorProductTest_SandboxMatrix.cpp \ No newline at end of file diff --git a/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.h b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.h new file mode 100644 index 0000000000000000000000000000000000000000..7b06af0f3a8049d22e81fa480c59c21d106696c0 --- /dev/null +++ b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.h @@ -0,0 +1,45 @@ +/*************************************************************************** + SparseMatrixVectorProductTest_SandbxMatrix.h - description + ------------------- + begin : Apr 22, 2021 + copyright : (C) 2021 by Tomas Oberhuber et al. + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#include <iostream> +#include <TNL/Matrices/Sandbox/SparseSandboxMatrix.h> + +#ifdef HAVE_GTEST +#include <gtest/gtest.h> + +const char* saveAndLoadFileName = "test_SparseMatrixTest_CSRScalar_segments"; + +// types for which MatrixTest is instantiated +using MatrixTypes = ::testing::Types +< + TNL::Matrices::Sandbox::SparseSandboxMatrix< int, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< long, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< float, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< double, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< int, TNL::Devices::Host, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< long, TNL::Devices::Host, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< float, TNL::Devices::Host, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< double, TNL::Devices::Host, long, TNL::Matrices::GeneralMatrix > +#ifdef HAVE_CUDA + ,TNL::Matrices::Sandbox::SparseSandboxMatrix< int, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< long, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< float, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< double, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< int, TNL::Devices::Cuda, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< long, TNL::Devices::Cuda, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< float, TNL::Devices::Cuda, long, TNL::Matrices::GeneralMatrix >, + TNL::Matrices::Sandbox::SparseSandboxMatrix< double, TNL::Devices::Cuda, long, TNL::Matrices::GeneralMatrix > +#endif +>; + +#endif + +#include "SparseMatrixVectorProductTest.h" +#include "../main.h"