From e76d46268f5098814836ca11dca7b0eff24d5822 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Wed, 19 Jun 2024 16:04:38 +0900 Subject: [PATCH] Test of asymmetry --- Source/C/PSMatrix_c.h | 1 + Source/CPlusPlus/PSMatrix.cc | 5 +++++ Source/CPlusPlus/PSMatrix.h | 2 ++ Source/Fortran/PSMatrixAlgebraModule.F90 | 19 +++++++++++++++++- Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 | 13 ++++++++++++ UnitTests/test_psmatrixalgebra.py | 21 ++++++++++++++++++++ 6 files changed, 60 insertions(+), 1 deletion(-) diff --git a/Source/C/PSMatrix_c.h b/Source/C/PSMatrix_c.h index b351b217..2448272c 100644 --- a/Source/C/PSMatrix_c.h +++ b/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); diff --git a/Source/CPlusPlus/PSMatrix.cc b/Source/CPlusPlus/PSMatrix.cc index c2c296bd..255a7499 100644 --- a/Source/CPlusPlus/PSMatrix.cc +++ b/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 diff --git a/Source/CPlusPlus/PSMatrix.h b/Source/CPlusPlus/PSMatrix.h index eb2bd2ab..f2a14ff9 100644 --- a/Source/CPlusPlus/PSMatrix.h +++ b/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; diff --git a/Source/Fortran/PSMatrixAlgebraModule.F90 b/Source/Fortran/PSMatrixAlgebraModule.F90 index c2b13f86..9cb793e2 100644 --- a/Source/Fortran/PSMatrixAlgebraModule.F90 +++ b/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 diff --git a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 index 3fe86b2d..aff5b8ed 100644 --- a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 +++ b/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 diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py index 1d09c88b..aeb3cc3b 100644 --- a/UnitTests/test_psmatrixalgebra.py +++ b/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'''