From 733401dcfd44ac508f43d1f5731bfc795226ddf0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Fri, 17 Jan 2020 22:51:01 +0100
Subject: [PATCH] Fixed tridiagonal to sparse matrix assignment. Debugging
multidiagonal to sparse matrix assignment.
---
.../Matrices/MultidiagonalMatrixRowView.hpp | 2 +-
src/TNL/Matrices/MultidiagonalMatrixView.hpp | 3 +-
src/TNL/Matrices/SparseMatrix.hpp | 3 +-
src/UnitTests/Matrices/SparseMatrixCopyTest.h | 117 ++++++++++++++----
4 files changed, 98 insertions(+), 27 deletions(-)
diff --git a/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp b/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp
index 88aae3f150..855b8463aa 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp
+++ b/src/TNL/Matrices/MultidiagonalMatrixRowView.hpp
@@ -30,7 +30,7 @@ auto
MultidiagonalMatrixRowView< ValuesView, Indexer, DiagonalsShiftsView >::
getSize() const -> IndexType
{
- return indexer.getRowSize();
+ return diagonalsShifts.getSize();//indexer.getRowSize( rowIdx );
}
template< typename ValuesView, typename Indexer, typename DiagonalsShiftsView >
diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.hpp b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
index 2839c997aa..96312d03cf 100644
--- a/src/TNL/Matrices/MultidiagonalMatrixView.hpp
+++ b/src/TNL/Matrices/MultidiagonalMatrixView.hpp
@@ -393,12 +393,13 @@ forRows( IndexType first, IndexType last, Function& function ) const
const IndexType diagonalsCount = this->diagonalsShifts.getSize();
const IndexType columns = this->getColumns();
const auto indexer = this->indexer;
+ bool compute( true );
auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
for( IndexType localIdx = 0; localIdx < diagonalsCount; localIdx++ )
{
const IndexType columnIdx = rowIdx + diagonalsShifts_view[ localIdx ];
if( columnIdx >= 0 && columnIdx < columns )
- function( rowIdx, localIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx ) ] );
+ function( rowIdx, localIdx, columnIdx, values_view[ indexer.getGlobalIndex( rowIdx, localIdx ) ], compute );
}
};
Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index 6d3d9d8b8b..66813a732c 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -682,7 +682,6 @@ operator=( const RHSMatrix& matrix )
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();
@@ -697,7 +696,7 @@ operator=( const RHSMatrix& matrix )
auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx_, RHSIndexType columnIndex, const RHSRealType& value, bool& compute ) mutable {
RealType inValue( 0.0 );
IndexType localIdx( rowLocalIndexes_view[ rowIdx ] );
- if( columnIndex != paddingIndex )
+ if( value != 0.0 && columnIndex != paddingIndex )
{
IndexType thisGlobalIdx = segments_view.getGlobalIndex( rowIdx, localIdx++ );
columns_view[ thisGlobalIdx ] = columnIndex;
diff --git a/src/UnitTests/Matrices/SparseMatrixCopyTest.h b/src/UnitTests/Matrices/SparseMatrixCopyTest.h
index b285f5e058..7f29191b6f 100644
--- a/src/UnitTests/Matrices/SparseMatrixCopyTest.h
+++ b/src/UnitTests/Matrices/SparseMatrixCopyTest.h
@@ -421,7 +421,6 @@ void testConversion()
checkAntiTriDiagMatrix( antiTriDiag1 );
Matrix2 antiTriDiag2;
- //TNL::Matrices::copySparseMatrix( antiTriDiag2, antiTriDiag1 );
antiTriDiag2 = antiTriDiag1;
checkAntiTriDiagMatrix( antiTriDiag2 );
}
@@ -433,7 +432,6 @@ void testConversion()
checkUnevenRowSizeMatrix( unevenRowSize1 );
Matrix2 unevenRowSize2;
- //TNL::Matrices::copySparseMatrix( unevenRowSize2, unevenRowSize1 );
unevenRowSize2 = unevenRowSize1;
checkUnevenRowSizeMatrix( unevenRowSize2 );
}
@@ -451,21 +449,19 @@ void tridiagonalMatrixAssignment()
const IndexType rows( 10 ), columns( 10 );
TridiagonalHost hostMatrix( rows, columns );
- for( IndexType i = 0; i < columns; i++ )
- for( IndexType j = TNL::max( 0, i - 1 ); j < TNL::min( rows, i + 2 ); j++ )
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = TNL::max( 0, i - 1 ); j < TNL::min( columns, i + 2 ); j++ )
hostMatrix.setElement( i, j, i + j );
- std::cerr << hostMatrix << std::endl;
Matrix matrix;
matrix = hostMatrix;
- std::cerr << matrix << std::endl;
using RowCapacitiesType = typename Matrix::RowsCapacitiesType;
RowCapacitiesType rowCapacities;
matrix.getCompressedRowLengths( rowCapacities );
- RowCapacitiesType exactRowLengths{ 0, 3, 3, 3, 3, 3, 3, 3, 3, 2 };
+ RowCapacitiesType exactRowLengths{ 1, 3, 3, 3, 3, 3, 3, 3, 3, 2 };
EXPECT_EQ( rowCapacities, exactRowLengths );
- for( IndexType i = 0; i < columns; i++ )
- for( IndexType j = 0; j < rows; j++ )
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < columns; j++ )
{
if( abs( i - j ) > 1 )
EXPECT_EQ( matrix.getElement( i, j ), 0.0 );
@@ -476,15 +472,11 @@ void tridiagonalMatrixAssignment()
#ifdef HAVE_CUDA
TridiagonalCuda cudaMatrix( rows, columns );
cudaMatrix = hostMatrix;
- /*for( IndexType i = 0; i < columns; i++ )
- for( IndexType j = TNL::max( 0, i - 1 ); j < TNL::min( row, i + 1 ); j++ )
- cudaMatrix.setElement( i, j, i + j );*/
-
matrix = cudaMatrix;
matrix.getCompressedRowLengths( rowCapacities );
EXPECT_EQ( rowCapacities, exactRowLengths );
- for( IndexType i = 0; i < columns; i++ )
- for( IndexType j = 0; j < rows; j++ )
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < columns; j++ )
{
if( abs( i - j ) > 1 )
EXPECT_EQ( matrix.getElement( i, j ), 0.0 );
@@ -492,7 +484,58 @@ void tridiagonalMatrixAssignment()
EXPECT_EQ( matrix.getElement( i, j ), i + j );
}
#endif
+}
+template< typename Matrix >
+void multidiagonalMatrixAssignment()
+{
+ using RealType = typename Matrix::RealType;
+ using DeviceType = typename Matrix::DeviceType;
+ using IndexType = typename Matrix::IndexType;
+
+ using MultidiagonalHost = TNL::Matrices::Multidiagonal< RealType, TNL::Devices::Host, IndexType >;
+ using MultidiagonalCuda = TNL::Matrices::Multidiagonal< RealType, TNL::Devices::Cuda, IndexType >;
+ using DiagonalsShiftsType = typename MultidiagonalHost::DiagonalsShiftsType;
+ DiagonalsShiftsType diagonals{ -4, -2, 0, 1, 3, 5 };
+
+ const IndexType rows( 10 ), columns( 10 );
+ MultidiagonalHost hostMatrix( rows, columns, diagonals );
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < columns; j++ )
+ if( diagonals.containsValue( i - j ) )
+ hostMatrix.setElement( i, j, i + j );
+
+ Matrix matrix;
+ matrix = hostMatrix;
+ using RowCapacitiesType = typename Matrix::RowsCapacitiesType;
+ RowCapacitiesType rowCapacities;
+ matrix.getCompressedRowLengths( rowCapacities );
+ RowCapacitiesType exactRowLengths{ 1, 3, 3, 3, 3, 3, 3, 3, 3, 2 };
+ EXPECT_EQ( rowCapacities, exactRowLengths );
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < columns; j++ )
+ {
+ if( diagonals.containsValue( i - j ) )
+ EXPECT_EQ( matrix.getElement( i, j ), 0.0 );
+ else
+ EXPECT_EQ( matrix.getElement( i, j ), i + j );
+ }
+
+#ifdef HAVE_CUDA
+ MultidiagonalCuda cudaMatrix( rows, columns, diagonals );
+ cudaMatrix = hostMatrix;
+ matrix = cudaMatrix;
+ matrix.getCompressedRowLengths( rowCapacities );
+ EXPECT_EQ( rowCapacities, exactRowLengths );
+ for( IndexType i = 0; i < rows; i++ )
+ for( IndexType j = 0; j < columns; j++ )
+ {
+ if( diagonals.containsValue( i - j ) > 1 )
+ EXPECT_EQ( matrix.getElement( i, j ), 0.0 );
+ else
+ EXPECT_EQ( matrix.getElement( i, j ), i + j );
+ }
+#endif
}
template< typename Matrix >
@@ -530,10 +573,6 @@ void denseMatrixAssignment()
#ifdef HAVE_CUDA
DenseCuda cudaMatrix( rows, columns );
cudaMatrix = hostMatrix;
- /*for( IndexType i = 0; i < columns; i++ )
- for( IndexType j = 0; j <= i; j++ )
- cudaMatrix.setElement( i, j, i + j );*/
-
matrix = cudaMatrix;
matrix.getCompressedRowLengths( rowCapacities );
EXPECT_EQ( rowCapacities, exactRowLengths );
@@ -547,7 +586,7 @@ void denseMatrixAssignment()
}
#endif
}
-/*
+
TEST( SparseMatrixCopyTest, CSR_HostToHost )
{
testCopyAssignment< CSR_host, CSR_host >();
@@ -616,8 +655,8 @@ TEST( SparseMatrixCopyTest, SlicedEllpack_CudaToCuda )
}
#endif
-
-// test conversion between formats
+////
+// Test of conversion between formats
TEST( SparseMatrixCopyTest, CSR_to_Ellpack_host )
{
testConversion< CSR_host, E_host >();
@@ -679,7 +718,6 @@ TEST( SparseMatrixCopyTest, SlicedEllpack_to_Ellpack_cuda )
testConversion< SE_cuda, E_cuda >();
}
#endif
-*/
////
// Tridiagonal matrix assignment test
@@ -715,8 +753,41 @@ TEST( SparseMatrixCopyTest, TridiagonalMatrixAssignment_to_SlicedEllpack_cuda )
}
#endif // HAVE_CUDA
+////
+// Multidiagonal matrix assignment test
+TEST( SparseMatrixCopyTest, MultidiagonalMatrixAssignment_to_CSR_host )
+{
+ multidiagonalMatrixAssignment< CSR_host >();
+}
+
+TEST( SparseMatrixCopyTest, MultidiagonalMatrixAssignment_to_Ellpack_host )
+{
+ multidiagonalMatrixAssignment< E_host >();
+}
+
+TEST( SparseMatrixCopyTest, MultidiagonalMatrixAssignment_to_SlicedEllpack_host )
+{
+ multidiagonalMatrixAssignment< SE_host >();
+}
+
+#ifdef HAVE_CUDA
+TEST( SparseMatrixCopyTest, MultidiagonalMatrixAssignment_to_CSR_cuda )
+{
+ multidiagonalMatrixAssignment< CSR_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, MultidiagonalMatrixAssignment_to_Ellpack_cuda )
+{
+ multidiagonalMatrixAssignment< E_cuda >();
+}
+TEST( SparseMatrixCopyTest, MultidiagonalMatrixAssignment_to_SlicedEllpack_cuda )
+{
+ multidiagonalMatrixAssignment< SE_cuda >();
+}
+#endif // HAVE_CUDA
+////
// Dense matrix assignment test
TEST( SparseMatrixCopyTest, DenseMatrixAssignment_to_CSR_host )
{
--
GitLab