Implementintg symmetric sparse matrix. (537a8805) · Commits · TNL / tnl-dev

src/TNL/Matrices/SparseMatrix.hpp

+17 −0

Original line number	Diff line number	Diff line
		@@ -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;
		}

src/TNL/Matrices/SparseMatrixView.h

+2 −2

Original line number	Diff line number	Diff line
		@@ -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 );

src/TNL/Matrices/SparseMatrixView.hpp

+78 −23

Original line number	Diff line number	Diff line
		@@ -153,11 +153,37 @@ getNumberOfNonzeroMatrixElements() 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( 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,
		typename Device,
		@@ -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,10 +628,25 @@ void
		SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView >::
		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 )
		str << " Col:" << column << "->" << value << "\t";
		}
		str << std::endl;
		}
		}
		else
		for( IndexType row = 0; row < this->getRows(); row++ )
		{
		str <<"Row: " << row << " -> ";
		const IndexType rowLength = this->segments.getSegmentSize( row );
		const auto rowLength = this->segments.getSegmentSize( row );
		for( IndexType i = 0; i < rowLength; i++ )
		{
		const IndexType globalIdx = this->segments.getGlobalIndex( row, i );
		@@ -600,7 +655,7 @@ print( std::ostream& str ) const
		break;
		RealType value;
		if( isBinary() )
		value = 1.0;
		value = ( RealType ) 1.0;
		else
		value = this->values.getElement( globalIdx );
		str << " Col:" << column << "->" << value << "\t";

src/UnitTests/Matrices/SymmetricSparseMatrixTest.h

+7 −0

Original line number	Diff line number	Diff line
		@@ -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;

src/UnitTests/Matrices/SymmetricSparseMatrixTest.hpp

+28 −9

Original line number	Diff line number	Diff line
		@@ -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 );