From 235289d61f11f01562455803367eec98ce6869a3 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Wed, 17 Apr 2024 13:59:49 +0900 Subject: [PATCH] Local version of scale with a diagonal --- Source/C/SMatrix_c.h | 4 +++ Source/CPlusPlus/SMatrix.cc | 10 ++++++ Source/CPlusPlus/SMatrix.h | 8 +++++ Source/Fortran/SMatrixAlgebraModule.F90 | 33 +++++++++++++++++++ .../Fortran/sparse_includes/DiagonalScale.f90 | 9 +++++ Source/Wrapper/SMatrixAlgebraModule_wrp.F90 | 32 ++++++++++++++++++ UnitTests/test_matrix.py | 27 +++++++++++++++ 7 files changed, 123 insertions(+) create mode 100644 Source/Fortran/sparse_includes/DiagonalScale.f90 diff --git a/Source/C/SMatrix_c.h b/Source/C/SMatrix_c.h index dcdb4fdc..d2629770 100644 --- a/Source/C/SMatrix_c.h +++ b/Source/C/SMatrix_c.h @@ -34,6 +34,8 @@ void PrintMatrix_lsr_wrp(const int *ih_this); void PrintMatrixF_lsr_wrp(const int *ih_this, const char *file_name, const int *name_size); void MatrixToTripletList_lsr_wrp(const int *ih_this, int *ih_triplet_list); +void DiagonalScale_lsr_wrp(int *ih_mat, const int *ih_tlist, + const double *threhsold); void ConstructMatrixFromFile_lsc_wrp(int *ih_this, const char *file_name, const int *name_size); @@ -70,5 +72,7 @@ void PrintMatrix_lsc_wrp(const int *ih_this); void PrintMatrixF_lsc_wrp(const int *ih_this, const char *file_name, const int *name_size); void MatrixToTripletList_lsc_wrp(const int *ih_this, int *ih_triplet_list); +void DiagonalScale_lsc_wrp(int *ih_mat, const int *ih_tlist, + const double *threhsold); #endif diff --git a/Source/CPlusPlus/SMatrix.cc b/Source/CPlusPlus/SMatrix.cc index 521d6546..68620ae4 100644 --- a/Source/CPlusPlus/SMatrix.cc +++ b/Source/CPlusPlus/SMatrix.cc @@ -175,6 +175,16 @@ void Matrix_lsc::Gemm(const Matrix_lsc &matA, const Matrix_lsc &matB, memory_pool.ih_this); } +//////////////////////////////////////////////////////////////////////////////// +void Matrix_lsr::DiagonalScale(const TripletList_r &tlist, double threshold) { + DiagonalScale_lsr_wrp(ih_this, tlist.ih_this, &threshold); +} + +//////////////////////////////////////////////////////////////////////////////// +void Matrix_lsc::DiagonalScale(const TripletList_c &tlist, double threshold) { + DiagonalScale_lsc_wrp(ih_this, tlist.ih_this, &threshold); +} + //////////////////////////////////////////////////////////////////////////////// void Matrix_lsr::Transpose(const Matrix_lsr &matA) { TransposeMatrix_lsr_wrp(matA.ih_this, ih_this); diff --git a/Source/CPlusPlus/SMatrix.h b/Source/CPlusPlus/SMatrix.h index fd5079c1..40984449 100644 --- a/Source/CPlusPlus/SMatrix.h +++ b/Source/CPlusPlus/SMatrix.h @@ -76,6 +76,10 @@ class Matrix_lsr { void Gemm(const NTPoly::Matrix_lsr &matA, const NTPoly::Matrix_lsr &matB, bool isATransposed, bool isBTransposed, double alpha, double beta, double threshold, NTPoly::MatrixMemoryPool_r &memory_pool); + //! Scale a matrix using a diagonal matrix (triplet list form). + //!\param tlist the triplet list. + //!\param threshold for flushing small values. + void DiagonalScale(const NTPoly::TripletList_r &tlist, double threshold); public: //! Transpose a sparse matrix. @@ -170,6 +174,10 @@ class Matrix_lsc { void Gemm(const NTPoly::Matrix_lsc &matA, const NTPoly::Matrix_lsc &matB, bool isATransposed, bool isBTransposed, double alpha, double beta, double threshold, NTPoly::MatrixMemoryPool_c &memory_pool); + //! Scale a matrix using a diagonal matrix (triplet list form). + //!\param tlist the triplet list. + //!\param threshold for flushing small values. + void DiagonalScale(const NTPoly::TripletList_c &tlist, double threshold); public: //! Transpose a sparse matrix. diff --git a/Source/Fortran/SMatrixAlgebraModule.F90 b/Source/Fortran/SMatrixAlgebraModule.F90 index cd790e76..1493a300 100644 --- a/Source/Fortran/SMatrixAlgebraModule.F90 +++ b/Source/Fortran/SMatrixAlgebraModule.F90 @@ -25,6 +25,7 @@ MODULE SMatrixAlgebraModule PUBLIC :: MatrixColumnNorm PUBLIC :: MatrixNorm PUBLIC :: MatrixGrandSum + PUBLIC :: DiagonalScale !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INTERFACE ScaleMatrix MODULE PROCEDURE ScaleMatrix_lsr @@ -75,6 +76,10 @@ MODULE SMatrixAlgebraModule MODULE PROCEDURE DenseBranch_lsr MODULE PROCEDURE DenseBranch_lsc END INTERFACE DenseBranch + INTERFACE DiagonalScale + MODULE PROCEDURE DiagonalScale_lsr + MODULE PROCEDURE DiagonalScale_lsc + END INTERFACE DiagonalScale CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Will scale a sparse matrix by a constant. PURE SUBROUTINE ScaleMatrix_lsr(matA,constant) @@ -526,5 +531,33 @@ PURE SUBROUTINE PruneList_lsc(memorypool,alpha,threshold, & #include "sparse_includes/PruneList.f90" END SUBROUTINE PruneList_lsc +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Scale a matrix using a diagonal matrix (triplet list form). + SUBROUTINE DiagonalScale_lsr(mat, tlist, threshold_in) + !> The matrix to scale. + TYPE(Matrix_lsr), INTENT(INOUT) :: mat + !> The diagonal matrix. + TYPE(TripletList_r), INTENT(IN) :: tlist + !> For flushing values to zero. Default value is 0.0. + REAL(NTREAL), OPTIONAL, INTENT(IN) :: threshold_in + !! Intermediate Data + REAL(NTREAL) :: val + +#include "sparse_includes/DiagonalScale.f90" + END SUBROUTINE DiagonalScale_lsr +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Scale a matrix using a diagonal matrix (triplet list form). + SUBROUTINE DiagonalScale_lsc(mat, tlist, threshold_in) + !> The matrix to scale. + TYPE(Matrix_lsc), INTENT(INOUT) :: mat + !> The diagonal matrix. + TYPE(TripletList_c), INTENT(IN) :: tlist + !> For flushing values to zero. Default value is 0.0. + REAL(NTREAL), OPTIONAL, INTENT(IN) :: threshold_in + !! Intermediate Data + COMPLEX(NTCOMPLEX) :: val + +#include "sparse_includes/DiagonalScale.f90" + END SUBROUTINE DiagonalScale_lsc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE SMatrixAlgebraModule diff --git a/Source/Fortran/sparse_includes/DiagonalScale.f90 b/Source/Fortran/sparse_includes/DiagonalScale.f90 new file mode 100644 index 00000000..c1997b5b --- /dev/null +++ b/Source/Fortran/sparse_includes/DiagonalScale.f90 @@ -0,0 +1,9 @@ + +INTEGER :: col, II + +DO II = 1, tlist%CurrentSize + col = tlist%DATA(II)%index_column + val = tlist%DATA(II)%point_value + mat%values(mat%outer_index(col) + 1:mat%outer_index(col + 1)) = & + val * mat%values(mat%outer_index(col) + 1:mat%outer_index(col + 1)) +END DO diff --git a/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 index 0292018a..55fec023 100644 --- a/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 +++ b/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 @@ -18,12 +18,14 @@ MODULE SMatrixAlgebraModule_wrp PUBLIC :: DotMatrix_lsr_wrp PUBLIC :: PairwiseMultiplyMatrix_lsr_wrp PUBLIC :: MatrixMultiply_lsr_wrp + PUBLIC :: DiagonalScale_lsr_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PUBLIC :: ScaleMatrix_lsc_wrp PUBLIC :: IncrementMatrix_lsc_wrp PUBLIC :: DotMatrix_lsc_wrp PUBLIC :: PairwiseMultiplyMatrix_lsc_wrp PUBLIC :: MatrixMultiply_lsc_wrp + PUBLIC :: DiagonalScale_lsc_wrp CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Wrap the scale a sparse matrix by a constant routine. SUBROUTINE ScaleMatrix_lsr_wrp(ih_this, constant) & @@ -110,6 +112,21 @@ SUBROUTINE MatrixMultiply_lsr_wrp(ih_matA, ih_matB, ih_matC, IsATransposed, & & LOGICAL(IsATransposed), LOGICAL(IsBTransposed), alpha, & & beta, threshold, h_blocked_memory_pool%DATA) END SUBROUTINE MatrixMultiply_lsr_wrp +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Scale a matrix using a diagonal matrix (triplet list form). + SUBROUTINE DiagonalScale_lsr_wrp(ih_mat, ih_tlist, threshold) & + & BIND(c,name="DiagonalScale_lsr_wrp") + INTEGER(kind=c_int), INTENT(INOUT) :: ih_mat(SIZE_wrp) + INTEGER(kind=c_int), INTENT(IN) :: ih_tlist(SIZE_wrp) + REAL(NTREAL), INTENT(in) :: threshold + TYPE(Matrix_lsr_wrp) :: h_mat + TYPE(TripletList_r_wrp) :: h_tlist + + h_mat = TRANSFER(ih_mat, h_mat) + h_tlist = TRANSFER(ih_tlist, h_tlist) + + CALL DiagonalScale(h_mat%DATA, h_tlist%DATA, threshold) + END SUBROUTINE DiagonalScale_lsr_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Wrap the scale a sparse matrix by a constant routine. SUBROUTINE ScaleMatrix_lsc_wrp(ih_this, constant) & @@ -201,5 +218,20 @@ SUBROUTINE MatrixMultiply_lsc_wrp(ih_matA, ih_matB, ih_matC, IsATransposed, & & LOGICAL(IsATransposed), LOGICAL(IsBTransposed), alpha, & & beta, threshold, h_blocked_memory_pool%DATA) END SUBROUTINE MatrixMultiply_lsc_wrp +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Scale a matrix using a diagonal matrix (triplet list form). + SUBROUTINE DiagonalScale_lsc_wrp(ih_mat, ih_tlist, threshold) & + & BIND(c,name="DiagonalScale_lsc_wrp") + INTEGER(kind=c_int), INTENT(INOUT) :: ih_mat(SIZE_wrp) + INTEGER(kind=c_int), INTENT(IN) :: ih_tlist(SIZE_wrp) + REAL(NTREAL), INTENT(in) :: threshold + TYPE(Matrix_lsc_wrp) :: h_mat + TYPE(TripletList_c_wrp) :: h_tlist + + h_mat = TRANSFER(ih_mat, h_mat) + h_tlist = TRANSFER(ih_tlist, h_tlist) + + CALL DiagonalScale(h_mat%DATA, h_tlist%DATA, threshold) + END SUBROUTINE DiagonalScale_lsc_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE SMatrixAlgebraModule_wrp diff --git a/UnitTests/test_matrix.py b/UnitTests/test_matrix.py index 585daa34..2d874988 100644 --- a/UnitTests/test_matrix.py +++ b/UnitTests/test_matrix.py @@ -57,6 +57,8 @@ class TestLocalMatrix(unittest.TestCase): SMatrix = nt.Matrix_lsr MatrixMemoryPool = nt.MatrixMemoryPool_r complex = False + TripletList = nt.TripletList_r + Triplet = nt.Triplet_r def _compare_mat(self, val1, val2): from helpers import THRESHOLD @@ -397,12 +399,37 @@ def test_get_column(self): ResultMat = mmread(self.file2) self._compare_mat(CheckMat, ResultMat) + def test_scalediag(self): + '''Test routines to scale by a diagonal matrix.''' + from copy import deepcopy + for param in self.parameters: + matrix = param.create_matrix(complex=self.complex) + mmwrite(self.file1, matrix) + CheckMat = deepcopy(matrix) + + tlist = self.TripletList() + for i in range(matrix.shape[1]): + t = self.Triplet() + t.index_column = i + 1 + t.index_row = i + 1 + t.point_value = i + tlist.Append(t) + CheckMat[:, i] *= i + ntmatrix = self.SMatrix(self.file1) + ntmatrix.DiagonalScale(tlist, 0) + ntmatrix.WriteToMatrixMarket(self.file2) + + ResultMat = mmread(self.file2) + self._compare_mat(CheckMat, ResultMat) + class TestLocalMatrix_c(TestLocalMatrix): '''Specialization for complex matrices''' SMatrix = nt.Matrix_lsc MatrixMemoryPool = nt.MatrixMemoryPool_c complex = True + TripletList = nt.TripletList_c + Triplet = nt.Triplet_c def test_conjugatetranspose(self): '''Test routines to compute the conjugate transpose of a matrix.'''