Symmetry helpers

Source/C/PSMatrix_c.h

@@ -60,9 +60,10 @@ void MatrixMultiply_ps_wrp(const int *ih_matA, const int *ih_matB, int *ih_matC,
                            const double *threshold_in, int *ih_memory_pool_in);
 void ScaleMatrix_ps_wrp(int *ih_this, const double *constant);
 double MatrixNorm_ps_wrp(const int *ih_this);
+double MeasureAsymmetry_ps_wrp(const int *ih_this);
 void MatrixTrace_ps_wrp(const int *ih_this, double *trace_val);
 int IsIdentity_ps_wrp(const int *ih_this);
+void SymmetrizeMatrix_ps_wrp(int *ih_this);
 void MatrixDiagonalScale_psr_wrp(int *ih_mat, const int *ih_tlist);
 void MatrixDiagonalScale_psc_wrp(int *ih_mat, const int *ih_tlist);
-
 #endif

Source/CPlusPlus/PSMatrix.cc

@@ -237,6 +237,14 @@ double Matrix_ps::Trace() const {
   return temp_val;
 }
 
+//////////////////////////////////////////////////////////////////////////////
+double Matrix_ps::MeasureAsymmetry() const {
+  return MeasureAsymmetry_ps_wrp(ih_this);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+void Matrix_ps::Symmetrize() { SymmetrizeMatrix_ps_wrp(ih_this); }
+
 //////////////////////////////////////////////////////////////////////////////
 Matrix_ps::~Matrix_ps() { DestructMatrix_ps_wrp(ih_this); }
 } // namespace NTPoly

Source/CPlusPlus/PSMatrix.h

