From e76d46268f5098814836ca11dca7b0eff24d5822 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Wed, 19 Jun 2024 16:04:38 +0900 Subject: [PATCH 1/6] 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''' From d215bee26862731124c511ef62dc4fbd769715bf Mon Sep 17 00:00:00 2001 From: William Dawson Date: Wed, 19 Jun 2024 16:05:58 +0900 Subject: [PATCH 2/6] Explicit destruct --- Source/Fortran/PSMatrixAlgebraModule.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/Fortran/PSMatrixAlgebraModule.F90 b/Source/Fortran/PSMatrixAlgebraModule.F90 index 9cb793e2..8e05ecd3 100644 --- a/Source/Fortran/PSMatrixAlgebraModule.F90 +++ b/Source/Fortran/PSMatrixAlgebraModule.F90 @@ -540,6 +540,7 @@ FUNCTION MeasureAsymmetry(this) RESULT(norm_value) CALL ConjugateMatrix(tmat) CALL IncrementMatrix(this, tmat, alpha_in=-1.0_NTREAL) norm_value = MatrixNorm(tmat) + CALL DestructMatrix(tmat) END FUNCTION MeasureAsymmetry !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From f365108f6d15301596d527e38ea3c9b63f20a5fd Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 20 Jun 2024 00:29:56 +0900 Subject: [PATCH 3/6] Add symmetrize routine --- Source/C/PSMatrix_c.h | 1 + Source/CPlusPlus/PSMatrix.cc | 5 ++++ Source/CPlusPlus/PSMatrix.h | 2 ++ Source/Fortran/PSMatrixAlgebraModule.F90 | 16 +++++++++++++ Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 | 11 +++++++++ UnitTests/test_psmatrixalgebra.py | 25 +++++++++++++++++++- 6 files changed, 59 insertions(+), 1 deletion(-) diff --git a/Source/C/PSMatrix_c.h b/Source/C/PSMatrix_c.h index 2448272c..3b0e93f8 100644 --- a/Source/C/PSMatrix_c.h +++ b/Source/C/PSMatrix_c.h @@ -63,5 +63,6 @@ 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); #endif diff --git a/Source/CPlusPlus/PSMatrix.cc b/Source/CPlusPlus/PSMatrix.cc index 255a7499..ea4c89e4 100644 --- a/Source/CPlusPlus/PSMatrix.cc +++ b/Source/CPlusPlus/PSMatrix.cc @@ -232,6 +232,11 @@ 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 diff --git a/Source/CPlusPlus/PSMatrix.h b/Source/CPlusPlus/PSMatrix.h index f2a14ff9..0fbf20a5 100644 --- a/Source/CPlusPlus/PSMatrix.h +++ b/Source/CPlusPlus/PSMatrix.h @@ -169,6 +169,8 @@ class Matrix_ps { double MeasureAsymmetry() const; //! compute the trace of a matrix. double Trace() const; + //! symmetrize the matrix + void Symmetrize(); public: //! Destructor. diff --git a/Source/Fortran/PSMatrixAlgebraModule.F90 b/Source/Fortran/PSMatrixAlgebraModule.F90 index 8e05ecd3..72d0b4a0 100644 --- a/Source/Fortran/PSMatrixAlgebraModule.F90 +++ b/Source/Fortran/PSMatrixAlgebraModule.F90 @@ -35,6 +35,7 @@ MODULE PSMatrixAlgebraModule PUBLIC :: MatrixTrace PUBLIC :: SimilarityTransform PUBLIC :: MeasureAsymmetry + PUBLIC :: SymmetrizeMatrix !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INTERFACE MatrixSigma MODULE PROCEDURE MatrixSigma_ps @@ -543,6 +544,21 @@ FUNCTION MeasureAsymmetry(this) RESULT(norm_value) 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 diff --git a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 index aff5b8ed..19deea30 100644 --- a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 +++ b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 @@ -21,6 +21,7 @@ MODULE PSMatrixAlgebraModule_wrp PUBLIC :: MatrixNorm_ps_wrp PUBLIC :: MatrixTrace_ps_wrp PUBLIC :: MeasureAsymmetry_ps_wrp + PUBLIC :: SymmetrizeMatrix_ps_wrp CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Matrix B = alpha*Matrix A + Matrix B (AXPY) SUBROUTINE IncrementMatrix_ps_wrp(ih_matA, ih_matB, alpha_in,threshold_in) & @@ -156,5 +157,15 @@ FUNCTION MeasureAsymmetry_ps_wrp(ih_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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE PSMatrixAlgebraModule_wrp diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py index aeb3cc3b..f9a1882c 100644 --- a/UnitTests/test_psmatrixalgebra.py +++ b/UnitTests/test_psmatrixalgebra.py @@ -285,7 +285,7 @@ def test_reverse(self): self.check_result() def test_asymmetry(self): - '''Test routines to permute a matrix.''' + '''Test routines to measure asymmetry of a matrix.''' from numpy import inf from scipy.linalg import norm for param in self.parameters: @@ -305,6 +305,29 @@ def test_asymmetry(self): self.check_floats(ref, comp) + def test_symmetrize(self): + '''Test routines to symmetrize 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) + + 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() + class TestPSMatrixAlgebra_r(TestPSMatrixAlgebra, unittest.TestCase): '''Special routines for real algebra''' From 2d193faf9c8698f441ade23d221a186cd9eb3cbf Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 20 Jun 2024 10:51:21 +0900 Subject: [PATCH 4/6] Lint --- UnitTests/test_psmatrixalgebra.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py index 2dd90363..a8529b5d 100644 --- a/UnitTests/test_psmatrixalgebra.py +++ b/UnitTests/test_psmatrixalgebra.py @@ -307,8 +307,6 @@ def test_asymmetry(self): def test_symmetrize(self): '''Test routines to symmetrize 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) From daa77a09e3990e6900ba302cf5364976d016d97f Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 20 Jun 2024 13:52:05 +0900 Subject: [PATCH 5/6] lint --- Source/CPlusPlus/PSMatrix.cc | 4 +--- Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/Source/CPlusPlus/PSMatrix.cc b/Source/CPlusPlus/PSMatrix.cc index 7fb74575..473dde8e 100644 --- a/Source/CPlusPlus/PSMatrix.cc +++ b/Source/CPlusPlus/PSMatrix.cc @@ -243,9 +243,7 @@ double Matrix_ps::MeasureAsymmetry() const { } ////////////////////////////////////////////////////////////////////////////// -void Matrix_ps::Symmetrize() { - SymmetrizeMatrix_ps_wrp(ih_this); -} +void Matrix_ps::Symmetrize() { SymmetrizeMatrix_ps_wrp(ih_this); } ////////////////////////////////////////////////////////////////////////////// Matrix_ps::~Matrix_ps() { DestructMatrix_ps_wrp(ih_this); } diff --git a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 index bad51c8f..2c3bd386 100644 --- a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 +++ b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 @@ -162,7 +162,7 @@ END FUNCTION MeasureAsymmetry_ps_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Symmetrize a matrix SUBROUTINE SymmetrizeMatrix_ps_wrp(ih_matA) & - & BIND(c,NAME="SymmetrizeMatrix_ps_wrp") + & BIND(c,NAME="SymmetrizeMatrix_ps_wrp") INTEGER(KIND=c_int), INTENT(INOUT) :: ih_matA(SIZE_wrp) TYPE(Matrix_ps_wrp) :: h_matA From 4564a7d313fef532c6047ba786ede610f8cd580f Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 20 Jun 2024 21:40:51 +0900 Subject: [PATCH 6/6] Fixes to actual testing --- UnitTests/test_psmatrixalgebra.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py index a8529b5d..8726a1a2 100644 --- a/UnitTests/test_psmatrixalgebra.py +++ b/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.''' @@ -302,6 +303,7 @@ def test_asymmetry(self): diff = matrix1 - matrix1.H ref = norm(diff.todense(), ord=inf) comp = ntmatrix1.MeasureAsymmetry() + comm.barrier() self.check_floats(ref, comp) @@ -324,6 +326,8 @@ def test_symmetrize(self): 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