Test of asymmetry

Source/C/PSMatrix_c.h

@@ -60,6 +60,7 @@ 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);
 

Source/CPlusPlus/PSMatrix.cc

@@ -227,6 +227,11 @@ double Matrix_ps::Trace() const {
   return temp_val;
 }
 
+//////////////////////////////////////////////////////////////////////////////
+double Matrix_ps::MeasureAsymmetry() const {
+  return MeasureAsymmetry_ps_wrp(ih_this);
+}
+
 //////////////////////////////////////////////////////////////////////////////
 Matrix_ps::~Matrix_ps() { DestructMatrix_ps_wrp(ih_this); }
 } // namespace NTPoly

Source/CPlusPlus/PSMatrix.h

@@ -165,6 +165,8 @@ 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;
 

Source/Fortran/PSMatrixAlgebraModule.F90

@@ -13,7 +13,7 @@ MODULE PSMatrixAlgebraModule
        & ConstructMatrixMemoryPool
   USE PSMatrixModule, ONLY : Matrix_ps, ConstructEmptyMatrix, CopyMatrix, &
        & DestructMatrix, ConvertMatrixToComplex, ConjugateMatrix, &
-       & MergeMatrixLocalBlocks, IsIdentity
+       & MergeMatrixLocalBlocks, IsIdentity, TransposeMatrix
   USE SMatrixAlgebraModule, ONLY : MatrixMultiply, MatrixGrandSum, &
        & PairwiseMultiplyMatrix, IncrementMatrix, ScaleMatrix, &
        & MatrixColumnNorm
@@ -34,6 +34,7 @@ MODULE PSMatrixAlgebraModule
   PUBLIC :: ScaleMatrix
   PUBLIC :: MatrixTrace
   PUBLIC :: SimilarityTransform
+  PUBLIC :: MeasureAsymmetry
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   INTERFACE MatrixSigma
      MODULE PROCEDURE MatrixSigma_ps
@@ -525,6 +526,22 @@ 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)
+
+  END FUNCTION MeasureAsymmetry
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> 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,7 @@ MODULE PSMatrixAlgebraModule_wrp
   PUBLIC :: ScaleMatrix_ps_wrp
   PUBLIC :: MatrixNorm_ps_wrp
   PUBLIC :: MatrixTrace_ps_wrp
+  PUBLIC :: MeasureAsymmetry_ps_wrp
 CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> Matrix B = alpha*Matrix A + Matrix B (AXPY)
   SUBROUTINE IncrementMatrix_ps_wrp(ih_matA, ih_matB, alpha_in,threshold_in) &
@@ -143,5 +144,17 @@ SUBROUTINE MatrixTrace_ps_wrp(ih_this, trace_value) &
     h_this = TRANSFER(ih_this,h_this)
     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
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 END MODULE PSMatrixAlgebraModule_wrp

UnitTests/test_psmatrixalgebra.py

@@ -284,6 +284,27 @@ def test_reverse(self):
 
             self.check_result()
 
+    def test_asymmetry(self):
+        '''Test routines to permute 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()
+
+            self.check_floats(ref, comp)
+
 
 class TestPSMatrixAlgebra_r(TestPSMatrixAlgebra, unittest.TestCase):
     '''Special routines for real algebra'''