@@ -165,8 +165,12 @@ class Matrix_ps {
   void Scale(double constant);
   //! compute the norm of a matrix.
   double Norm() const;
+  //! compute the asymmetry (norm(A - A.T)) of a matrix.
+  double MeasureAsymmetry() const;
   //! compute the trace of a matrix.
   double Trace() const;
+  //! symmetrize the matrix
+  void Symmetrize();
   //!\param tlist the triplet list.
   //!\param threshold for flushing small values.
   void DiagonalScale(const NTPoly::TripletList_r &tlist);

Source/Fortran/PSMatrixAlgebraModule.F90

@@ -13,7 +13,8 @@ MODULE PSMatrixAlgebraModule
        & ConstructMatrixMemoryPool
   USE PSMatrixModule, ONLY : Matrix_ps, ConstructEmptyMatrix, CopyMatrix, &
        & DestructMatrix, ConvertMatrixToComplex, ConjugateMatrix, &
-       & MergeMatrixLocalBlocks, IsIdentity, SplitMatrixToLocalBlocks
+       & MergeMatrixLocalBlocks, IsIdentity, TransposeMatrix, &
+       & SplitMatrixToLocalBlocks
   USE SMatrixAlgebraModule, ONLY : MatrixMultiply, MatrixGrandSum, &
        & PairwiseMultiplyMatrix, IncrementMatrix, ScaleMatrix, &
        & MatrixColumnNorm, MatrixDiagonalScale
@@ -37,6 +38,8 @@ MODULE PSMatrixAlgebraModule
   PUBLIC :: ScaleMatrix
   PUBLIC :: MatrixTrace
   PUBLIC :: SimilarityTransform
+  PUBLIC :: MeasureAsymmetry
+  PUBLIC :: SymmetrizeMatrix
   PUBLIC :: MatrixDiagonalScale
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   INTERFACE MatrixSigma
@@ -561,6 +564,38 @@ SUBROUTINE MatrixTrace_psr(this, trace_value)
 #undef TLIST
     END IF
   END SUBROUTINE MatrixTrace_psr
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Measure the asymmetry of a matrix NORM(A - A.T)
+  FUNCTION MeasureAsymmetry(this) RESULT(norm_value)
+    !> The matrix to measure
+    TYPE(Matrix_ps), INTENT(IN) :: this
+    !> The norm value of the full distributed sparse matrix.
+    REAL(NTREAL) :: norm_value
+    !! Local variables
+    TYPE(Matrix_ps) :: tmat
+
+    CALL TransposeMatrix(this, tmat)
+    CALL ConjugateMatrix(tmat)
+    CALL IncrementMatrix(this, tmat, alpha_in=-1.0_NTREAL)
+    norm_value = MatrixNorm(tmat)
+    CALL DestructMatrix(tmat)
+
+  END FUNCTION MeasureAsymmetry
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Make the matrix symmetric
+  SUBROUTINE SymmetrizeMatrix(this)
+    !> The matrix to symmetrize
+    TYPE(Matrix_ps), INTENT(INOUT) :: this
+    !! Local variables
+    TYPE(Matrix_ps) :: tmat
+
+    CALL TransposeMatrix(this, tmat)
+    CALL ConjugateMatrix(tmat)
+    CALL IncrementMatrix(tmat, this, alpha_in=1.0_NTREAL)
+    CALL ScaleMatrix(this, 0.5_NTREAL)
+    CALL DestructMatrix(tmat)
+
+  END SUBROUTINE SymmetrizeMatrix
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> Transform a matrix B = P * A * P^-1
   !! This routine will check if P is the identity matrix, and if so

Source/Wrapper/PSMatrixAlgebraModule_wrp.F90

@@ -20,6 +20,8 @@ MODULE PSMatrixAlgebraModule_wrp
   PUBLIC :: ScaleMatrix_ps_wrp
   PUBLIC :: MatrixNorm_ps_wrp
   PUBLIC :: MatrixTrace_ps_wrp
+  PUBLIC :: MeasureAsymmetry_ps_wrp
+  PUBLIC :: SymmetrizeMatrix_ps_wrp
   PUBLIC :: MatrixDiagonalScale_psr_wrp
   PUBLIC :: MatrixDiagonalScale_psc_wrp
 CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -146,6 +148,27 @@ SUBROUTINE MatrixTrace_ps_wrp(ih_this, trace_value) &
     CALL MatrixTrace(h_this%DATA, trace_value)
   END SUBROUTINE MatrixTrace_ps_wrp
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Compute the norm of a distributed sparse matrix along the rows.
+  FUNCTION MeasureAsymmetry_ps_wrp(ih_this) &
+       & BIND(c,name="MeasureAsymmetry_ps_wrp") &
+       & RESULT(norm_value)
+    INTEGER(kind=c_int), INTENT(IN) :: ih_this(SIZE_wrp)
+    REAL(NTREAL) :: norm_value
+    TYPE(Matrix_ps_wrp) :: h_this
+
+    h_this = TRANSFER(ih_this,h_this)
+    norm_value = MeasureAsymmetry(h_this%DATA)
+  END FUNCTION MeasureAsymmetry_ps_wrp
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Symmetrize a matrix
+  SUBROUTINE SymmetrizeMatrix_ps_wrp(ih_matA) &
+       & BIND(c,NAME="SymmetrizeMatrix_ps_wrp")
+    INTEGER(KIND=c_int), INTENT(INOUT) :: ih_matA(SIZE_wrp)
+    TYPE(Matrix_ps_wrp) :: h_matA
+
+    h_matA = TRANSFER(ih_matA,h_matA)
+    CALL SymmetrizeMatrix(h_matA%DATA)
+  END SUBROUTINE SymmetrizeMatrix_ps_wrp
   !> Scale a matrix using a diagonal matrix (triplet list form).
   SUBROUTINE MatrixDiagonalScale_psr_wrp(ih_mat, ih_tlist) &
        & BIND(c,name="MatrixDiagonalScale_psr_wrp")

UnitTests/test_psmatrixalgebra.py

@@ -121,7 +121,8 @@ def check_result(self):
     def check_floats(self, val1, val2):
         from helpers import THRESHOLD
         normval = abs(val1 - val2)
-        self.assertLessEqual(normval, THRESHOLD)
+        global_norm = comm.bcast(normval, root=0)
+        self.assertLessEqual(global_norm, THRESHOLD)
 
     def test_addition_pg(self):
         '''Test routines to add together matrices with an explicit grid.'''
@@ -284,6 +285,49 @@ def test_reverse(self):
 
             self.check_result()
 
+    def test_asymmetry(self):
+        '''Test routines to measure asymmetry of a matrix.'''
+        from numpy import inf
+        from scipy.linalg import norm
+        for param in self.parameters:
+            matrix1 = param.create_matrix(snum=1, complex=self.complex1)
+            self.write_matrix(matrix1, self.input_file1)
+
+            comm.barrier()
+
+            if param.sparsity > 0.0:
+                ntmatrix1 = nt.Matrix_ps(self.input_file1, False)
+            else:
+                ntmatrix1 = nt.Matrix_ps(param.rows)
+
+            diff = matrix1 - matrix1.H
+            ref = norm(diff.todense(), ord=inf)
+            comp = ntmatrix1.MeasureAsymmetry()
+            comm.barrier()
+
+            self.check_floats(ref, comp)
+
+    def test_symmetrize(self):
+        '''Test routines to symmetrize a matrix.'''
+        for param in self.parameters:
+            matrix1 = param.create_matrix(snum=1, complex=self.complex1)
+            self.write_matrix(matrix1, self.input_file1)
+
+            comm.barrier()
+
+            if param.sparsity > 0.0:
+                ntmatrix1 = nt.Matrix_ps(self.input_file1, False)
+            else:
+                ntmatrix1 = nt.Matrix_ps(param.rows)
+
+            self.CheckMat = 0.5 * (matrix1 + matrix1.H)
+            ntmatrix1 = nt.Matrix_ps(self.input_file1, False)
+            ntmatrix1.Symmetrize()
+            ntmatrix1.WriteToMatrixMarket(self.result_file)
+            comm.barrier()
+
+            self.check_result()
+
     def test_scalediag(self):
         '''Test routines to scale by a diagonal matrix.'''
         from copy import deepcopy