diff --git a/src/TNL/Matrices/SparseMatrix.h b/src/TNL/Matrices/SparseMatrix.h
index b510636f55d101492552fce1f2830d898f6c26cd..268fba6d39a7b26dc36431fbac63a318a7796782 100644
--- a/src/TNL/Matrices/SparseMatrix.h
+++ b/src/TNL/Matrices/SparseMatrix.h
@@ -34,10 +34,17 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
using IndexType = Index;
using RealAllocatorType = RealAllocator;
using IndexAllocatorType = IndexAllocator;
- using CompressedRowLengthsVectorView = Containers::VectorView< IndexType, DeviceType, IndexType >;
- using ConstCompressedRowLengthsVectorView = typename CompressedRowLengthsVectorView::ConstViewType;
+ using RowsCapacitiesType = Containers::Vector< IndexType, DeviceType, IndexType >;
+ using RowsCapacitiesView = Containers::VectorView< IndexType, DeviceType, IndexType >;
+ using ConstRowsCapacitiesView = typename RowsCapacitiesView::ConstViewType;
using ValuesVectorType = typename Matrix< Real, Device, Index, RealAllocator >::ValuesVector;
using ColumnsVectorType = Containers::Vector< IndexType, DeviceType, IndexType, IndexAllocatorType >;
+
+ // TODO: remove this - it is here only for compatibility with original matrix implementation
+ typedef Containers::Vector< IndexType, DeviceType, IndexType > CompressedRowLengthsVector;
+ typedef Containers::VectorView< IndexType, DeviceType, IndexType > CompressedRowLengthsVectorView;
+ typedef typename CompressedRowLengthsVectorView::ConstViewType ConstCompressedRowLengthsVectorView;
+
SparseMatrix( const RealAllocatorType& realAllocator = RealAllocatorType(),
const IndexAllocatorType& indexAllocator = IndexAllocatorType() );
@@ -158,6 +165,12 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
const RealType& matrixMultiplicator = 1.0 );
*/
+ template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+ void rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const;
+
+ template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+ void allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const;
+
template< typename Vector1, typename Vector2 >
bool performSORIteration( const Vector1& b,
const IndexType row,
diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index 1c243bcea886d27b093daef4a1befe948953398a..e26693f6a12953d59a5950b9def449121b68de75 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -517,6 +517,44 @@ vectorProduct( const InVector& inVector,
this->segments.segmentsReduction( 0, this->getRows(), fetch, reduction, keeper, ( RealType ) 0.0 );
}
+template< typename Real,
+ template< typename, typename > class Segments,
+ typename Device,
+ typename Index,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Fetch, typename Reduce, typename Keep, typename FetchValue >
+void
+SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchValue& zero ) const
+{
+ const auto columns_view = this->columnIndexes.getConstView();
+ const auto values_view = this->values.getConstView();
+ const IndexType paddingIndex_ = this->getPaddingIndex();
+ auto fetch_ = [=] __cuda_callable__ ( IndexType rowIdx, IndexType globalIdx ) mutable -> decltype( fetch( IndexType(), IndexType(), RealType() ) ) {
+ IndexType columnIdx = columns_view[ globalIdx ];
+ if( columnIdx != paddingIndex_ )
+ return fetch( rowIdx, columnIdx, values_view[ globalIdx ] );
+ return zero;
+ };
+ this->segments.segmentsReduction( first, last, fetch_, reduce, keep, zero );
+}
+
+template< typename Real,
+ template< typename, typename > class Segments,
+ typename Device,
+ typename Index,
+ typename RealAllocator,
+ typename IndexAllocator >
+ template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
+{
+ this->rowsReduction( 0, this->getRows(), fetch, reduce, keep, zero );
+}
+
+
/*template< typename Real,
template< typename, typename > class Segments,
typename Device,
@@ -576,7 +614,11 @@ SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >&
SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
operator=( const SparseMatrix& matrix )
{
-
+ Matrix< Real, Device, Index >::operator=( matrix );
+ this->columnIndexes = matrix.columnIndexes;
+ this->segments = matrix.segments;
+ this->indexAlloctor = matrix.indexAllocator;
+ this->realAllocator = matrix.realAllocator;
}
// cross-device copy assignment
@@ -596,7 +638,12 @@ SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >&
SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
operator=( const SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2, IndexAllocator2 >& matrix )
{
-
+ if( std::is_same< Device, Device2 >::value )
+ {
+
+ }
+
+
}
template< typename Real,
diff --git a/src/UnitTests/Matrices/SparseMatrixTest.hpp b/src/UnitTests/Matrices/SparseMatrixTest.hpp
index 5dcd96ebc7c5d0cac1e6a7f1c5b031f27efe79a6..b366f4e2fea5b8688ee89fed3cf4d6111cbf484b 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest.hpp
+++ b/src/UnitTests/Matrices/SparseMatrixTest.hpp
@@ -908,18 +908,18 @@ void test_VectorProduct()
EXPECT_EQ( outVector_4.getElement( 7 ), 330 );
-/*
- * Sets up the following 8x8 sparse matrix:
- *
- * / 1 2 3 0 4 5 0 1 \ 6
- * | 0 6 0 7 0 0 0 1 | 3
- * | 0 8 9 0 10 0 0 1 | 4
- * | 0 11 12 13 14 0 0 1 | 5
- * | 0 15 0 0 0 0 0 1 | 2
- * | 0 16 17 18 19 20 21 1 | 7
- * | 22 23 24 25 26 27 28 1 | 8
- * \ 29 30 31 32 33 34 35 36 / 8
- */
+ /*
+ * Sets up the following 8x8 sparse matrix:
+ *
+ * / 1 2 3 0 4 5 0 1 \ 6
+ * | 0 6 0 7 0 0 0 1 | 3
+ * | 0 8 9 0 10 0 0 1 | 4
+ * | 0 11 12 13 14 0 0 1 | 5
+ * | 0 15 0 0 0 0 0 1 | 2
+ * | 0 16 17 18 19 20 21 1 | 7
+ * | 22 23 24 25 26 27 28 1 | 8
+ * \ 29 30 31 32 33 34 35 36 / 8
+ */
const IndexType m_rows_5 = 8;
const IndexType m_cols_5 = 8;
@@ -970,20 +970,18 @@ void test_VectorProduct()
for( IndexType i = 0; i < 7; i++ ) // 1s at the end of rows
m_5.setElement( i, 7, 1);
-
+
VectorType inVector_5;
inVector_5.setSize( m_cols_5 );
- for( IndexType i = 0; i < inVector_5.getSize(); i++ )
+ for( IndexType i = 0; i < inVector_5.getSize(); i++ )
inVector_5.setElement( i, 2 );
VectorType outVector_5;
outVector_5.setSize( m_rows_5 );
for( IndexType j = 0; j < outVector_5.getSize(); j++ )
outVector_5.setElement( j, 0 );
-
-
+
m_5.vectorProduct( inVector_5, outVector_5 );
-
EXPECT_EQ( outVector_5.getElement( 0 ), 32 );
EXPECT_EQ( outVector_5.getElement( 1 ), 28 );
@@ -995,6 +993,109 @@ void test_VectorProduct()
EXPECT_EQ( outVector_5.getElement( 7 ), 520 );
}
+template< typename Matrix >
+void test_RowsReduction()
+{
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
+
+ /*
+ * Sets up the following 8x8 sparse matrix:
+ *
+ * / 1 2 3 0 4 5 0 1 \ 6
+ * | 0 6 0 7 0 0 0 1 | 3
+ * | 0 8 9 0 10 0 0 1 | 4
+ * | 0 11 12 13 14 0 0 1 | 5
+ * | 0 15 0 0 0 0 0 1 | 2
+ * | 0 16 17 18 19 20 21 1 | 7
+ * | 22 23 24 25 26 27 28 1 | 8
+ * \ 29 30 31 32 33 34 35 36 / 8
+ */
+
+ const IndexType rows = 8;
+ const IndexType cols = 8;
+
+ Matrix m;
+ m.setDimensions( rows, cols );
+ typename Matrix::RowsCapacitiesType rowsCapacities( rows );
+ //rowLengths.setSize( rows );
+ rowsCapacities.setElement(0, 6);
+ rowsCapacities.setElement(1, 3);
+ rowsCapacities.setElement(2, 4);
+ rowsCapacities.setElement(3, 5);
+ rowsCapacities.setElement(4, 2);
+ rowsCapacities.setElement(5, 7);
+ rowsCapacities.setElement(6, 8);
+ rowsCapacities.setElement(7, 8);
+ m.setCompressedRowLengths( rowsCapacities );
+
+ RealType value = 1;
+ for( IndexType i = 0; i < 3; i++ ) // 0th row
+ m.setElement( 0, i, value++ );
+
+ m.setElement( 0, 4, value++ ); // 0th row
+ m.setElement( 0, 5, value++ );
+
+ m.setElement( 1, 1, value++ ); // 1st row
+ m.setElement( 1, 3, value++ );
+
+ for( IndexType i = 1; i < 3; i++ ) // 2nd row
+ m.setElement( 2, i, value++ );
+
+ m.setElement( 2, 4, value++ ); // 2nd row
+
+ for( IndexType i = 1; i < 5; i++ ) // 3rd row
+ m.setElement( 3, i, value++ );
+
+ m.setElement( 4, 1, value++ ); // 4th row
+
+ for( IndexType i = 1; i < 7; i++ ) // 5th row
+ m.setElement( 5, i, value++ );
+
+ for( IndexType i = 0; i < 7; i++ ) // 6th row
+ m.setElement( 6, i, value++ );
+
+ for( IndexType i = 0; i < 8; i++ ) // 7th row
+ m.setElement( 7, i, value++ );
+
+ for( IndexType i = 0; i < 7; i++ ) // 1s at the end of rows
+ m.setElement( i, 7, 1);
+
+ ////
+ // Compute number of non-zero elements in rows.
+ typename Matrix::RowsCapacitiesType rowLengths( rows );
+ auto rowLengths_view = rowLengths.getView();
+ auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+ return ( value != 0.0 );
+ };
+ auto reduce = [] __cuda_callable__ ( IndexType& aux, const IndexType a ) {
+ aux += a;
+ };
+ auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+ rowLengths_view[ rowIdx ] = value;
+ };
+ m.allRowsReduction( fetch, reduce, keep, 0 );
+ EXPECT_EQ( rowsCapacities, rowLengths );
+
+ ////
+ // Compute max norm
+ TNL::Containers::Vector< RealType, DeviceType, IndexType > rowSums( rows );
+ auto rowSums_view = rowSums.getView();
+ auto max_fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+ return abs( value );
+ };
+ auto max_reduce = [] __cuda_callable__ ( IndexType& aux, const IndexType a ) {
+ aux += a;
+ };
+ auto max_keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+ rowSums_view[ rowIdx ] = value;
+ };
+ m.allRowsReduction( max_fetch, max_reduce, max_keep, 0 );
+ const RealType maxNorm = TNL::max( rowSums );
+ EXPECT_EQ( maxNorm, 260 ) ; // 29+30+31+32+33+34+35+36
+}
+
template< typename Matrix >
void test_PerformSORIteration()
{
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_CSR_segments.h b/src/UnitTests/Matrices/SparseMatrixTest_CSR_segments.h
index bf4e452fa1af77388efb1ba2da76b5843d04e067..0718e3a69a8978d322288d9ee65ec1a471b6ee22 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_CSR_segments.h
+++ b/src/UnitTests/Matrices/SparseMatrixTest_CSR_segments.h
@@ -122,6 +122,13 @@ TYPED_TEST( CSRMatrixTest, vectorProductTest )
test_VectorProduct< CSRMatrixType >();
}
+TYPED_TEST( CSRMatrixTest, rowsReduction )
+{
+ using CSRMatrixType = typename TestFixture::CSRMatrixType;
+
+ test_RowsReduction< CSRMatrixType >();
+}
+
TYPED_TEST( CSRMatrixTest, saveAndLoadTest )
{
using CSRMatrixType = typename TestFixture::CSRMatrixType;
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_Ellpack_segments.h b/src/UnitTests/Matrices/SparseMatrixTest_Ellpack_segments.h
index edfe0bc2839383a08c9fa65bd69c0c733753abab..2c0514c0aba69e5e26abdf1434b7a2d333c25b63 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_Ellpack_segments.h
+++ b/src/UnitTests/Matrices/SparseMatrixTest_Ellpack_segments.h
@@ -133,6 +133,13 @@ TYPED_TEST( EllpackMatrixTest, vectorProductTest )
test_VectorProduct< EllpackMatrixType >();
}
+TYPED_TEST( EllpackMatrixTest, rowsReduction )
+{
+ using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+
+ test_RowsReduction< EllpackMatrixType >();
+}
+
TYPED_TEST( EllpackMatrixTest, saveAndLoadTest )
{
using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack_segments.h b/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack_segments.h
index 8d17b8be7a4fed6e14c97e3848e891ef1bbd4e19..5efcb1eae0a01e9d4267c8e07cdcd396ecb1015e 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack_segments.h
+++ b/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack_segments.h
@@ -133,6 +133,13 @@ TYPED_TEST( SlicedEllpackMatrixTest, vectorProductTest )
test_VectorProduct< SlicedEllpackMatrixType >();
}
+TYPED_TEST( SlicedEllpackMatrixTest, rowsReduction )
+{
+ using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+
+ test_RowsReduction< SlicedEllpackMatrixType >();
+}
+
TYPED_TEST( SlicedEllpackMatrixTest, saveAndLoadTest )
{
using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;