diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index 03273c98b2128dfff41feeda0bed1aafc9c54536..3eccc7211dc37b6d42178aff66a512a020450e3e 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -11,6 +11,7 @@
#pragma once
#include <functional>
+#include <sstream>
#include <TNL/Algorithms/Reduction.h>
#include <TNL/Matrices/SparseMatrix.h>
@@ -107,11 +108,27 @@ SparseMatrix( const IndexType rows,
{
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 ) ]++;
+ }
SparseMatrix< Real, Devices::Host, Index, MatrixType, Segments > hostMatrix( rows, columns );
hostMatrix.setCompressedRowLengths( 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;
}
diff --git a/src/TNL/Matrices/SparseMatrixView.h b/src/TNL/Matrices/SparseMatrixView.h
index 2756c80d7d2d3806c8e93ddfae12583cbeaabd5f..54d4f17660c85f08c2f7b4c74b46a180e588ff08 100644
--- a/src/TNL/Matrices/SparseMatrixView.h
+++ b/src/TNL/Matrices/SparseMatrixView.h
@@ -98,8 +98,8 @@ class SparseMatrixView : public MatrixView< Real, Device, Index >
const IndexType column,
const RealType& value );
- void addElement( const IndexType row,
- const IndexType column,
+ void addElement( IndexType row,
+ IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator = 1.0 );
diff --git a/src/TNL/Matrices/SparseMatrixView.hpp b/src/TNL/Matrices/SparseMatrixView.hpp
index 16a8bc62f355c0331809951974cee569bfdaa175..62300217f83d32ff3d8734a450e394f4a522faa7 100644
--- a/src/TNL/Matrices/SparseMatrixView.hpp
+++ b/src/TNL/Matrices/SparseMatrixView.hpp
@@ -153,10 +153,36 @@ getNumberOfNonzeroMatrixElements() const
{
const auto columns_view = this->columnIndexes.getConstView();
const IndexType paddingIndex = this->getPaddingIndex();
- auto fetch = [=] __cuda_callable__ ( const IndexType i ) -> IndexType {
- return ( columns_view[ i ] != paddingIndex );
- };
- return Algorithms::Reduction< DeviceType >::reduce( this->columnIndexes.getSize(), std::plus<>{}, fetch, 0 );
+ if( ! isSymmetric() )
+ {
+ auto fetch = [=] __cuda_callable__ ( const IndexType i ) -> IndexType {
+ return ( columns_view[ i ] != paddingIndex );
+ };
+ return Algorithms::Reduction< DeviceType >::reduce( this->columnIndexes.getSize(), std::plus<>{}, fetch, 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();
+ const auto columnIndexesView = this->columnIndexes.getConstView();
+ auto fetch = [=] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType globalIdx, bool& compute ) -> RealType {
+ const IndexType column = columnIndexesView[ globalIdx ];
+ compute = ( column != paddingIndex );
+ if( ! compute )
+ return 0.0;
+ return 1 + ( column != row && column < rows && row < columns ); // the addition is for non-diagonal elements
+ };
+ auto reduction = [] __cuda_callable__ ( RealType& sum, const RealType& value ) {
+ sum += value;
+ };
+ auto keeper = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable {
+ row_sums_view[ row ] = value;
+ };
+ this->segments.segmentsReduction( 0, this->getRows(), fetch, reduction, keeper, ( RealType ) 0.0 );
+ return sum( row_sums );
+ }
}
template< typename Real,
@@ -206,8 +232,8 @@ template< typename Real,
template< typename, typename > class SegmentsView >
void
SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView >::
-addElement( const IndexType row,
- const IndexType column,
+addElement( IndexType row,
+ IndexType column,
const RealType& value,
const RealType& thisElementMultiplicator )
{
@@ -216,6 +242,13 @@ addElement( const IndexType row,
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." );
+ }
+
const IndexType rowSize = this->segments.getSegmentSize( row );
IndexType col( this->getPaddingIndex() );
IndexType i;
@@ -276,14 +309,21 @@ template< typename Real,
template< typename, typename > class SegmentsView >
Real
SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView >::
-getElement( const IndexType row,
- const IndexType column ) const
+getElement( IndexType row,
+ IndexType column ) const
{
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;
+ }
+
const IndexType rowSize = this->segments.getSegmentSize( row );
for( IndexType i = 0; i < rowSize; i++ )
{
@@ -588,25 +628,40 @@ void
SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView >::
print( std::ostream& str ) const
{
- for( IndexType row = 0; row < this->getRows(); row++ )
+ if( isSymmetric() )
{
- str <<"Row: " << row << " -> ";
- const IndexType rowLength = this->segments.getSegmentSize( row );
- for( IndexType i = 0; i < rowLength; i++ )
+ for( IndexType row = 0; row < this->getRows(); row++ )
{
- const IndexType globalIdx = this->segments.getGlobalIndex( row, i );
- const IndexType column = this->columnIndexes.getElement( globalIdx );
- if( column == this->getPaddingIndex() )
- break;
- RealType value;
- if( isBinary() )
- value = 1.0;
- else
- value = this->values.getElement( globalIdx );
- str << " Col:" << column << "->" << value << "\t";
+ str <<"Row: " << row << " -> ";
+ for( IndexType column = 0; column < this->getColumns(); column++ )
+ {
+ auto value = this->getElement( row, column );
+ if( value )
+ str << " Col:" << column << "->" << value << "\t";
+ }
+ str << std::endl;
}
- str << std::endl;
}
+ else
+ for( IndexType row = 0; row < this->getRows(); row++ )
+ {
+ str <<"Row: " << row << " -> ";
+ const auto rowLength = this->segments.getSegmentSize( row );
+ for( IndexType i = 0; i < rowLength; i++ )
+ {
+ const IndexType globalIdx = this->segments.getGlobalIndex( 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 );
+ str << " Col:" << column << "->" << value << "\t";
+ }
+ str << std::endl;
+ }
}
template< typename Real,
diff --git a/src/UnitTests/Matrices/SymmetricSparseMatrixTest.h b/src/UnitTests/Matrices/SymmetricSparseMatrixTest.h
index 5582b138ddf3173ecbde7ecf9cee9e73a6acf69c..02fd8c585366f4da12d1218a28adca717dd2cdf2 100644
--- a/src/UnitTests/Matrices/SymmetricSparseMatrixTest.h
+++ b/src/UnitTests/Matrices/SymmetricSparseMatrixTest.h
@@ -45,6 +45,13 @@ TYPED_TEST( MatrixTest, setLikeTest )
test_SetLike< MatrixType, MatrixType >();
}
+TYPED_TEST( MatrixTest, getNumberOfNonzeroMatrixElements )
+{
+ using MatrixType = typename TestFixture::MatrixType;
+
+ test_GetNumberOfNonzeroMatrixElements< MatrixType >();
+}
+
TYPED_TEST( MatrixTest, resetTest )
{
using MatrixType = typename TestFixture::MatrixType;
diff --git a/src/UnitTests/Matrices/SymmetricSparseMatrixTest.hpp b/src/UnitTests/Matrices/SymmetricSparseMatrixTest.hpp
index 75b121060a93c760442d7a4dc61e909c357a6630..7c82784228195ee372bd14b63ea73574c6359240 100644
--- a/src/UnitTests/Matrices/SymmetricSparseMatrixTest.hpp
+++ b/src/UnitTests/Matrices/SymmetricSparseMatrixTest.hpp
@@ -73,11 +73,11 @@ void test_SetCompressedRowLengths()
| 24 25 26 |
\ 27 28 30 /
*/
- const IndexType rows = 10;
+ const IndexType rows = 11;
const IndexType cols = 11;
Matrix m( rows, cols );
- typename Matrix::CompressedRowLengthsVector rowLengths { 1, 2, 3, 3, 3, 3, 3, 3, 3, 3 };
+ typename Matrix::CompressedRowLengthsVector rowLengths { 1, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3 };
m.setCompressedRowLengths( rowLengths );
// Insert values into the rows.
@@ -120,9 +120,25 @@ void test_SetCompressedRowLengths()
m.setElement( 7, 3, value++ );
m.setElement( 7, 7, value++ );
+ // 8th row - lower part
+ m.setElement( 8, 2, value++ );
+ m.setElement( 8, 4, value++ );
+ m.setElement( 8, 8, value++ );
+
+ // 8th row - lower part
+ m.setElement( 9, 2, value++ );
+ m.setElement( 9, 4, value++ );
+ m.setElement( 9, 9, value++ );
+
+ // 8th row - lower part
+ m.setElement( 10, 2, value++ );
+ m.setElement( 10, 5, value++ );
+ m.setElement( 10, 10, value++ );
+
rowLengths = 0;
m.getCompressedRowLengths( rowLengths );
- typename Matrix::CompressedRowLengthsVector correctRowLengths{ 1, 2, 3, 3, 3, 3, 3, 3, 3, 3 };
+
+ typename Matrix::CompressedRowLengthsVector correctRowLengths{ 1, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3 };
EXPECT_EQ( rowLengths, correctRowLengths );
}
@@ -170,8 +186,8 @@ void test_GetNumberOfNonzeroMatrixElements()
49
*/
- const IndexType rows = 10;
- const IndexType cols = 10;
+ const IndexType rows = 11;
+ const IndexType cols = 11;
Matrix m( rows, cols, {
{ 0, 0, 1 },
@@ -276,7 +292,7 @@ void test_GetRow()
EXPECT_EQ( m.getElement( 1, 0 ), 2 );
EXPECT_EQ( m.getElement( 1, 1 ), 3 );
- EXPECT_EQ( m.getElement( 1, 2 ), 4 );
+ EXPECT_EQ( m.getElement( 1, 2 ), 5 );
EXPECT_EQ( m.getElement( 1, 3 ), 8 );
EXPECT_EQ( m.getElement( 1, 4 ), 10 );
EXPECT_EQ( m.getElement( 1, 5 ), 13 );
@@ -421,7 +437,7 @@ void test_SetElement()
Matrix m( { 1, 1, 1, 4, 1, 1, 7, 1, 1, 1 }, 10 );
RealType value = 1;
- for( IndexType i = 0; i < 4; i++ )
+ for( IndexType i = 0; i < 3; i++ )
m.setElement( i, i, value++ );
for( IndexType i = 0; i < 4; i++ )
@@ -574,7 +590,7 @@ void test_AddElement()
{ 2, 1, 4 }, { 2, 2, 5 },
{ 3, 2, 6 }, { 3, 3, 7 },
{ 4, 3, 8 }, { 4, 4, 9 },
- { 5, 5, 10 } } );
+ { 5, 4, 10 } } );
// Check the set elements
EXPECT_EQ( m.getElement( 0, 0 ), 1 );
@@ -626,9 +642,12 @@ void test_AddElement()
*/
for( IndexType i = 1; i < rows; i++ )
+ {
m.addElement( i, i - 1, 1.0, 2.0 );
+ m.addElement( i, i, 0.0, 2.0 );
+ }
-
+ std::cerr << m << std::endl;
EXPECT_EQ( m.getElement( 0, 0 ), 2 );
EXPECT_EQ( m.getElement( 0, 1 ), 5 );
EXPECT_EQ( m.getElement( 0, 2 ), 0 );