Conversion to a fixed sparsity pattern. (#147)

Source/C/CMakeLists.txt

@@ -8,6 +8,7 @@ set(Csrc
     InverseSolvers_c.h
     LinearSolvers_c.h
     LoadBalancer_c.h
+    MatrixConversion_c.h
     MatrixMemoryPool_c.h
     PMatrixMemoryPool_c.h
     PSMatrix_c.h

Source/C/MatrixConversion_c.h

@@ -0,0 +1,6 @@
+#ifndef MATRIXCONVERSIONMODULE_ch
+#define MATRIXCONVERSIONMODULE_ch
+
+void SnapMatrixToSparsityPattern_wrp(int *ih_matA, const int *ih_matB);
+
+#endif

Source/C/PSMatrix_c.h

@@ -30,6 +30,7 @@ void FillMatrixPermutation_ps_wrp(int *ih_this, const int *ih_permutation,
 void FillMatrixIdentity_ps_wrp(int *ih_this);
 void GetMatrixActualDimension_ps_wrp(const int *ih_this, int *size);
 void GetMatrixLogicalDimension_ps_wrp(const int *ih_this, int *size);
+void GetMatrixSize_ps_wrp(const int *ih_this, int *size);
 void GetMatrixTripletList_psr_wrp(const int *ih_this, int *ih_triplet_list);
 void GetMatrixTripletList_psc_wrp(const int *ih_this, int *ih_triplet_list);
 void GetMatrixBlock_psr_wrp(const int *ih_this, int *ih_triplet_list,

Source/CPlusPlus/CMakeLists.txt

@@ -9,6 +9,7 @@ set(Csrc
     InverseSolvers.cc
     LinearSolvers.cc
     LoadBalancer.cc
+    MatrixConversion.cc
     MatrixMapper.cc
     MatrixMemoryPool.cc
     Permutation.cc
@@ -36,6 +37,7 @@ set(Chead
     InverseSolvers.h
     LinearSolvers.h
     LoadBalancer.h
+    MatrixConversion.h
     MatrixMapper.h
     MatrixMemoryPool.h
     Permutation.h

Source/CPlusPlus/MatrixConversion.cc

@@ -0,0 +1,17 @@
+#include "MatrixConversion.h"
+#include "PSMatrix.h"
+using namespace NTPoly;
+
+////////////////////////////////////////////////////////////////////////////////
+extern "C" {
+#include "MatrixConversion_c.h"
+}
+
+////////////////////////////////////////////////////////////////////////////////
+namespace NTPoly {
+//////////////////////////////////////////////////////////////////////////////
+void MatrixConversion::SnapMatrixToSparsityPattern(Matrix_ps& mata,
+	                                               const Matrix_ps& matb) {
+  SnapMatrixToSparsityPattern_wrp(mata.ih_this, matb.ih_this);
+}
+} // namespace NTPoly

Source/CPlusPlus/MatrixConversion.h

@@ -0,0 +1,22 @@
+#ifndef MATRIXCONVERSION_h
+#define MATRIXCONVERSION_h
+
+#include "Wrapper.h"
+
+////////////////////////////////////////////////////////////////////////////////
+namespace NTPoly {
+class Matrix_ps;
+class MatrixConversion {
+public:
+  //! Some codes use a fixed sparsity pattern for a matrix instead of filtering
+  //! small values. Using this routine, the matrix is filled to have the same
+  //! pattern as the second matrix argument. Zeros of the sparsity pattern are
+  //! left in, whereas values outside the sparsity are removed. This can 
+  //! faciliate conversion between formats.
+  //! \param mata the matrix to modify.
+  //! \param matb the matrix which defines the sparsity pattern.
+  static void SnapMatrixToSparsityPattern(Matrix_ps& mata, 
+                                          const Matrix_ps& matb);
+};
+} // namespace NTPoly
+#endif

Source/CPlusPlus/PSMatrix.cc

@@ -105,6 +105,13 @@ int Matrix_ps::GetActualDimension() const {
   return temp;
 }
 
+//////////////////////////////////////////////////////////////////////////////
+int Matrix_ps::GetSize() const {
+  int temp;
+  GetMatrixSize_ps_wrp(ih_this, &temp);
+  return temp;
+}
+
 //////////////////////////////////////////////////////////////////////////////
 void Matrix_ps::GetTripletList(TripletList_r &triplet_list) const {
   GetMatrixTripletList_psr_wrp(ih_this, triplet_list.ih_this);

Source/CPlusPlus/PSMatrix.h

@@ -14,6 +14,7 @@ class TripletList_r;
 class TripletList_c;
 class ProcessGrid;
 class MatrixMapper;
+class MatrixConversion;
 ////////////////////////////////////////////////////////////////////////////////
 //! A Module For Performing Distributed Sparse Matrix Operations.
 class Matrix_ps {
@@ -67,6 +68,8 @@ class Matrix_ps {
   int GetActualDimension() const;
   //! the logical dimension is scaled so each process has an even slice.
   int GetLogicalDimension() const;
+  //! Get the total number of non-zero entries in the matrix
+  int GetSize() const;
   //! Extracts a triplet list of the data that is stored on this process.
   //! Data is returned with absolute coordinates.
   //! \param triplet_list the list to fill.
@@ -174,6 +177,7 @@ class Matrix_ps {
   template <class T> friend class TripletList;
   friend class PMatrixMemoryPool;
   friend class MatrixMapper;
+  friend class MatrixConversion;
 };
 } // namespace NTPoly
 #endif

Source/Fortran/CMakeLists.txt

@@ -14,6 +14,7 @@ set(Fsrc
     LinearSolversModule.F90
     LoadBalancerModule.F90
     LoggingModule.F90
+    MatrixConversionModule.F90
     MatrixMarketModule.F90
     MatrixMapsModule.F90
     MatrixMemoryPoolModule.F90

Source/Fortran/MatrixConversionModule.F90

@@ -0,0 +1,72 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!> This module contains helper routines for converting an NTPoly matrix
+!> to data structures used in other programs.
+MODULE MatrixConversionModule
+  USE DataTypesModule, ONLY : NTREAL
+  USE MatrixMapsModule, ONLY : MapMatrix_psr
+  USE PSMatrixModule, ONLY : Matrix_ps, ConvertMatrixToReal, CopyMatrix, &
+       & DestructMatrix
+  USE PSMatrixAlgebraModule, ONLY : PairwiseMultiplyMatrix, ScaleMatrix, &
+       & IncrementMatrix
+  IMPLICIT NONE
+  PRIVATE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  PUBLIC :: SnapMatrixToSparsityPattern
+CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Some codes use a fixed sparsity pattern for a matrix instead of filtering
+  !> small values. Using this routine, the matrix is filled to have the same
+  !> pattern as the second matrix argument. Zeros of the sparsity pattern are
+  !> left in, whereas values outside the sparsity are removed. This can 
+  !> faciliate conversion between formats.
+  SUBROUTINE SnapMatrixToSparsityPattern(mat, pattern)
+    !> The matrix to modify.
+    TYPE(Matrix_ps), INTENT(INOUT) :: mat
+    !> The matrix which defines the sparsity pattern.
+    TYPE(Matrix_ps), INTENT(IN) :: pattern
+    !! Local Variables
+    TYPE(Matrix_ps) :: filtered
+    TYPE(Matrix_ps) :: pattern_1s
+    TYPE(Matrix_ps) :: pattern_0s
+    TYPE(Matrix_ps) :: pattern_real
+    INTEGER :: II
+
+    !! First we need to make sure that the sparsity pattern is all 1s.
+    IF (pattern%is_complex) THEN
+       CALL ConvertMatrixToReal(pattern, pattern_real)
+    ELSE
+       CALL CopyMatrix(pattern, pattern_real)
+    END IF
+    CALL MapMatrix_psr(pattern_real, pattern_1s, SetMatrixToOne)
+
+    !! Then all zeros
+    CALL CopyMatrix(pattern_1s, pattern_0s)
+    CALL ScaleMatrix(pattern_0s, 0.0_NTREAL)
+
+    !! Here we add in the zero values that were missing from the original
+    !! matrix. The secret here is that if you use a negative threshold, we
+    !! never filter a value.
+    CALL IncrementMatrix(pattern_0s, mat, threshold_in=-1.0_NTREAL)
+
+    !! Next, we zero out values outside of the sparsity pattern.
+    CALL CopyMatrix(mat, filtered)
+    CALL PairwiseMultiplyMatrix(pattern_1s, filtered, mat)
+
+    !! Cleanup
+    CALL DestructMatrix(pattern_1s)
+    CALL DestructMatrix(pattern_0s)
+    CALL DestructMatrix(pattern_real)
+    CALL DestructMatrix(filtered)
+
+  END SUBROUTINE SnapMatrixToSparsityPattern
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  FUNCTION SetMatrixToOne(row, column, val) RESULT(valid)
+    INTEGER, INTENT(INOUT), OPTIONAL :: row
+    INTEGER, INTENT(INOUT), OPTIONAL :: column
+    REAL(NTREAL), INTENT(INOUT), OPTIONAL :: val
+    LOGICAL :: valid
+
+    val = 1.0_NTREAL
+    valid = .TRUE.
+  END FUNCTION SetMatrixToOne
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+END MODULE

Source/Swig/NTPolySwig.i

@@ -18,6 +18,7 @@
 #include "InverseSolvers.h"
 #include "LinearSolvers.h"
 #include "LoadBalancer.h"
+#include "MatrixConversion.h"
 #include "MatrixMapper.h"
 #include "MatrixMemoryPool.h"
 #include "Permutation.h"
@@ -53,6 +54,7 @@ using namespace NTPoly;
 %include "InverseSolvers.h"
 %include "LinearSolvers.h"
 %include "LoadBalancer.h"
+%include "MatrixConversion.h"
 %include "MatrixMapper.h"
 %include "MatrixMemoryPool.h"
 %include "Permutation.h"

Source/Wrapper/CMakeLists.txt

@@ -9,6 +9,7 @@ set(Wsrc
     InverseSolversModule_wrp.F90
     LinearSolversModule_wrp.F90
     LoadBalancerModule_wrp.F90
+    MatrixConversionModule_wrp.F90
     MatrixMemoryPoolModule_wrp.F90
     PermutationModule_wrp.F90
     PMatrixMemoryPoolModule_wrp.F90

Source/Wrapper/MatrixConversionModule_wrp.F90

@@ -0,0 +1,31 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!> This module contains helper routines for converting an NTPoly matrix
+!> to data structures used in other programs.
+MODULE MatrixConversionModule_wrp
+  USE MatrixConversionModule
+  USE PSMatrixModule_wrp, ONLY : Matrix_ps_wrp
+  USE WrapperModule, ONLY : SIZE_wrp
+  USE ISO_C_BINDING, ONLY : C_INT
+  IMPLICIT NONE
+  PRIVATE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  PUBLIC :: SnapMatrixToSparsityPattern_wrp
+CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Some codes use a fixed sparsity pattern for a matrix instead of filtering
+  !> small values. Using this routine, the matrix is filled to have the same
+  !> pattern as the second matrix argument. Zeros of the sparsity pattern are
+  !> left in, whereas values outside the sparsity are removed. This can 
+  !> faciliate conversion between formats.
+  SUBROUTINE SnapMatrixToSparsityPattern_wrp(ih_matA, ih_matB) &
+       & BIND(C,NAME="SnapMatrixToSparsityPattern_wrp")
+    INTEGER(KIND=C_INT), INTENT(INOUT) :: ih_matA(SIZE_wrp)
+    INTEGER(KIND=C_INT), INTENT(IN) :: ih_matB(SIZE_wrp)
+    TYPE(Matrix_ps_wrp) :: h_matA
+    TYPE(Matrix_ps_wrp) :: h_matB
+
+    h_matA = TRANSFER(ih_matA,h_matA)
+    h_matB = TRANSFER(ih_matB,h_matB)
+    CALL SnapMatrixToSparsityPattern(h_matA%data, h_matB%data)
+  END SUBROUTINE SnapMatrixToSparsityPattern_wrp
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+END MODULE

Source/Wrapper/PSMatrixModule_wrp.F90

@@ -29,6 +29,7 @@ MODULE PSMatrixModule_wrp
   PUBLIC :: FillMatrixIdentity_ps_wrp
   PUBLIC :: GetMatrixActualDimension_ps_wrp
   PUBLIC :: GetMatrixLogicalDimension_ps_wrp
+  PUBLIC :: GetMatrixSize_ps_wrp
   PUBLIC :: GetMatrixTripletList_psr_wrp
   PUBLIC :: GetMatrixTripletList_psc_wrp
   PUBLIC :: GetMatrixBlock_psr_wrp
@@ -289,6 +290,17 @@ SUBROUTINE GetMatrixLogicalDimension_ps_wrp(ih_this, mat_dimension) &
     h_this = TRANSFER(ih_this,h_this)
     mat_dimension = GetMatrixLogicalDimension(h_this%data)
   END SUBROUTINE GetMatrixLogicalDimension_ps_wrp
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Get the total number of non-zero entries in the distributed sparse matrix.
+  SUBROUTINE GetMatrixSize_ps_wrp(ih_this, matsize) &
+       & BIND(c,NAME="GetMatrixSize_ps_wrp")
+    INTEGER(KIND=c_int), INTENT(IN) :: ih_this(SIZE_wrp)
+    INTEGER(KIND=c_int), INTENT(OUT) :: matsize
+    TYPE(Matrix_ps_wrp) :: h_this
+
+    h_this = TRANSFER(ih_this,h_this)
+    matsize = GetMatrixSize(h_this%data)
+  END SUBROUTINE GetMatrixSize_ps_wrp
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> Extracts a triplet list of the data that is stored on this process.
   SUBROUTINE GetMatrixTripletList_psr_wrp(ih_this, ih_triplet_list) &

UnitTests/test_psmatrix.py

@@ -427,6 +427,43 @@ def __call__(self):
 
             self.check_result()
 
+    def test_snap(self):
+        '''Test the sparsity pattern setting routine.'''
+        from scipy.sparse import dok_matrix, csr_matrix
+
+        for param in self.parameters:
+            matrix1 = param.create_matrix(self.complex)
+            matrix2 = param.create_matrix(self.complex)
+            self.write_matrix(matrix1, self.input_file1)
+            self.write_matrix(matrix2, self.input_file2)
+
+            # Convert to DOK format
+            matrix1 = dok_matrix(matrix1)
+            matrix2 = dok_matrix(matrix2)
+            self.CheckMat = dok_matrix(matrix1*0)
+
+            # Iterate over keys of mat2.
+            for idx in matrix2.keys():
+                if idx in matrix1:
+                    self.CheckMat[idx] = matrix1[idx]
+                else:
+                    self.CheckMat[idx] = 0
+
+            ntmatrix1 = nt.Matrix_ps(self.input_file1, False)
+            ntmatrix2 = nt.Matrix_ps(self.input_file2, False)
+
+            nt.MatrixConversion.SnapMatrixToSparsityPattern(ntmatrix1,
+                                                            ntmatrix2)
+            ntmatrix1.WriteToMatrixMarket(self.result_file)
+            comm.barrier()
+
+            self.check_result()
+
+            # We have to check that we have the same number of values.
+            nnz = ntmatrix1.GetSize()
+            if self.my_rank == 0:
+                self.assertEqual(len(matrix2.keys()), nnz)
+
 
 class TestPSMatrix_c(TestPSMatrix):
     '''Specialization for complex matrices'''