From 15f4dfeef13acfacc29352a3b7e549bbaab71771 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jakub=20Klinkovsk=C3=BD?= <klinkjak@fjfi.cvut.cz>
Date: Thu, 4 Oct 2018 23:24:07 +0200
Subject: [PATCH] Added copyAdjacencyStructure function
---
src/TNL/Matrices/SparseOperations.h | 9 ++++
src/TNL/Matrices/SparseOperations_impl.h | 69 ++++++++++++++++++++++++
2 files changed, 78 insertions(+)
diff --git a/src/TNL/Matrices/SparseOperations.h b/src/TNL/Matrices/SparseOperations.h
index e375d596ad..6da801ee88 100644
--- a/src/TNL/Matrices/SparseOperations.h
+++ b/src/TNL/Matrices/SparseOperations.h
@@ -22,6 +22,15 @@ namespace Matrices {
template< typename Matrix1, typename Matrix2 >
void copySparseMatrix( Matrix1& A, const Matrix2& B );
+// NOTE: if `has_symmetric_pattern`, the sparsity pattern of `A` is assumed
+// to be symmetric and it is just copied to `B`. Otherwise, the sparsity
+// pattern of `A^T + A` is copied to `B`.
+template< typename Matrix, typename AdjacencyMatrix >
+void
+copyAdjacencyStructure( const Matrix& A, AdjacencyMatrix& B,
+ bool has_symmetric_pattern = false,
+ bool ignore_diagonal = true );
+
} // namespace Matrices
} // namespace TNL
diff --git a/src/TNL/Matrices/SparseOperations_impl.h b/src/TNL/Matrices/SparseOperations_impl.h
index 6f9adf9314..d083fb21a4 100644
--- a/src/TNL/Matrices/SparseOperations_impl.h
+++ b/src/TNL/Matrices/SparseOperations_impl.h
@@ -13,6 +13,7 @@
#pragma once
#include <type_traits>
+#include <stdexcept>
#include <TNL/Pointers/DevicePointer.h>
@@ -178,5 +179,73 @@ copySparseMatrix( Matrix1& A, const Matrix2& B )
copySparseMatrix_impl( A, B );
}
+
+template< typename Matrix, typename AdjacencyMatrix >
+void
+copyAdjacencyStructure( const Matrix& A, AdjacencyMatrix& B,
+ bool has_symmetric_pattern,
+ bool ignore_diagonal )
+{
+ static_assert( std::is_same< typename Matrix::DeviceType, Devices::Host >::value,
+ "The function is not implemented for CUDA matrices - it would require atomic insertions "
+ "of elements into the sparse format." );
+ static_assert( std::is_same< typename Matrix::DeviceType, typename AdjacencyMatrix::DeviceType >::value,
+ "The matrices must be allocated on the same device." );
+ static_assert( std::is_same< typename Matrix::IndexType, typename AdjacencyMatrix::IndexType >::value,
+ "The matrices must have the same IndexType." );
+// static_assert( std::is_same< typename AdjacencyMatrix::RealType, bool >::value,
+// "The RealType of the adjacency matrix must be bool." );
+
+ using IndexType = typename Matrix::IndexType;
+
+ if( A.getRows() != A.getColumns() ) {
+ throw std::logic_error( "The matrix is not square: " + std::to_string( A.getRows() ) + " rows, "
+ + std::to_string( A.getColumns() ) + " columns." );
+ }
+
+ const IndexType N = A.getRows();
+ B.setDimensions( N, N );
+
+ // set row lengths
+ typename AdjacencyMatrix::CompressedRowLengthsVector rowLengths;
+ rowLengths.setSize( N );
+ rowLengths.setValue( 0 );
+ for( IndexType i = 0; i < A.getRows(); i++ ) {
+ const int maxLength = A.getRowLength( i );
+ const auto row = A.getRow( i );
+ IndexType length = 0;
+ for( int c_j = 0; c_j < maxLength; c_j++ ) {
+ const IndexType j = row.getElementColumn( c_j );
+ if( j >= A.getColumns() )
+ break;
+ length++;
+ if( ! has_symmetric_pattern && i != j )
+ if( A.getElement( j, i ) == 0 )
+ rowLengths[ j ]++;
+ }
+ if( ignore_diagonal )
+ length--;
+ rowLengths[ i ] += length;
+ }
+ B.setCompressedRowLengths( rowLengths );
+
+ // set non-zeros
+ for( IndexType i = 0; i < A.getRows(); i++ ) {
+ const int maxLength = A.getRowLength( i );
+ const auto row = A.getRow( i );
+ for( int c_j = 0; c_j < maxLength; c_j++ ) {
+ const IndexType j = row.getElementColumn( c_j );
+ if( j >= A.getColumns() )
+ break;
+ if( ! ignore_diagonal || i != j )
+ if( A.getElement( i, j ) != 0 ) {
+ B.setElement( i, j, true );
+ if( ! has_symmetric_pattern )
+ B.setElement( j, i, true );
+ }
+ }
+ }
+}
+
} // namespace Matrices
} // namespace TNL
--
GitLab