diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h
index 27a9dfaa83fb06c2539269e8f21b1f61fa5585ca..2e4d55643cc7995cc1fcb7a9e2f9445d286b9deb 100644
--- a/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h
+++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h
@@ -1054,11 +1054,12 @@ public:
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 );
- */
+ void
+ getTransposition( const SparseSandboxMatrix< Real2, Device, Index2, MatrixType >& matrix,
+ const RealType& matrixMultiplicator = 1.0 );
template< typename Vector1, typename Vector2 >
bool
diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp
index a5fc69142ab937eb3c4335c26b31032d697d93fe..f84f79f920953136d6b0cb96f2841e1176195536 100644
--- a/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp
+++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp
@@ -521,7 +521,8 @@ SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAlloca
this->sequentialForRows( 0, this->getRows(), function );
}
-/*template< typename Real,
+/*
+template< typename Real,
template< typename, typename, typename > class Segments,
typename Device,
typename Index,
@@ -532,23 +533,67 @@ IndexAllocator2 > void SparseSandboxMatrix< Real, Device, Index, MatrixType, Rea
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 Real, typename Device, typename Index, typename MatrixType, 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 )
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::getTransposition(
+ const SparseSandboxMatrix< Real2, Device, Index2, MatrixType >& matrix,
+ const RealType& matrixMultiplicator )
{
+ // set transposed dimensions
+ setDimensions( matrix.getColumns(), matrix.getRows() );
-}*/
+ // stage 1: compute row capacities for the transposition
+ RowsCapacitiesType capacities;
+ capacities.resize( this->getRows(), 0 );
+ auto capacities_view = capacities.getView();
+ using MatrixRowView = typename SparseSandboxMatrix< Real2, Device, Index2, MatrixType >::ConstRowView;
+ matrix.forAllRows(
+ [ = ] __cuda_callable__( const MatrixRowView& row ) mutable
+ {
+ for( IndexType c = 0; c < row.getSize(); c++ ) {
+ // row index of the transpose = column index of the input
+ const IndexType& transRowIdx = row.getColumnIndex( c );
+ if( transRowIdx < 0 )
+ continue;
+ // increment the capacity for the row in the transpose
+ Algorithms::AtomicOperations< DeviceType >::add( capacities_view[ row.getColumnIndex( c ) ], IndexType( 1 ) );
+ }
+ } );
+
+ // set the row capacities
+ setRowCapacities( capacities );
+ capacities.reset();
+
+ // index of the first unwritten element per row
+ RowsCapacitiesType offsets;
+ offsets.resize( this->getRows(), 0 );
+ auto offsets_view = offsets.getView();
+
+ // stage 2: copy and transpose the data
+ auto trans_view = getView();
+ matrix.forAllRows(
+ [ = ] __cuda_callable__( const MatrixRowView& row ) mutable
+ {
+ // row index of the input = column index of the transpose
+ const IndexType& rowIdx = row.getRowIndex();
+ for( IndexType c = 0; c < row.getSize(); c++ ) {
+ // row index of the transpose = column index of the input
+ const IndexType& transRowIdx = row.getColumnIndex( c );
+ if( transRowIdx < 0 )
+ continue;
+ // local index in the row of the transpose
+ const IndexType transLocalIdx =
+ Algorithms::AtomicOperations< DeviceType >::add( offsets_view[ transRowIdx ], IndexType( 1 ) );
+ // get the row in the transposed matrix and set the value
+ auto transRow = trans_view.getRow( transRowIdx );
+ transRow.setElement( transLocalIdx, rowIdx, row.getValue( c ) * matrixMultiplicator );
+ }
+ } );
+}
template< typename Real, typename Device, typename Index, typename MatrixType, typename RealAllocator, typename IndexAllocator >
template< typename Vector1, typename Vector2 >
diff --git a/src/TNL/Matrices/SparseMatrix.h b/src/TNL/Matrices/SparseMatrix.h
index 76aee6dfb4498c782f9d99c6281e9c264749b5b8..b47ef9c853d991774ac384e2f5cf2f2fb25466e8 100644
--- a/src/TNL/Matrices/SparseMatrix.h
+++ b/src/TNL/Matrices/SparseMatrix.h
@@ -1088,12 +1088,13 @@ public:
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 Real2, typename Index2, template< typename, typename, typename > class Segments2 >
+ void
+ getTransposition( const SparseMatrix< Real2, Device, Index2, MatrixType, Segments2 >& matrix,
+ const ComputeRealType& matrixMultiplicator = 1.0 );
+
/**
* \brief Copy-assignment of exactly the same matrix type.
*
diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index 9892b948c85f9f8333d76e15da6e5f2456ba4d18..a2d497f1e51cb96c9a9bfb4bb49bbbe631ecb8f3 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -870,7 +870,8 @@ SparseMatrix< Real, Device, Index, MatrixType, Segments, ComputeReal, RealAlloca
this->sequentialForRows( (IndexType) 0, this->getRows(), function );
}
-/*template< typename Real,
+/*
+template< typename Real,
template< typename, typename, typename > class Segments,
typename Device,
typename Index,
@@ -882,23 +883,75 @@ addMatrix( const SparseMatrix< Real2, Segments2, Device, Index2, RealAllocator2,
const RealType& matrixMultiplicator,
const RealType& thisMatrixMultiplicator )
{
-
}
+*/
template< typename Real,
- template< typename, typename, typename > class Segments,
typename Device,
typename Index,
+ typename MatrixType,
+ template< typename, typename, typename >
+ class Segments,
+ typename ComputeReal,
typename RealAllocator,
typename IndexAllocator >
-template< typename Real2, typename Index2 >
+template< typename Real2, typename Index2, template< typename, typename, typename > class Segments2 >
void
-SparseMatrix< Real, Device, Index, MatrixType, Segments, ComputeReal, RealAllocator, IndexAllocator >::
-getTransposition( const SparseMatrix< Real2, Device, Index2 >& matrix,
- const RealType& matrixMultiplicator )
+SparseMatrix< Real, Device, Index, MatrixType, Segments, ComputeReal, RealAllocator, IndexAllocator >::getTransposition(
+ const SparseMatrix< Real2, Device, Index2, MatrixType, Segments2 >& matrix,
+ const ComputeRealType& matrixMultiplicator )
{
+ // set transposed dimensions
+ setDimensions( matrix.getColumns(), matrix.getRows() );
+
+ // stage 1: compute row capacities for the transposition
+ RowsCapacitiesType capacities;
+ capacities.resize( this->getRows(), 0 );
+ auto capacities_view = capacities.getView();
+ using MatrixRowView = typename SparseMatrix< Real2, Device, Index2, MatrixType, Segments2 >::ConstRowView;
+ matrix.forAllRows(
+ [ = ] __cuda_callable__( const MatrixRowView& row ) mutable
+ {
+ for( IndexType c = 0; c < row.getSize(); c++ ) {
+ // row index of the transpose = column index of the input
+ const IndexType& transRowIdx = row.getColumnIndex( c );
+ if( transRowIdx == row.getPaddingIndex() )
+ continue;
+ // increment the capacity for the row in the transpose
+ Algorithms::AtomicOperations< DeviceType >::add( capacities_view[ transRowIdx ], IndexType( 1 ) );
+ }
+ } );
+
+ // set the row capacities
+ setRowCapacities( capacities );
+ capacities.reset();
+
+ // index of the first unwritten element per row
+ RowsCapacitiesType offsets;
+ offsets.resize( this->getRows(), 0 );
+ auto offsets_view = offsets.getView();
-}*/
+ // stage 2: copy and transpose the data
+ auto trans_view = getView();
+ matrix.forAllRows(
+ [ = ] __cuda_callable__( const MatrixRowView& row ) mutable
+ {
+ // row index of the input = column index of the transpose
+ const IndexType& rowIdx = row.getRowIndex();
+ for( IndexType c = 0; c < row.getSize(); c++ ) {
+ // row index of the transpose = column index of the input
+ const IndexType& transRowIdx = row.getColumnIndex( c );
+ if( transRowIdx == row.getPaddingIndex() )
+ continue;
+ // local index in the row of the transpose
+ const IndexType transLocalIdx =
+ Algorithms::AtomicOperations< DeviceType >::add( offsets_view[ transRowIdx ], IndexType( 1 ) );
+ // get the row in the transposed matrix and set the value
+ auto transRow = trans_view.getRow( transRowIdx );
+ transRow.setElement( transLocalIdx, rowIdx, row.getValue( c ) * matrixMultiplicator );
+ }
+ } );
+}
// copy assignment
template< typename Real,
diff --git a/src/TNL/Matrices/SparseMatrixView.hpp b/src/TNL/Matrices/SparseMatrixView.hpp
index 9420363b82be57eee4995bffa783ce02fac09114..9ee552fbd6e88be04f9247148690707a78f22f1c 100644
--- a/src/TNL/Matrices/SparseMatrixView.hpp
+++ b/src/TNL/Matrices/SparseMatrixView.hpp
@@ -820,35 +820,6 @@ SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView, ComputeReal >::
this->sequentialForRows( (IndexType) 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 SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView, ComputeReal >:: addMatrix( const
-SparseMatrixView< 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
-SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView, ComputeReal >::
-getTransposition( const SparseMatrixView< Real2, Device, Index2 >& matrix,
- const RealType& matrixMultiplicator )
-{
-
-}*/
-
template< typename Real,
typename Device,
typename Index,
diff --git a/src/UnitTests/Matrices/SparseMatrixTest.h b/src/UnitTests/Matrices/SparseMatrixTest.h
index 2281d229709f53214d14c420dd92ff67d093801c..c4086e634807b842b8aa67d23e7fac303854a958 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest.h
+++ b/src/UnitTests/Matrices/SparseMatrixTest.h
@@ -99,10 +99,17 @@ TYPED_TEST( MatrixTest, reduceRows )
test_reduceRows< MatrixType >();
}
-TYPED_TEST( MatrixTest, saveAndLoadTest )
+TYPED_TEST( MatrixTest, saveAndLoad )
{
using MatrixType = typename TestFixture::MatrixType;
test_SaveAndLoad< MatrixType >( saveAndLoadFileName );
}
+
+TYPED_TEST( MatrixTest, getTransposition )
+{
+ using MatrixType = typename TestFixture::MatrixType;
+
+ test_getTransposition< MatrixType >();
+}
#endif
diff --git a/src/UnitTests/Matrices/SparseMatrixTest.hpp b/src/UnitTests/Matrices/SparseMatrixTest.hpp
index 0e7884859d63a07c6fd97cd41354fbda590dde79..6a74e06be6eb7eb5dae680ed8db2a99f7cc6840b 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest.hpp
+++ b/src/UnitTests/Matrices/SparseMatrixTest.hpp
@@ -1356,4 +1356,68 @@ void test_SaveAndLoad( const char* filename )
EXPECT_EQ( std::remove( filename ), 0 );
}
+template< typename Matrix >
+void test_getTransposition()
+{
+ using RealType = typename Matrix::RealType;
+ using IndexType = typename Matrix::IndexType;
+
+ /*
+ * Sets up the following 5x4 sparse matrix:
+ *
+ * / 1 2 3 0 \
+ * | 0 4 0 5 |
+ * | 6 7 8 0 |
+ * | 0 9 10 11 |
+ * \ 12 0 0 0 /
+ */
+
+ const IndexType m_rows = 5;
+ const IndexType m_cols = 4;
+
+ Matrix matrix( m_rows, m_cols );
+ typename Matrix::RowsCapacitiesType capacities( m_rows, 4 );
+ matrix.setRowCapacities( capacities );
+
+ RealType value = 1;
+ for( IndexType i = 0; i < m_cols - 1; i++ ) // 0th row
+ matrix.setElement( 0, i, value++ );
+
+ matrix.setElement( 1, 1, value++ ); // 1st row
+ matrix.setElement( 1, 3, value++ );
+
+ for( IndexType i = 0; i < m_cols - 1; i++ ) // 2nd row
+ matrix.setElement( 2, i, value++ );
+
+ for( IndexType i = 1; i < m_cols; i++ ) // 3rd row
+ matrix.setElement( 3, i, value++ );
+
+ matrix.setElement( 4, 0, value++ ); // 4th row
+
+ Matrix transposed;
+ transposed.getTransposition( matrix );
+
+ EXPECT_EQ( transposed.getElement( 0, 0 ), 1 );
+ EXPECT_EQ( transposed.getElement( 1, 0 ), 2 );
+ EXPECT_EQ( transposed.getElement( 2, 0 ), 3 );
+ EXPECT_EQ( transposed.getElement( 3, 0 ), 0 );
+
+ EXPECT_EQ( transposed.getElement( 0, 1 ), 0 );
+ EXPECT_EQ( transposed.getElement( 1, 1 ), 4 );
+ EXPECT_EQ( transposed.getElement( 2, 1 ), 0 );
+ EXPECT_EQ( transposed.getElement( 3, 1 ), 5 );
+
+ EXPECT_EQ( transposed.getElement( 0, 2 ), 6 );
+ EXPECT_EQ( transposed.getElement( 1, 2 ), 7 );
+ EXPECT_EQ( transposed.getElement( 2, 2 ), 8 );
+ EXPECT_EQ( transposed.getElement( 3, 2 ), 0 );
+
+ EXPECT_EQ( transposed.getElement( 0, 3 ), 0 );
+ EXPECT_EQ( transposed.getElement( 1, 3 ), 9 );
+ EXPECT_EQ( transposed.getElement( 2, 3 ), 10 );
+ EXPECT_EQ( transposed.getElement( 3, 3 ), 11 );
+
+ EXPECT_EQ( transposed.getElement( 0, 4 ), 12 );
+}
+
#endif