From 235289d61f11f01562455803367eec98ce6869a3 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Wed, 17 Apr 2024 13:59:49 +0900 Subject: [PATCH 01/16] 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.''' From 42a94c6505251b738ae0bda21b33df84c0c0d243 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Wed, 17 Apr 2024 14:12:01 +0900 Subject: [PATCH 02/16] Remove unused parameter --- Source/C/SMatrix_c.h | 6 ++---- Source/CPlusPlus/SMatrix.cc | 6 +++--- Source/CPlusPlus/SMatrix.h | 2 +- Source/Fortran/SMatrixAlgebraModule.F90 | 8 ++------ Source/Wrapper/SMatrixAlgebraModule_wrp.F90 | 10 ++++------ UnitTests/test_matrix.py | 2 +- 6 files changed, 13 insertions(+), 21 deletions(-) diff --git a/Source/C/SMatrix_c.h b/Source/C/SMatrix_c.h index d2629770..4a49f801 100644 --- a/Source/C/SMatrix_c.h +++ b/Source/C/SMatrix_c.h @@ -34,8 +34,7 @@ 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 DiagonalScale_lsr_wrp(int *ih_mat, const int *ih_tlist); void ConstructMatrixFromFile_lsc_wrp(int *ih_this, const char *file_name, const int *name_size); @@ -72,7 +71,6 @@ 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); +void DiagonalScale_lsc_wrp(int *ih_mat, const int *ih_tlist); #endif diff --git a/Source/CPlusPlus/SMatrix.cc b/Source/CPlusPlus/SMatrix.cc index 68620ae4..48dbbe8d 100644 --- a/Source/CPlusPlus/SMatrix.cc +++ b/Source/CPlusPlus/SMatrix.cc @@ -176,13 +176,13 @@ void Matrix_lsc::Gemm(const Matrix_lsc &matA, const Matrix_lsc &matB, } //////////////////////////////////////////////////////////////////////////////// -void Matrix_lsr::DiagonalScale(const TripletList_r &tlist, double threshold) { - DiagonalScale_lsr_wrp(ih_this, tlist.ih_this, &threshold); +void Matrix_lsr::DiagonalScale(const TripletList_r &tlist) { + DiagonalScale_lsr_wrp(ih_this, tlist.ih_this); } //////////////////////////////////////////////////////////////////////////////// void Matrix_lsc::DiagonalScale(const TripletList_c &tlist, double threshold) { - DiagonalScale_lsc_wrp(ih_this, tlist.ih_this, &threshold); + DiagonalScale_lsc_wrp(ih_this, tlist.ih_this); } //////////////////////////////////////////////////////////////////////////////// diff --git a/Source/CPlusPlus/SMatrix.h b/Source/CPlusPlus/SMatrix.h index 40984449..21f90116 100644 --- a/Source/CPlusPlus/SMatrix.h +++ b/Source/CPlusPlus/SMatrix.h @@ -79,7 +79,7 @@ class Matrix_lsr { //! 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); + void DiagonalScale(const NTPoly::TripletList_r &tlist); public: //! Transpose a sparse matrix. diff --git a/Source/Fortran/SMatrixAlgebraModule.F90 b/Source/Fortran/SMatrixAlgebraModule.F90 index 1493a300..ace4bbc9 100644 --- a/Source/Fortran/SMatrixAlgebraModule.F90 +++ b/Source/Fortran/SMatrixAlgebraModule.F90 @@ -533,13 +533,11 @@ PURE SUBROUTINE PruneList_lsc(memorypool,alpha,threshold, & END SUBROUTINE PruneList_lsc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE DiagonalScale_lsr(mat, tlist, threshold_in) + SUBROUTINE DiagonalScale_lsr(mat, tlist) !> 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 @@ -547,13 +545,11 @@ SUBROUTINE DiagonalScale_lsr(mat, tlist, threshold_in) END SUBROUTINE DiagonalScale_lsr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE DiagonalScale_lsc(mat, tlist, threshold_in) + SUBROUTINE DiagonalScale_lsc(mat, tlist) !> 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 diff --git a/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 index 55fec023..405bd326 100644 --- a/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 +++ b/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 @@ -114,18 +114,17 @@ SUBROUTINE MatrixMultiply_lsr_wrp(ih_matA, ih_matB, ih_matC, IsATransposed, & END SUBROUTINE MatrixMultiply_lsr_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE DiagonalScale_lsr_wrp(ih_mat, ih_tlist, threshold) & + SUBROUTINE DiagonalScale_lsr_wrp(ih_mat, ih_tlist) & & 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) + CALL DiagonalScale(h_mat%DATA, h_tlist%DATA) END SUBROUTINE DiagonalScale_lsr_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Wrap the scale a sparse matrix by a constant routine. @@ -220,18 +219,17 @@ SUBROUTINE MatrixMultiply_lsc_wrp(ih_matA, ih_matB, ih_matC, IsATransposed, & END SUBROUTINE MatrixMultiply_lsc_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE DiagonalScale_lsc_wrp(ih_mat, ih_tlist, threshold) & + SUBROUTINE DiagonalScale_lsc_wrp(ih_mat, ih_tlist) & & 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) + CALL DiagonalScale(h_mat%DATA, h_tlist%DATA) END SUBROUTINE DiagonalScale_lsc_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE SMatrixAlgebraModule_wrp diff --git a/UnitTests/test_matrix.py b/UnitTests/test_matrix.py index 2d874988..90ee3d14 100644 --- a/UnitTests/test_matrix.py +++ b/UnitTests/test_matrix.py @@ -416,7 +416,7 @@ def test_scalediag(self): tlist.Append(t) CheckMat[:, i] *= i ntmatrix = self.SMatrix(self.file1) - ntmatrix.DiagonalScale(tlist, 0) + ntmatrix.DiagonalScale(tlist) ntmatrix.WriteToMatrixMarket(self.file2) ResultMat = mmread(self.file2) From 74a54254f2b566b56f358dd1dc96b52e7028f0b9 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Wed, 17 Apr 2024 14:36:31 +0900 Subject: [PATCH 03/16] Fix API --- Source/C/SMatrix_c.h | 4 ++-- Source/CPlusPlus/SMatrix.cc | 4 ++-- Source/Fortran/SMatrixAlgebraModule.F90 | 18 +++++++++--------- Source/Wrapper/SMatrixAlgebraModule_wrp.F90 | 20 ++++++++++---------- 4 files changed, 23 insertions(+), 23 deletions(-) diff --git a/Source/C/SMatrix_c.h b/Source/C/SMatrix_c.h index 4a49f801..9307c931 100644 --- a/Source/C/SMatrix_c.h +++ b/Source/C/SMatrix_c.h @@ -34,7 +34,7 @@ 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); +void MatrixDiagonalScale_lsr_wrp(int *ih_mat, const int *ih_tlist); void ConstructMatrixFromFile_lsc_wrp(int *ih_this, const char *file_name, const int *name_size); @@ -71,6 +71,6 @@ 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); +void MatrixDiagonalScale_lsc_wrp(int *ih_mat, const int *ih_tlist); #endif diff --git a/Source/CPlusPlus/SMatrix.cc b/Source/CPlusPlus/SMatrix.cc index 48dbbe8d..588677f6 100644 --- a/Source/CPlusPlus/SMatrix.cc +++ b/Source/CPlusPlus/SMatrix.cc @@ -177,12 +177,12 @@ void Matrix_lsc::Gemm(const Matrix_lsc &matA, const Matrix_lsc &matB, //////////////////////////////////////////////////////////////////////////////// void Matrix_lsr::DiagonalScale(const TripletList_r &tlist) { - DiagonalScale_lsr_wrp(ih_this, tlist.ih_this); + MatrixDiagonalScale_lsr_wrp(ih_this, tlist.ih_this); } //////////////////////////////////////////////////////////////////////////////// void Matrix_lsc::DiagonalScale(const TripletList_c &tlist, double threshold) { - DiagonalScale_lsc_wrp(ih_this, tlist.ih_this); + MatrixDiagonalScale_lsc_wrp(ih_this, tlist.ih_this); } //////////////////////////////////////////////////////////////////////////////// diff --git a/Source/Fortran/SMatrixAlgebraModule.F90 b/Source/Fortran/SMatrixAlgebraModule.F90 index ace4bbc9..ce6415c4 100644 --- a/Source/Fortran/SMatrixAlgebraModule.F90 +++ b/Source/Fortran/SMatrixAlgebraModule.F90 @@ -25,7 +25,7 @@ MODULE SMatrixAlgebraModule PUBLIC :: MatrixColumnNorm PUBLIC :: MatrixNorm PUBLIC :: MatrixGrandSum - PUBLIC :: DiagonalScale + PUBLIC :: MatrixDiagonalScale !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INTERFACE ScaleMatrix MODULE PROCEDURE ScaleMatrix_lsr @@ -76,10 +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 + INTERFACE MatrixDiagonalScale + MODULE PROCEDURE MatrixDiagonalScale_lsr + MODULE PROCEDURE MatrixDiagonalScale_lsc + END INTERFACE MatrixDiagonalScale CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Will scale a sparse matrix by a constant. PURE SUBROUTINE ScaleMatrix_lsr(matA,constant) @@ -533,7 +533,7 @@ PURE SUBROUTINE PruneList_lsc(memorypool,alpha,threshold, & END SUBROUTINE PruneList_lsc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE DiagonalScale_lsr(mat, tlist) + SUBROUTINE MatrixDiagonalScale_lsr(mat, tlist) !> The matrix to scale. TYPE(Matrix_lsr), INTENT(INOUT) :: mat !> The diagonal matrix. @@ -542,10 +542,10 @@ SUBROUTINE DiagonalScale_lsr(mat, tlist) REAL(NTREAL) :: val #include "sparse_includes/DiagonalScale.f90" - END SUBROUTINE DiagonalScale_lsr + END SUBROUTINE MatrixDiagonalScale_lsr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE DiagonalScale_lsc(mat, tlist) + SUBROUTINE MatrixDiagonalScale_lsc(mat, tlist) !> The matrix to scale. TYPE(Matrix_lsc), INTENT(INOUT) :: mat !> The diagonal matrix. @@ -554,6 +554,6 @@ SUBROUTINE DiagonalScale_lsc(mat, tlist) COMPLEX(NTCOMPLEX) :: val #include "sparse_includes/DiagonalScale.f90" - END SUBROUTINE DiagonalScale_lsc + END SUBROUTINE MatrixDiagonalScale_lsc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE SMatrixAlgebraModule diff --git a/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 index 405bd326..03b5d0ed 100644 --- a/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 +++ b/Source/Wrapper/SMatrixAlgebraModule_wrp.F90 @@ -18,14 +18,14 @@ MODULE SMatrixAlgebraModule_wrp PUBLIC :: DotMatrix_lsr_wrp PUBLIC :: PairwiseMultiplyMatrix_lsr_wrp PUBLIC :: MatrixMultiply_lsr_wrp - PUBLIC :: DiagonalScale_lsr_wrp + PUBLIC :: MatrixDiagonalScale_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 + PUBLIC :: MatrixDiagonalScale_lsc_wrp CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Wrap the scale a sparse matrix by a constant routine. SUBROUTINE ScaleMatrix_lsr_wrp(ih_this, constant) & @@ -114,8 +114,8 @@ SUBROUTINE MatrixMultiply_lsr_wrp(ih_matA, ih_matB, ih_matC, IsATransposed, & END SUBROUTINE MatrixMultiply_lsr_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE DiagonalScale_lsr_wrp(ih_mat, ih_tlist) & - & BIND(c,name="DiagonalScale_lsr_wrp") + SUBROUTINE MatrixDiagonalScale_lsr_wrp(ih_mat, ih_tlist) & + & BIND(c,name="MatrixDiagonalScale_lsr_wrp") INTEGER(kind=c_int), INTENT(INOUT) :: ih_mat(SIZE_wrp) INTEGER(kind=c_int), INTENT(IN) :: ih_tlist(SIZE_wrp) TYPE(Matrix_lsr_wrp) :: h_mat @@ -124,8 +124,8 @@ SUBROUTINE DiagonalScale_lsr_wrp(ih_mat, ih_tlist) & h_mat = TRANSFER(ih_mat, h_mat) h_tlist = TRANSFER(ih_tlist, h_tlist) - CALL DiagonalScale(h_mat%DATA, h_tlist%DATA) - END SUBROUTINE DiagonalScale_lsr_wrp + CALL MatrixDiagonalScale(h_mat%DATA, h_tlist%DATA) + END SUBROUTINE MatrixDiagonalScale_lsr_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Wrap the scale a sparse matrix by a constant routine. SUBROUTINE ScaleMatrix_lsc_wrp(ih_this, constant) & @@ -219,8 +219,8 @@ SUBROUTINE MatrixMultiply_lsc_wrp(ih_matA, ih_matB, ih_matC, IsATransposed, & END SUBROUTINE MatrixMultiply_lsc_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE DiagonalScale_lsc_wrp(ih_mat, ih_tlist) & - & BIND(c,name="DiagonalScale_lsc_wrp") + SUBROUTINE MatrixDiagonalScale_lsc_wrp(ih_mat, ih_tlist) & + & BIND(c,name="MatrixDiagonalScale_lsc_wrp") INTEGER(kind=c_int), INTENT(INOUT) :: ih_mat(SIZE_wrp) INTEGER(kind=c_int), INTENT(IN) :: ih_tlist(SIZE_wrp) TYPE(Matrix_lsc_wrp) :: h_mat @@ -229,7 +229,7 @@ SUBROUTINE DiagonalScale_lsc_wrp(ih_mat, ih_tlist) & h_mat = TRANSFER(ih_mat, h_mat) h_tlist = TRANSFER(ih_tlist, h_tlist) - CALL DiagonalScale(h_mat%DATA, h_tlist%DATA) - END SUBROUTINE DiagonalScale_lsc_wrp + CALL MatrixDiagonalScale(h_mat%DATA, h_tlist%DATA) + END SUBROUTINE MatrixDiagonalScale_lsc_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE SMatrixAlgebraModule_wrp From 15045779197f6c79670142b065e53780e3848421 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Wed, 17 Apr 2024 17:32:11 +0900 Subject: [PATCH 04/16] First attempt at parallel, have to test somewhere else... --- Source/C/PSMatrix_c.h | 1 + Source/CPlusPlus/PSMatrix.cc | 5 +++ Source/CPlusPlus/PSMatrix.h | 5 ++- Source/Fortran/PSMatrixAlgebraModule.F90 | 42 +++++++++++++++++-- Source/Fortran/PSMatrixModule.F90 | 4 +- .../ScaleDiagonal.f90 | 24 +++++++++++ Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 | 15 +++++++ UnitTests/test_psmatrixalgebra.py | 31 ++++++++++++++ 8 files changed, 121 insertions(+), 6 deletions(-) create mode 100644 Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 diff --git a/Source/C/PSMatrix_c.h b/Source/C/PSMatrix_c.h index b351b217..c07bd7a0 100644 --- a/Source/C/PSMatrix_c.h +++ b/Source/C/PSMatrix_c.h @@ -62,5 +62,6 @@ void ScaleMatrix_ps_wrp(int *ih_this, const double *constant); double MatrixNorm_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 MatrixDiagonalScale_ps_wrp(int *ih_mat, const int *ih_tlist); #endif diff --git a/Source/CPlusPlus/PSMatrix.cc b/Source/CPlusPlus/PSMatrix.cc index c2c296bd..d421c4df 100644 --- a/Source/CPlusPlus/PSMatrix.cc +++ b/Source/CPlusPlus/PSMatrix.cc @@ -217,6 +217,11 @@ void Matrix_ps::Scale(double constant) { ScaleMatrix_ps_wrp(ih_this, &constant); } +////////////////////////////////////////////////////////////////////////////// +void Matrix_ps::DiagonalScale(const TripletList_r &tlist) { + MatrixDiagonalScale_ps_wrp(ih_this, tlist.ih_this); +} + ////////////////////////////////////////////////////////////////////////////// double Matrix_ps::Norm() const { return MatrixNorm_ps_wrp(ih_this); } diff --git a/Source/CPlusPlus/PSMatrix.h b/Source/CPlusPlus/PSMatrix.h index eb2bd2ab..d6e3bc90 100644 --- a/Source/CPlusPlus/PSMatrix.h +++ b/Source/CPlusPlus/PSMatrix.h @@ -160,13 +160,16 @@ class Matrix_ps { void Gemm(const Matrix_ps &matA, const Matrix_ps &matB, PMatrixMemoryPool &memory_pool, double alpha = 1.0, double beta = 0.0, double threshold = 0.0); - //! scale the matrix by a constatn. + //! scale the matrix by a constant. //! constant the value to scale by. void Scale(double constant); //! compute the norm of a matrix. double Norm() const; //! compute the trace of a matrix. double Trace() const; + //!\param tlist the triplet list. + //!\param threshold for flushing small values. + void DiagonalScale(const NTPoly::TripletList_r &tlist); public: //! Destructor. diff --git a/Source/Fortran/PSMatrixAlgebraModule.F90 b/Source/Fortran/PSMatrixAlgebraModule.F90 index c2b13f86..6ecea2da 100644 --- a/Source/Fortran/PSMatrixAlgebraModule.F90 +++ b/Source/Fortran/PSMatrixAlgebraModule.F90 @@ -13,13 +13,16 @@ MODULE PSMatrixAlgebraModule & ConstructMatrixMemoryPool USE PSMatrixModule, ONLY : Matrix_ps, ConstructEmptyMatrix, CopyMatrix, & & DestructMatrix, ConvertMatrixToComplex, ConjugateMatrix, & - & MergeMatrixLocalBlocks, IsIdentity + & MergeMatrixLocalBlocks, IsIdentity, SplitMatrixToLocalBlocks USE SMatrixAlgebraModule, ONLY : MatrixMultiply, MatrixGrandSum, & & PairwiseMultiplyMatrix, IncrementMatrix, ScaleMatrix, & - & MatrixColumnNorm + & MatrixColumnNorm, MatrixDiagonalScale USE SMatrixModule, ONLY : Matrix_lsr, Matrix_lsc, DestructMatrix, CopyMatrix,& & TransposeMatrix, ComposeMatrixColumns, MatrixToTripletList - USE TripletListModule, ONLY : TripletList_r, TripletList_c + USE TripletListModule, ONLY : TripletList_r, TripletList_c, & + & ConstructTripletList, AppendToTripletList, DestructTripletList, & + & GetTripletAt + USE TripletModule, ONLY : Triplet_r, Triplet_c USE NTMPIModule IMPLICIT NONE PRIVATE @@ -34,6 +37,7 @@ MODULE PSMatrixAlgebraModule PUBLIC :: ScaleMatrix PUBLIC :: MatrixTrace PUBLIC :: SimilarityTransform + PUBLIC :: MatrixDiagonalScale !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INTERFACE MatrixSigma MODULE PROCEDURE MatrixSigma_ps @@ -62,6 +66,10 @@ MODULE PSMatrixAlgebraModule MODULE PROCEDURE ScaleMatrix_psr MODULE PROCEDURE ScaleMatrix_psc END INTERFACE ScaleMatrix + INTERFACE MatrixDiagonalScale + MODULE PROCEDURE MatrixDiagonalScale_psr + MODULE PROCEDURE MatrixDiagonalScale_psc + END INTERFACE MatrixDiagonalScale INTERFACE MatrixTrace MODULE PROCEDURE MatrixTrace_psr END INTERFACE MatrixTrace @@ -491,6 +499,34 @@ RECURSIVE SUBROUTINE ScaleMatrix_psc(this, constant) END IF END SUBROUTINE ScaleMatrix_psc +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Will scale a distributed sparse matrix by a constant. + SUBROUTINE MatrixDiagonalScale_psr(this, tlist) + !> Matrix to scale. + TYPE(Matrix_ps), INTENT(INOUT) :: this + !> A constant scale factor. + TYPE(TripletList_r), INTENT(IN) :: tlist + !! Local Data + TYPE(Matrix_lsr) :: lmat + TYPE(TripletList_r) :: filtered + TYPE(Triplet_r) :: trip + +#include "distributed_algebra_includes/ScaleDiagonal.f90" + END SUBROUTINE MatrixDiagonalScale_psr +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Will scale a distributed sparse matrix by a constant. + RECURSIVE SUBROUTINE MatrixDiagonalScale_psc(this, tlist) + !> Matrix to scale. + TYPE(Matrix_ps), INTENT(INOUT) :: this + !> A constant scale factor. + TYPE(TripletList_c), INTENT(IN) :: tlist + !! Local Data + TYPE(Matrix_lsc) :: lmat + TYPE(TripletList_c) :: filtered + TYPE(Triplet_c) :: trip + +#include "distributed_algebra_includes/ScaleDiagonal.f90" + END SUBROUTINE MatrixDiagonalScale_psc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Compute the trace of the matrix. SUBROUTINE MatrixTrace_psr(this, trace_value) diff --git a/Source/Fortran/PSMatrixModule.F90 b/Source/Fortran/PSMatrixModule.F90 index 19812a28..423254e8 100644 --- a/Source/Fortran/PSMatrixModule.F90 +++ b/Source/Fortran/PSMatrixModule.F90 @@ -1443,7 +1443,7 @@ SUBROUTINE TransposeMatrix_psr(AMat, TransMat) !! Local Variables TYPE(TripletList_r) :: tlist TYPE(TripletList_r) :: new_list - TYPE(Triplet_r) :: trip, trip_t + TYPE(Triplet_r) :: trip #include "distributed_includes/TransposeMatrix.f90" @@ -1458,7 +1458,7 @@ SUBROUTINE TransposeMatrix_psc(AMat, TransMat) !! Local Variables TYPE(TripletList_c) :: tlist TYPE(TripletList_c) :: new_list - TYPE(Triplet_c) :: trip, trip_t + TYPE(Triplet_c) :: trip #include "distributed_includes/TransposeMatrix.f90" diff --git a/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 b/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 new file mode 100644 index 00000000..ee3b4db5 --- /dev/null +++ b/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 @@ -0,0 +1,24 @@ +INTEGER :: II, row + +!! Merge to the local block +CALL MergeMatrixLocalBlocks(this, lmat) + +!! Filter out the triplets that aren't stored locally +CALL ConstructTripletList(filtered) +DO II = 1, tlist%CurrentSize + CALL GetTripletAt(tlist, II, trip) + row = trip%index_row + IF (row .GE. this%start_row .AND. row .LT. this%end_row) THEN + trip%index_row = trip%index_row - this%start_row + 1 + trip%index_column = trip%index_row + CALL AppendToTripletList(filtered, trip) + END IF +END DO + +!! Scale +CALL MatrixDiagonalScale(lmat, filtered) + +!! Split +CALL SplitMatrixToLocalBlocks(this, lmat) +CALL DestructMatrix(lmat) +CALL DestructTripletList(filtered) \ No newline at end of file diff --git a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 index 3fe86b2d..e333d6ac 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 :: MatrixDiagonalScale_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,19 @@ 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 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Scale a matrix using a diagonal matrix (triplet list form). + SUBROUTINE MatrixDiagonalScale_ps_wrp(ih_mat, ih_tlist) & + & BIND(c,name="MatrixDiagonalScale_ps_wrp") + INTEGER(kind=c_int), INTENT(INOUT) :: ih_mat(SIZE_wrp) + INTEGER(kind=c_int), INTENT(IN) :: ih_tlist(SIZE_wrp) + TYPE(Matrix_ps_wrp) :: h_mat + TYPE(TripletList_r_wrp) :: h_tlist + + h_mat = TRANSFER(ih_mat, h_mat) + h_tlist = TRANSFER(ih_tlist, h_tlist) + + CALL MatrixDiagonalScale(h_mat%DATA, h_tlist%DATA) + END SUBROUTINE MatrixDiagonalScale_ps_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE PSMatrixAlgebraModule_wrp diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py index 1d09c88b..f653aa6f 100644 --- a/UnitTests/test_psmatrixalgebra.py +++ b/UnitTests/test_psmatrixalgebra.py @@ -284,6 +284,29 @@ def test_reverse(self): self.check_result() + 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.PSMatrix(self.file1) + ntmatrix.DiagonalScale(tlist) + ntmatrix.WriteToMatrixMarket(self.file2) + + ResultMat = mmread(self.file2) + self._compare_mat(CheckMat, ResultMat) + class TestPSMatrixAlgebra_r(TestPSMatrixAlgebra, unittest.TestCase): '''Special routines for real algebra''' @@ -291,6 +314,8 @@ class TestPSMatrixAlgebra_r(TestPSMatrixAlgebra, unittest.TestCase): complex1 = False # Whether the second matrix is complex or not complex2 = False + TripletList = nt.TripletList_r + Triplet = nt.Triplet_r def test_dot(self): '''Test routines to add together matrices.''' @@ -327,6 +352,8 @@ class TestPSMatrixAlgebra_c(TestPSMatrixAlgebra, unittest.TestCase): complex1 = True # Whether the second matrix is complex or not complex2 = True + TripletList = nt.TripletList_c + Triplet = nt.Triplet_c def test_dot(self): '''Test routines to add together matrices.''' @@ -366,6 +393,8 @@ class TestPSMatrixAlgebra_rc(TestPSMatrixAlgebra, unittest.TestCase): complex1 = True # Whether the second matrix is complex or not complex2 = False + TripletList = nt.TripletList_r + Triplet = nt.Triplet_r def setUp(self): '''Set up a specific test.''' @@ -381,6 +410,8 @@ class TestPSMatrixAlgebra_cr(TestPSMatrixAlgebra, unittest.TestCase): complex1 = False # Whether the second matrix is complex or not complex2 = True + TripletList = nt.TripletList_c + Triplet = nt.Triplet_c def setUp(self): '''Set up a specific test.''' From 2e13731515d954fb7deee14d69a8164f33f2e208 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 08:50:23 +0900 Subject: [PATCH 05/16] Arg fix --- Source/CPlusPlus/SMatrix.cc | 2 +- Source/CPlusPlus/SMatrix.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/CPlusPlus/SMatrix.cc b/Source/CPlusPlus/SMatrix.cc index 588677f6..4feda63b 100644 --- a/Source/CPlusPlus/SMatrix.cc +++ b/Source/CPlusPlus/SMatrix.cc @@ -181,7 +181,7 @@ void Matrix_lsr::DiagonalScale(const TripletList_r &tlist) { } //////////////////////////////////////////////////////////////////////////////// -void Matrix_lsc::DiagonalScale(const TripletList_c &tlist, double threshold) { +void Matrix_lsc::DiagonalScale(const TripletList_c &tlist) { MatrixDiagonalScale_lsc_wrp(ih_this, tlist.ih_this); } diff --git a/Source/CPlusPlus/SMatrix.h b/Source/CPlusPlus/SMatrix.h index 21f90116..0accdb88 100644 --- a/Source/CPlusPlus/SMatrix.h +++ b/Source/CPlusPlus/SMatrix.h @@ -177,7 +177,7 @@ class Matrix_lsc { //! 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); + void DiagonalScale(const NTPoly::TripletList_c &tlist); public: //! Transpose a sparse matrix. From 1a99c113f1f4a35788965b825df3dc01dbb3211c Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 09:02:44 +0900 Subject: [PATCH 06/16] Test setup, now get to fixing it --- UnitTests/test_psmatrixalgebra.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py index f653aa6f..5be607f7 100644 --- a/UnitTests/test_psmatrixalgebra.py +++ b/UnitTests/test_psmatrixalgebra.py @@ -288,8 +288,8 @@ 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) + matrix = param.create_matrix(complex=self.complex1) + self.write_matrix(matrix, self.input_file1) CheckMat = deepcopy(matrix) tlist = self.TripletList() @@ -300,9 +300,9 @@ def test_scalediag(self): t.point_value = i tlist.Append(t) CheckMat[:, i] *= i - ntmatrix = self.PSMatrix(self.file1) + ntmatrix = nt.Matrix_ps(self.input_file1) ntmatrix.DiagonalScale(tlist) - ntmatrix.WriteToMatrixMarket(self.file2) + ntmatrix.WriteToMatrixMarket(self.input_file2) ResultMat = mmread(self.file2) self._compare_mat(CheckMat, ResultMat) From 4a9ed1192f89e0dc914091a9539a136d1c645e2b Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 10:12:13 +0900 Subject: [PATCH 07/16] More fixing for the test --- UnitTests/test_psmatrixalgebra.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py index 5be607f7..fea5551c 100644 --- a/UnitTests/test_psmatrixalgebra.py +++ b/UnitTests/test_psmatrixalgebra.py @@ -290,7 +290,7 @@ def test_scalediag(self): for param in self.parameters: matrix = param.create_matrix(complex=self.complex1) self.write_matrix(matrix, self.input_file1) - CheckMat = deepcopy(matrix) + self.CheckMat = deepcopy(matrix) tlist = self.TripletList() for i in range(matrix.shape[1]): @@ -299,13 +299,13 @@ def test_scalediag(self): t.index_row = i + 1 t.point_value = i tlist.Append(t) - CheckMat[:, i] *= i + self.CheckMat[:, i] *= i ntmatrix = nt.Matrix_ps(self.input_file1) ntmatrix.DiagonalScale(tlist) - ntmatrix.WriteToMatrixMarket(self.input_file2) + ntmatrix.WriteToMatrixMarket(self.result_file) - ResultMat = mmread(self.file2) - self._compare_mat(CheckMat, ResultMat) + comm.barrier() + self.check_result() class TestPSMatrixAlgebra_r(TestPSMatrixAlgebra, unittest.TestCase): From d02dc11beed8978f2239bf8a022f3f8007caa4b2 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 10:15:32 +0900 Subject: [PATCH 08/16] Fix partitioning --- .../distributed_algebra_includes/ScaleDiagonal.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 b/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 index ee3b4db5..32492c14 100644 --- a/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 +++ b/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 @@ -1,4 +1,4 @@ -INTEGER :: II, row +INTEGER :: II, col !! Merge to the local block CALL MergeMatrixLocalBlocks(this, lmat) @@ -7,10 +7,10 @@ CALL ConstructTripletList(filtered) DO II = 1, tlist%CurrentSize CALL GetTripletAt(tlist, II, trip) - row = trip%index_row - IF (row .GE. this%start_row .AND. row .LT. this%end_row) THEN - trip%index_row = trip%index_row - this%start_row + 1 - trip%index_column = trip%index_row + col = trip%index_column + IF (col .GE. this%start_column .AND. col .LT. this%end_column) THEN + trip%index_column = trip%index_column - this%start_column + 1 + trip%index_row = trip%index_column CALL AppendToTripletList(filtered, trip) END IF END DO @@ -21,4 +21,4 @@ !! Split CALL SplitMatrixToLocalBlocks(this, lmat) CALL DestructMatrix(lmat) -CALL DestructTripletList(filtered) \ No newline at end of file +CALL DestructTripletList(filtered) From 5053a6ae23fc059e320f81a157f5f4d31bb71d60 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 11:30:13 +0900 Subject: [PATCH 09/16] Gather strategy seems to work --- Source/Fortran/FermiOperatorModule.F90 | 34 ++++++++++++-------- Source/Fortran/TripletListModule.F90 | 43 ++++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 13 deletions(-) diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 index b43d6f0b..768478a4 100644 --- a/Source/Fortran/FermiOperatorModule.F90 +++ b/Source/Fortran/FermiOperatorModule.F90 @@ -7,17 +7,19 @@ MODULE FermiOperatorModule USE LoggingModule, ONLY : WriteElement, WriteHeader, & & EnterSubLog, ExitSubLog, WriteListElement, WriteComment USE PSMatrixAlgebraModule, ONLY : MatrixMultiply, SimilarityTransform, & - & IncrementMatrix, MatrixNorm, ScaleMatrix, DotMatrix, MatrixTrace + & IncrementMatrix, MatrixNorm, ScaleMatrix, DotMatrix, MatrixTrace, & + & MatrixDiagonalScale USE PSMatrixModule, ONLY : Matrix_ps, ConstructEmptyMatrix, & & FillMatrixFromTripletList, GetMatrixTripletList, & - & TransposeMatrix, ConjugateMatrix, DestructMatrix, & - & FillMatrixIdentity, PrintMatrixInformation, CopyMatrix, GetMatrixSize + & TransposeMatrix, ConjugateMatrix, DestructMatrix, FilterMatrix, & + & FillMatrixIdentity, PrintMatrixInformation, CopyMatrix, GetMatrixSize, printmatrix USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, & & DestructMatrixMemoryPool USE SolverParametersModule, ONLY : SolverParameters_t, & & PrintParameters, DestructSolverParameters, CopySolverParameters, & & ConstructSolverParameters - USE TripletListModule, ONLY : TripletList_r, DestructTripletList + USE TripletListModule, ONLY : TripletList_r, DestructTripletList, & + & AllGatherTripletList USE NTMPIModule IMPLICIT NONE PRIVATE @@ -54,7 +56,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & TYPE(Matrix_ps) :: WD TYPE(Matrix_ps) :: vecs, vecsT, vals, Temp TYPE(MatrixMemoryPool_p) :: pool - TYPE(TripletList_r) :: tlist + TYPE(TripletList_r) :: tlist, gathered_list REAL(NTREAL) :: chemical_potential, energy_value REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: eigs, occ REAL(NTREAL) :: sval, sv, occ_temp @@ -148,7 +150,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & CALL ExitSubLog END IF - !! Map + !! Map - note that we store the square root of the occupation numbers energy_value = 0.0_NTREAL DO II = 1, tlist%CurrentSize IF (.NOT. do_smearing) THEN @@ -159,7 +161,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & occ_temp = trace - FLOOR(trace) energy_value = energy_value + & & occ_temp * tlist%DATA(II)%point_value - tlist%DATA(II)%point_value = occ_temp + tlist%DATA(II)%point_value = SQRT(occ_temp) ELSE tlist%DATA(II)%point_value = 0.0_NTREAL ENDIF @@ -168,21 +170,26 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & ! occ_temp = 0.5_NTREAL * (1.0_NTREAL - ERF(inv_temp * sval)) occ_temp = 1.0_NTREAL / (1.0_NTREAL + EXP(inv_temp * sval)) energy_value = energy_value + occ_temp * tlist%DATA(II)%point_value - tlist%DATA(II)%point_value = occ_temp + IF (occ_temp .LT. 0) THEN ! for safety + tlist%DATA(II)%point_value = 0 + ELSE + tlist%DATA(II)%point_value = SQRT(occ_temp) + END IF END IF END DO CALL MPI_ALLREDUCE(MPI_IN_PLACE, energy_value, 1, MPINTREAL, MPI_SUM, & & H%process_grid%within_slice_comm, ierr) + + !! Gather the triplet lists + CALL AllGatherTripletList(tlist, H%process_grid%column_comm, gathered_list) - !! Fill - CALL ConstructEmptyMatrix(vals, H) - CALL FillMatrixFromTripletList(vals, tlist, preduplicated_in = .TRUE.) + !! Scale the eigenvectors + CALL MatrixDiagonalScale(vecs, gathered_list) !! Multiply Back Together - CALL MatrixMultiply(vecs, vals, temp, threshold_in = params%threshold) CALL TransposeMatrix(vecs, vecsT) CALL ConjugateMatrix(vecsT) - CALL MatrixMultiply(temp, vecsT, WD, & + CALL MatrixMultiply(vecs, vecsT, WD, & & threshold_in = params%threshold) !! Compute the density matrix in the non-orthogonalized basis @@ -208,6 +215,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & CALL DestructMatrix(vals) CALL DestructMatrix(temp) CALL DestructTripletList(tlist) + CALL DestructTripletList(gathered_list) CALL DestructMatrixMemoryPool(pool) IF (ALLOCATED(occ)) THEN DEALLOCATE(occ) diff --git a/Source/Fortran/TripletListModule.F90 b/Source/Fortran/TripletListModule.F90 index 29c0de68..e9a19548 100644 --- a/Source/Fortran/TripletListModule.F90 +++ b/Source/Fortran/TripletListModule.F90 @@ -39,6 +39,7 @@ MODULE TripletListModule PUBLIC :: SymmetrizeTripletList PUBLIC :: GetTripletListSize PUBLIC :: RedistributeTripletLists + PUBLIC :: AllGatherTripletList PUBLIC :: ShiftTripletList PUBLIC :: ConvertTripletListType PUBLIC :: MergeTripletLists @@ -91,6 +92,10 @@ MODULE TripletListModule MODULE PROCEDURE RedistributeTripletLists_r MODULE PROCEDURE RedistributeTripletLists_c END INTERFACE RedistributeTripletLists + INTERFACE AllGatherTripletList + MODULE PROCEDURE AllGatherTripletList_r + MODULE PROCEDURE AllGatherTripletList_c + END INTERFACE AllGatherTripletList INTERFACE ShiftTripletList MODULE PROCEDURE ShiftTripletList_r MODULE PROCEDURE ShiftTripletList_c @@ -365,6 +370,44 @@ SUBROUTINE RedistributeTripletLists_c(triplet_lists, comm, local_data_out) #undef MPIDATATYPE END SUBROUTINE RedistributeTripletLists_c +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Gather triplet lists from a set of processors. + SUBROUTINE AllGatherTripletList_r(triplet_in, comm, gathered_out) + !> Locally held triplet list + TYPE(TripletList_r), INTENT(IN) :: triplet_in + !> The mpi communicator to gather along. + INTEGER, INTENT(IN) :: comm + !> The resulting gathered triplet list. + TYPE(TripletList_r), INTENT(INOUT) :: gathered_out + !! Local data (type specific) + REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: send_buffer_val + REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: recv_buffer_val + TYPE(Triplet_r) :: temp_triplet + +#define MPIDATATYPE MPINTREAL +#include "triplet_includes/GatherTripletList.f90" +#undef MPIDATATYPE + + END SUBROUTINE AllGatherTripletList_r +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Gather triplet lists from a set of processors. + SUBROUTINE AllGatherTripletList_c(triplet_in, comm, gathered_out) + !> Locally held triplet list + TYPE(TripletList_c), INTENT(IN) :: triplet_in + !> The mpi communicator to gather along. + INTEGER, INTENT(IN) :: comm + !> The resulting gathered triplet list. + TYPE(TripletList_c), INTENT(INOUT) :: gathered_out + !! Local data (type specific) + COMPLEX(NTCOMPLEX), DIMENSION(:), ALLOCATABLE :: send_buffer_val + COMPLEX(NTCOMPLEX), DIMENSION(:), ALLOCATABLE :: recv_buffer_val + TYPE(Triplet_c) :: temp_triplet + +#define MPIDATATYPE MPINTCOMPLEX +#include "triplet_includes/GatherTripletList.f90" +#undef MPIDATATYPE + + END SUBROUTINE AllGatherTripletList_c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Shift the rows and columns of a triplet list by set values. !> Frequently, we have a triplet list that comes from the global matrix which From 5c554bd7bc7b03d62fb566be5c1286b678598d95 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 11:30:32 +0900 Subject: [PATCH 10/16] forgot a file --- .../triplet_includes/GatherTripletList.f90 | 66 +++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 Source/Fortran/triplet_includes/GatherTripletList.f90 diff --git a/Source/Fortran/triplet_includes/GatherTripletList.f90 b/Source/Fortran/triplet_includes/GatherTripletList.f90 new file mode 100644 index 00000000..2c19518f --- /dev/null +++ b/Source/Fortran/triplet_includes/GatherTripletList.f90 @@ -0,0 +1,66 @@ +!! Local Data - Send/Recv Buffers +INTEGER, DIMENSION(:), ALLOCATABLE :: send_buffer_row +INTEGER, DIMENSION(:), ALLOCATABLE :: send_buffer_col +INTEGER, DIMENSION(:), ALLOCATABLE :: recv_buffer_row +INTEGER, DIMENSION(:), ALLOCATABLE :: recv_buffer_col + +!! Sizes help +INTEGER, DIMENSION(:), ALLOCATABLE :: recvcounts +INTEGER, DIMENSION(:), ALLOCATABLE :: displ +INTEGER :: gather_size + +!! Temporary variables +INTEGER :: num_processes, II, ierr + +!! Figure out the comm size +CALL MPI_COMM_SIZE(comm, num_processes, ierr) + +!! Get the count +ALLOCATE(recvcounts(num_processes)) +CALL MPI_Allgather(triplet_in%CurrentSize, 1, MPI_INTEGER, recvcounts, & + & 1, MPI_INTEGER, comm, ierr) + +!! Get the displacements +gather_size = SUM(recvcounts) +ALLOCATE(displ(num_processes)) +displ(1) = 0 +DO II = 2, num_processes + displ(II) = displ(II - 1) + recvcounts(II - 1) +END DO + +!! Prepare the send buffers +ALLOCATE(send_buffer_row(triplet_in%CurrentSize)) +ALLOCATE(send_buffer_col(triplet_in%CurrentSize)) +ALLOCATE(send_buffer_val(triplet_in%CurrentSize)) +DO II = 1, triplet_in%CurrentSize + CALL GetTripletAt(triplet_in, II, temp_triplet) + send_buffer_row(II) = temp_triplet%index_row + send_buffer_col(II) = temp_triplet%index_column + send_buffer_val(II) = temp_triplet%point_value +END DO + +!! Gather Call +ALLOCATE(recv_buffer_row(gather_size)) +ALLOCATE(recv_buffer_col(gather_size)) +ALLOCATE(recv_buffer_val(gather_size)) +CALL MPI_Allgatherv(send_buffer_row, triplet_in%CurrentSize, MPI_INTEGER, & + & recv_buffer_row, recvcounts, displ, MPI_INTEGER, comm, ierr) +CALL MPI_Allgatherv(send_buffer_col, triplet_in%CurrentSize, MPI_INTEGER, & + & recv_buffer_col, recvcounts, displ, MPI_INTEGER, comm, ierr) +CALL MPI_Allgatherv(send_buffer_val, triplet_in%CurrentSize, MPIDATATYPE, & + & recv_buffer_val, recvcounts, displ, MPIDATATYPE, comm, ierr) + +!! Unpack +CALL ConstructTripletList(gathered_out, gather_size) +DO II = 1, gather_size + gathered_out%DATA(II)%index_row = recv_buffer_row(II) + gathered_out%DATA(II)%index_column = recv_buffer_col(II) + gathered_out%DATA(II)%point_value = recv_buffer_val(II) +END DO + +!! Cleanup +DEALLOCATE(recvcounts) +DEALLOCATE(displ) +DEALLOCATE(send_buffer_row) +DEALLOCATE(send_buffer_col) +DEALLOCATE(send_buffer_val) From b35062cf29f1154dcf6788353c5a739a383e4c9d Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 11:42:06 +0900 Subject: [PATCH 11/16] Seems to be working --- Source/Fortran/FermiOperatorModule.F90 | 1 + Source/Fortran/distributed_includes/FilterMatrix.f90 | 12 ++---------- 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 index 768478a4..b8f88a03 100644 --- a/Source/Fortran/FermiOperatorModule.F90 +++ b/Source/Fortran/FermiOperatorModule.F90 @@ -185,6 +185,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & !! Scale the eigenvectors CALL MatrixDiagonalScale(vecs, gathered_list) + CALL FilterMatrix(vecs, params%threshold) !! Multiply Back Together CALL TransposeMatrix(vecs, vecsT) diff --git a/Source/Fortran/distributed_includes/FilterMatrix.f90 b/Source/Fortran/distributed_includes/FilterMatrix.f90 index 0373f513..2530b267 100644 --- a/Source/Fortran/distributed_includes/FilterMatrix.f90 +++ b/Source/Fortran/distributed_includes/FilterMatrix.f90 @@ -1,7 +1,5 @@ !! Local Data - TYPE(ProcessGrid_t) :: grid_temp - LOGICAL :: is_complex_temp - INTEGER :: II, size_temp + INTEGER :: II CALL GetMatrixTripletList(this, tlist) CALL ConstructTripletList(new_list) @@ -13,10 +11,4 @@ END IF END DO - size_temp = this%actual_matrix_dimension - grid_temp = this%process_grid - is_complex_temp = this%is_complex - CALL DestructMatrix(this) - CALL ConstructEmptyMatrix(this, size_temp, grid_temp, is_complex_temp) - - CALL FillMatrixFromTripletList(this, new_list) + CALL FillMatrixFromTripletList(this, new_list, preduplicated_in=.TRUE.) From 3f0a31a5764c3c4a4f8d7478a5e135b45ddbc963 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 15:26:48 +0900 Subject: [PATCH 12/16] Fix case of failing to upcast --- Source/Fortran/FermiOperatorModule.F90 | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 index b8f88a03..56ab94f5 100644 --- a/Source/Fortran/FermiOperatorModule.F90 +++ b/Source/Fortran/FermiOperatorModule.F90 @@ -1,7 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Computing The Density Matrix Using the Fermi Operator Expansion MODULE FermiOperatorModule - USE DataTypesModule, ONLY : NTREAL, MPINTREAL + USE DataTypesModule, ONLY : NTREAL, MPINTREAL, NTCOMPLEX USE EigenSolversModule, ONLY : EigenDecomposition USE LoadBalancerModule, ONLY : PermuteMatrix, UndoPermuteMatrix USE LoggingModule, ONLY : WriteElement, WriteHeader, & @@ -18,8 +18,8 @@ MODULE FermiOperatorModule USE SolverParametersModule, ONLY : SolverParameters_t, & & PrintParameters, DestructSolverParameters, CopySolverParameters, & & ConstructSolverParameters - USE TripletListModule, ONLY : TripletList_r, DestructTripletList, & - & AllGatherTripletList + USE TripletListModule, ONLY : TripletList_r, TripletList_c, & + & DestructTripletList, AllGatherTripletList, ConstructTripletList USE NTMPIModule IMPLICIT NONE PRIVATE @@ -57,6 +57,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & TYPE(Matrix_ps) :: vecs, vecsT, vals, Temp TYPE(MatrixMemoryPool_p) :: pool TYPE(TripletList_r) :: tlist, gathered_list + TYPE(TripletList_c) :: gathered_list_c REAL(NTREAL) :: chemical_potential, energy_value REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: eigs, occ REAL(NTREAL) :: sval, sv, occ_temp @@ -184,7 +185,19 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & CALL AllGatherTripletList(tlist, H%process_grid%column_comm, gathered_list) !! Scale the eigenvectors - CALL MatrixDiagonalScale(vecs, gathered_list) + IF (vecs%is_complex) THEN + CALL ConstructTripletList(gathered_list_c, gathered_list%CurrentSize) + DO II = 1, gathered_list%CurrentSize + gathered_list_c%DATA(II)%index_row = gathered_list%DATA(II)%index_row + gathered_list_c%DATA(II)%index_column = & + & gathered_list%DATA(II)%index_column + gathered_list_c%DATA(II)%point_value = & + & CMPLX(gathered_list%DATA(II)%point_value, KIND=NTCOMPLEX) + END DO + CALL MatrixDiagonalScale(vecs, gathered_list_c) + ELSE + CALL MatrixDiagonalScale(vecs, gathered_list) + END IF CALL FilterMatrix(vecs, params%threshold) !! Multiply Back Together @@ -217,6 +230,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & CALL DestructMatrix(temp) CALL DestructTripletList(tlist) CALL DestructTripletList(gathered_list) + CALL DestructTripletList(gathered_list_c) CALL DestructMatrixMemoryPool(pool) IF (ALLOCATED(occ)) THEN DEALLOCATE(occ) From 9826c0ac9ed542dbc93255dac1e66bbdd9d50411 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 16:27:43 +0900 Subject: [PATCH 13/16] Fix test conditions to handle the case an edge case --- Source/C/PSMatrix_c.h | 3 ++- Source/CPlusPlus/PSMatrix.cc | 7 +++++- Source/CPlusPlus/PSMatrix.h | 3 +++ Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 | 25 ++++++++++++++++---- UnitTests/test_psmatrixalgebra.py | 12 ++++++---- 5 files changed, 39 insertions(+), 11 deletions(-) diff --git a/Source/C/PSMatrix_c.h b/Source/C/PSMatrix_c.h index c07bd7a0..a211fe73 100644 --- a/Source/C/PSMatrix_c.h +++ b/Source/C/PSMatrix_c.h @@ -62,6 +62,7 @@ void ScaleMatrix_ps_wrp(int *ih_this, const double *constant); double MatrixNorm_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 MatrixDiagonalScale_ps_wrp(int *ih_mat, const int *ih_tlist); +void MatrixDiagonalScale_psr_wrp(int *ih_mat, const int *ih_tlist); +void MatrixDiagonalScale_psc_wrp(int *ih_mat, const int *ih_tlist); #endif diff --git a/Source/CPlusPlus/PSMatrix.cc b/Source/CPlusPlus/PSMatrix.cc index d421c4df..05f0d593 100644 --- a/Source/CPlusPlus/PSMatrix.cc +++ b/Source/CPlusPlus/PSMatrix.cc @@ -219,7 +219,12 @@ void Matrix_ps::Scale(double constant) { ////////////////////////////////////////////////////////////////////////////// void Matrix_ps::DiagonalScale(const TripletList_r &tlist) { - MatrixDiagonalScale_ps_wrp(ih_this, tlist.ih_this); + MatrixDiagonalScale_psr_wrp(ih_this, tlist.ih_this); +} + +////////////////////////////////////////////////////////////////////////////// +void Matrix_ps::DiagonalScale(const TripletList_c &tlist) { + MatrixDiagonalScale_psc_wrp(ih_this, tlist.ih_this); } ////////////////////////////////////////////////////////////////////////////// diff --git a/Source/CPlusPlus/PSMatrix.h b/Source/CPlusPlus/PSMatrix.h index d6e3bc90..8aba171a 100644 --- a/Source/CPlusPlus/PSMatrix.h +++ b/Source/CPlusPlus/PSMatrix.h @@ -170,6 +170,9 @@ class Matrix_ps { //!\param tlist the triplet list. //!\param threshold for flushing small values. void DiagonalScale(const NTPoly::TripletList_r &tlist); + //!\param tlist the triplet list. + //!\param threshold for flushing small values. + void DiagonalScale(const NTPoly::TripletList_c &tlist); public: //! Destructor. diff --git a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 index e333d6ac..f57aa77a 100644 --- a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 +++ b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 @@ -6,7 +6,7 @@ MODULE PSMatrixAlgebraModule_wrp USE PSMatrixAlgebraModule USE PSMatrixModule_wrp, ONLY : Matrix_ps_wrp USE PermutationModule_wrp, ONLY : Permutation_wrp - USE TripletListModule_wrp, ONLY : TripletList_r_wrp + USE TripletListModule_wrp, ONLY : TripletList_r_wrp, TripletList_c_wrp USE WrapperModule, ONLY : SIZE_wrp USE ISO_C_BINDING, ONLY : c_int, c_char, c_bool IMPLICIT NONE @@ -20,7 +20,8 @@ MODULE PSMatrixAlgebraModule_wrp PUBLIC :: ScaleMatrix_ps_wrp PUBLIC :: MatrixNorm_ps_wrp PUBLIC :: MatrixTrace_ps_wrp - PUBLIC :: MatrixDiagonalScale_ps_wrp + PUBLIC :: MatrixDiagonalScale_psr_wrp + PUBLIC :: MatrixDiagonalScale_psc_wrp CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Matrix B = alpha*Matrix A + Matrix B (AXPY) SUBROUTINE IncrementMatrix_ps_wrp(ih_matA, ih_matB, alpha_in,threshold_in) & @@ -146,8 +147,8 @@ SUBROUTINE MatrixTrace_ps_wrp(ih_this, trace_value) & END SUBROUTINE MatrixTrace_ps_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Scale a matrix using a diagonal matrix (triplet list form). - SUBROUTINE MatrixDiagonalScale_ps_wrp(ih_mat, ih_tlist) & - & BIND(c,name="MatrixDiagonalScale_ps_wrp") + SUBROUTINE MatrixDiagonalScale_psr_wrp(ih_mat, ih_tlist) & + & BIND(c,name="MatrixDiagonalScale_psr_wrp") INTEGER(kind=c_int), INTENT(INOUT) :: ih_mat(SIZE_wrp) INTEGER(kind=c_int), INTENT(IN) :: ih_tlist(SIZE_wrp) TYPE(Matrix_ps_wrp) :: h_mat @@ -157,6 +158,20 @@ SUBROUTINE MatrixDiagonalScale_ps_wrp(ih_mat, ih_tlist) & h_tlist = TRANSFER(ih_tlist, h_tlist) CALL MatrixDiagonalScale(h_mat%DATA, h_tlist%DATA) - END SUBROUTINE MatrixDiagonalScale_ps_wrp + END SUBROUTINE MatrixDiagonalScale_psr_wrp +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Scale a matrix using a diagonal matrix (triplet list form). + SUBROUTINE MatrixDiagonalScale_psc_wrp(ih_mat, ih_tlist) & + & BIND(c,name="MatrixDiagonalScale_psc_wrp") + INTEGER(kind=c_int), INTENT(INOUT) :: ih_mat(SIZE_wrp) + INTEGER(kind=c_int), INTENT(IN) :: ih_tlist(SIZE_wrp) + TYPE(Matrix_ps_wrp) :: h_mat + TYPE(TripletList_c_wrp) :: h_tlist + + h_mat = TRANSFER(ih_mat, h_mat) + h_tlist = TRANSFER(ih_tlist, h_tlist) + + CALL MatrixDiagonalScale(h_mat%DATA, h_tlist%DATA) + END SUBROUTINE MatrixDiagonalScale_psc_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE PSMatrixAlgebraModule_wrp diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py index fea5551c..052dff8c 100644 --- a/UnitTests/test_psmatrixalgebra.py +++ b/UnitTests/test_psmatrixalgebra.py @@ -292,6 +292,10 @@ def test_scalediag(self): self.write_matrix(matrix, self.input_file1) self.CheckMat = deepcopy(matrix) + # Need a guard because otherwise we write a real matrix by mistake. + if self.complex1 and matrix.nnz == 0: + continue + tlist = self.TripletList() for i in range(matrix.shape[1]): t = self.Triplet() @@ -393,8 +397,8 @@ class TestPSMatrixAlgebra_rc(TestPSMatrixAlgebra, unittest.TestCase): complex1 = True # Whether the second matrix is complex or not complex2 = False - TripletList = nt.TripletList_r - Triplet = nt.Triplet_r + TripletList = nt.TripletList_c + Triplet = nt.Triplet_c def setUp(self): '''Set up a specific test.''' @@ -410,8 +414,8 @@ class TestPSMatrixAlgebra_cr(TestPSMatrixAlgebra, unittest.TestCase): complex1 = False # Whether the second matrix is complex or not complex2 = True - TripletList = nt.TripletList_c - Triplet = nt.Triplet_c + TripletList = nt.TripletList_r + Triplet = nt.Triplet_r def setUp(self): '''Set up a specific test.''' From 306286782a3269b7572b2faaf807b91a8e204db7 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 17:28:51 +0900 Subject: [PATCH 14/16] Try to simplify the code --- Source/Fortran/FermiOperatorModule.F90 | 44 +++++++------------ Source/Fortran/PSMatrixModule.F90 | 29 ++++++++++++ Source/Fortran/TripletListModule.F90 | 21 +++++++++ .../GatherMatrixTripletList.f90 | 3 ++ 4 files changed, 69 insertions(+), 28 deletions(-) create mode 100644 Source/Fortran/distributed_includes/GatherMatrixTripletList.f90 diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 index 56ab94f5..5b4e1ddb 100644 --- a/Source/Fortran/FermiOperatorModule.F90 +++ b/Source/Fortran/FermiOperatorModule.F90 @@ -12,14 +12,15 @@ MODULE FermiOperatorModule USE PSMatrixModule, ONLY : Matrix_ps, ConstructEmptyMatrix, & & FillMatrixFromTripletList, GetMatrixTripletList, & & TransposeMatrix, ConjugateMatrix, DestructMatrix, FilterMatrix, & - & FillMatrixIdentity, PrintMatrixInformation, CopyMatrix, GetMatrixSize, printmatrix + & FillMatrixIdentity, PrintMatrixInformation, CopyMatrix, & + & GatherMatrixTripletList, GetMatrixSize USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, & & DestructMatrixMemoryPool USE SolverParametersModule, ONLY : SolverParameters_t, & & PrintParameters, DestructSolverParameters, CopySolverParameters, & & ConstructSolverParameters USE TripletListModule, ONLY : TripletList_r, TripletList_c, & - & DestructTripletList, AllGatherTripletList, ConstructTripletList + & DestructTripletList, CopyTripletList USE NTMPIModule IMPLICIT NONE PRIVATE @@ -56,15 +57,14 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & TYPE(Matrix_ps) :: WD TYPE(Matrix_ps) :: vecs, vecsT, vals, Temp TYPE(MatrixMemoryPool_p) :: pool - TYPE(TripletList_r) :: tlist, gathered_list - TYPE(TripletList_c) :: gathered_list_c + TYPE(TripletList_r) :: tlist + TYPE(TripletList_c) :: tlist_c REAL(NTREAL) :: chemical_potential, energy_value REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: eigs, occ REAL(NTREAL) :: sval, sv, occ_temp REAL(NTREAL) :: left, right, homo, lumo INTEGER :: num_eigs INTEGER :: II, JJ - INTEGER :: ierr !! Optional Parameters IF (PRESENT(solver_parameters_in)) THEN @@ -103,15 +103,15 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & & eigenvectors_in = vecs, solver_parameters_in = params) !! Gather the eigenvalues on to every process - CALL GetMatrixTripletList(vals, tlist) + CALL GatherMatrixTripletList(vals, tlist) + + !! Put them in an array for simplicity num_eigs = H%actual_matrix_dimension ALLOCATE(eigs(num_eigs)) eigs = 0 DO II = 1, tlist%CurrentSize - eigs(tlist%DATA(II)%index_column) = tlist%DATA(II)%point_value + eigs(II) = tlist%DATA(II)%point_value END DO - CALL MPI_ALLREDUCE(MPI_IN_PLACE, eigs, num_eigs, MPINTREAL, & - & MPI_SUM, H%process_grid%within_slice_comm, ierr) !! Compute MU By Bisection IF (do_smearing) THEN @@ -119,7 +119,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & left = MINVAL(eigs) right = MAXVAL(eigs) DO JJ = 1, 10*params%max_iterations - chemical_potential = left + (right - left)/2 + chemical_potential = left + (right - left) / 2 DO II = 1, num_eigs sval = eigs(II) - chemical_potential ! occ(II) = 0.5_NTREAL * (1.0_NTREAL - ERF(inv_temp * sval)) @@ -170,7 +170,8 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & sval = tlist%DATA(II)%point_value - chemical_potential ! occ_temp = 0.5_NTREAL * (1.0_NTREAL - ERF(inv_temp * sval)) occ_temp = 1.0_NTREAL / (1.0_NTREAL + EXP(inv_temp * sval)) - energy_value = energy_value + occ_temp * tlist%DATA(II)%point_value + energy_value = energy_value + & + & occ_temp * tlist%DATA(II)%point_value IF (occ_temp .LT. 0) THEN ! for safety tlist%DATA(II)%point_value = 0 ELSE @@ -178,25 +179,13 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & END IF END IF END DO - CALL MPI_ALLREDUCE(MPI_IN_PLACE, energy_value, 1, MPINTREAL, MPI_SUM, & - & H%process_grid%within_slice_comm, ierr) - - !! Gather the triplet lists - CALL AllGatherTripletList(tlist, H%process_grid%column_comm, gathered_list) !! Scale the eigenvectors IF (vecs%is_complex) THEN - CALL ConstructTripletList(gathered_list_c, gathered_list%CurrentSize) - DO II = 1, gathered_list%CurrentSize - gathered_list_c%DATA(II)%index_row = gathered_list%DATA(II)%index_row - gathered_list_c%DATA(II)%index_column = & - & gathered_list%DATA(II)%index_column - gathered_list_c%DATA(II)%point_value = & - & CMPLX(gathered_list%DATA(II)%point_value, KIND=NTCOMPLEX) - END DO - CALL MatrixDiagonalScale(vecs, gathered_list_c) + CALL CopyTripletList(tlist, tlist_c) + CALL MatrixDiagonalScale(vecs, tlist_c) ELSE - CALL MatrixDiagonalScale(vecs, gathered_list) + CALL MatrixDiagonalScale(vecs, tlist) END IF CALL FilterMatrix(vecs, params%threshold) @@ -229,8 +218,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & CALL DestructMatrix(vals) CALL DestructMatrix(temp) CALL DestructTripletList(tlist) - CALL DestructTripletList(gathered_list) - CALL DestructTripletList(gathered_list_c) + CALL DestructTripletList(tlist_c) CALL DestructMatrixMemoryPool(pool) IF (ALLOCATED(occ)) THEN DEALLOCATE(occ) diff --git a/Source/Fortran/PSMatrixModule.F90 b/Source/Fortran/PSMatrixModule.F90 index 423254e8..75f60b1e 100644 --- a/Source/Fortran/PSMatrixModule.F90 +++ b/Source/Fortran/PSMatrixModule.F90 @@ -87,6 +87,7 @@ MODULE PSMatrixModule PUBLIC :: CommSplitMatrix PUBLIC :: ResizeMatrix PUBLIC :: GatherMatrixToProcess + PUBLIC :: GatherMatrixTripletList PUBLIC :: IsIdentity !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INTERFACE ConstructEmptyMatrix @@ -180,6 +181,10 @@ MODULE PSMatrixModule MODULE PROCEDURE GatherMatrixToProcess_psc_id MODULE PROCEDURE GatherMatrixToProcess_psc_all END INTERFACE GatherMatrixToProcess + INTERFACE GatherMatrixTripletList + MODULE PROCEDURE GatherMatrixTripletList_r + MODULE PROCEDURE GatherMatrixTripletList_c + END INTERFACE GatherMatrixTripletList CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Construct an empty sparse, distributed, matrix. SUBROUTINE ConstructEmptyMatrix_ps(this, matrix_dim, process_grid_in, & @@ -1843,5 +1848,29 @@ FUNCTION IsIdentity_psc(this) RESULT(is_identity) #include "distributed_includes/IsIdentity.f90" #undef ISCOMPLEX END FUNCTION IsIdentity_psc +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> This gathers the entire matrix into a triplet list on all processes. + SUBROUTINE GatherMatrixTripletList_r(this, tlist) + !> The matrix to gather. + TYPE(Matrix_ps), INTENT(IN) :: this + !> The full matrix, stored in a local matrix. + TYPE(TripletList_r), INTENT(INOUT) :: tlist + !! Local Variables + TYPE(Matrix_lsr) :: lmat + +#include "distributed_includes/GatherMatrixTripletList.f90" + END SUBROUTINE GatherMatrixTripletList_r +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> This gathers the entire matrix into a triplet list on all processes. + SUBROUTINE GatherMatrixTripletList_c(this, tlist) + !> The matrix to gather. + TYPE(Matrix_ps), INTENT(IN) :: this + !> The full matrix, stored in a local matrix. + TYPE(TripletList_c), INTENT(INOUT) :: tlist + !! Local Variables + TYPE(Matrix_lsc) :: lmat + +#include "distributed_includes/GatherMatrixTripletList.f90" + END SUBROUTINE GatherMatrixTripletList_c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE PSMatrixModule diff --git a/Source/Fortran/TripletListModule.F90 b/Source/Fortran/TripletListModule.F90 index e9a19548..c74f7ee3 100644 --- a/Source/Fortran/TripletListModule.F90 +++ b/Source/Fortran/TripletListModule.F90 @@ -51,6 +51,7 @@ MODULE TripletListModule INTERFACE CopyTripletList MODULE PROCEDURE CopyTripletList_r MODULE PROCEDURE CopyTripletList_c + MODULE PROCEDURE CopyTripletList_rc END INTERFACE CopyTripletList INTERFACE DestructTripletList MODULE PROCEDURE DestructTripletList_r @@ -149,6 +150,7 @@ PURE SUBROUTINE DestructTripletList_c(this) END SUBROUTINE DestructTripletList_c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Copy a triplet list (real). SUBROUTINE CopyTripletList_r(tripA, tripB) !> The triplet list to copy. TYPE(TripletList_r), INTENT(IN) :: tripA @@ -158,6 +160,7 @@ SUBROUTINE CopyTripletList_r(tripA, tripB) #include "triplet_includes/CopyTripletList.f90" END SUBROUTINE CopyTripletList_r !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Copy a triplet list (complex). SUBROUTINE CopyTripletList_c(tripA, tripB) !> The triplet list to copy. TYPE(TripletList_c), INTENT(IN) :: tripA @@ -166,6 +169,24 @@ SUBROUTINE CopyTripletList_c(tripA, tripB) #include "triplet_includes/CopyTripletList.f90" END SUBROUTINE CopyTripletList_c +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Copy and upcast a triplet list (real -> complex). + SUBROUTINE CopyTripletList_rc(tripA, tripB) + !> The triplet list to copy. + TYPE(TripletList_r), INTENT(IN) :: tripA + !> tripB = tripA + TYPE(TripletList_c), INTENT(INOUT) :: tripB + !! Local varaibles + INTEGER II + + CALL ConstructTripletList(tripB, tripA%CurrentSize) + DO II = 1, tripA%CurrentSize + tripB%DATA(II)%index_row = tripA%DATA(II)%index_row + tripB%DATA(II)%index_column = tripA%DATA(II)%index_column + tripB%DATA(II)%point_value = & + & CMPLX(tripA%DATA(II)%point_value, KIND=NTCOMPLEX) + END DO + END SUBROUTINE CopyTripletList_rc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Increase the size of a triplet list. PURE SUBROUTINE ResizeTripletList_r(this, size) diff --git a/Source/Fortran/distributed_includes/GatherMatrixTripletList.f90 b/Source/Fortran/distributed_includes/GatherMatrixTripletList.f90 new file mode 100644 index 00000000..f4fa0949 --- /dev/null +++ b/Source/Fortran/distributed_includes/GatherMatrixTripletList.f90 @@ -0,0 +1,3 @@ +CALL GatherMatrixToProcess(this, lmat) +CALL MatrixToTripletList(lmat, tlist) +CALL DestructMatrix(lmat) \ No newline at end of file From f5b45ccefe0f79119532208d6d85533426da71d7 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 17:40:45 +0900 Subject: [PATCH 15/16] Move the approach into the eigensolvers module --- Source/Fortran/EigenSolversModule.F90 | 28 +++++++++++++++----------- Source/Fortran/FermiOperatorModule.F90 | 2 +- 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/Source/Fortran/EigenSolversModule.F90 b/Source/Fortran/EigenSolversModule.F90 index 103e751c..fc820627 100644 --- a/Source/Fortran/EigenSolversModule.F90 +++ b/Source/Fortran/EigenSolversModule.F90 @@ -13,10 +13,10 @@ MODULE EigenSolversModule USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, & & DestructMatrixMemoryPool USE PSMatrixAlgebraModule, ONLY : IncrementMatrix, MatrixMultiply, & - & ScaleMatrix + & ScaleMatrix, MatrixDiagonalScale USE PSMatrixModule, ONLY : Matrix_ps, GatherMatrixToProcess, & & FillMatrixFromTripletList, ConstructEmptyMatrix, ConvertMatrixToReal, & - & DestructMatrix, CopyMatrix, GetMatrixTripletList, TransposeMatrix, & + & DestructMatrix, CopyMatrix, GatherMatrixTripletList, TransposeMatrix, & & ConjugateMatrix, FillMatrixIdentity, WriteMatrixToMatrixMarket USE SolverParametersModule, ONLY : SolverParameters_t, PrintParameters, & & DestructSolverParameters, ConstructSolverParameters, & @@ -24,7 +24,7 @@ MODULE EigenSolversModule USE SMatrixModule, ONLY : Matrix_lsr, Matrix_lsc, MatrixToTripletList, & & DestructMatrix USE TripletListModule, ONLY : TripletList_r, TripletList_c, & - & ConstructTripletList, DestructTripletList + & ConstructTripletList, CopyTripletList, DestructTripletList USE NTMPIModule IMPLICIT NONE PRIVATE @@ -107,6 +107,7 @@ END FUNCTION func !! Local Variables TYPE(Matrix_ps) :: vecs, vecsT, vals, temp TYPE(TripletList_r) :: tlist + TYPE(TripletList_c) :: tlist_c INTEGER :: II !! Optional Parameters @@ -121,26 +122,29 @@ END FUNCTION func & eigenvectors_in = vecs) !! Convert to a triplet list, map the triplet list, fill. - CALL GetMatrixTripletList(vals, tlist) + CALL GatherMatrixTripletList(vals, tlist) CALL DestructMatrix(vals) DO II = 1, tlist%CurrentSize tlist%DATA(II)%point_value = func(tlist%DATA(II)%point_value) END DO - !! Fill - CALL ConstructEmptyMatrix(ResultMat, this) - CALL FillMatrixFromTripletList(ResultMat, tlist, preduplicated_in = .TRUE.) - CALL DestructTripletList(tlist) - !! Multiply Back Together - CALL MatrixMultiply(vecs, ResultMat, temp, threshold_in = params%threshold) CALL TransposeMatrix(vecs, vecsT) - CALL DestructMatrix(vecs) CALL ConjugateMatrix(vecsT) - CALL MatrixMultiply(temp, vecsT, ResultMat, threshold_in = params%threshold) + IF (this%is_complex) THEN + CALL CopyTripletList(tlist, tlist_c) + CALL MatrixDiagonalScale(vecs, tlist_c) + ELSE + CALL MatrixDiagonalScale(vecs, tlist) + END IF + CALL MatrixMultiply(vecs, vecsT, ResultMat, threshold_in = params%threshold) !! Cleanup + CALL DestructMatrix(vecs) + CALL DestructMatrix(vecsT) CALL DestructMatrix(temp) + CALL DestructTripletList(tlist) + CALL DestructTripletList(tlist_c) CALL DestructSolverParameters(params) END SUBROUTINE DenseMatrixFunction diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 index 5b4e1ddb..9a4a2049 100644 --- a/Source/Fortran/FermiOperatorModule.F90 +++ b/Source/Fortran/FermiOperatorModule.F90 @@ -182,7 +182,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & !! Scale the eigenvectors IF (vecs%is_complex) THEN - CALL CopyTripletList(tlist, tlist_c) + CALL CopyTripletList(tlist, tlist_c) CALL MatrixDiagonalScale(vecs, tlist_c) ELSE CALL MatrixDiagonalScale(vecs, tlist) From 1620c8e942e89b0e68df1a278d6dd8e7c3dcc470 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Thu, 18 Apr 2024 17:56:52 +0900 Subject: [PATCH 16/16] lint --- Source/Fortran/FermiOperatorModule.F90 | 2 +- Source/Fortran/TripletListModule.F90 | 2 +- .../ScaleDiagonal.f90 | 34 +++--- .../GatherMatrixTripletList.f90 | 6 +- .../Fortran/sparse_includes/DiagonalScale.f90 | 16 +-- .../triplet_includes/GatherTripletList.f90 | 114 +++++++++--------- 6 files changed, 87 insertions(+), 87 deletions(-) diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 index 9a4a2049..fafeb1c3 100644 --- a/Source/Fortran/FermiOperatorModule.F90 +++ b/Source/Fortran/FermiOperatorModule.F90 @@ -171,7 +171,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & ! occ_temp = 0.5_NTREAL * (1.0_NTREAL - ERF(inv_temp * sval)) occ_temp = 1.0_NTREAL / (1.0_NTREAL + EXP(inv_temp * sval)) energy_value = energy_value + & - & occ_temp * tlist%DATA(II)%point_value + & occ_temp * tlist%DATA(II)%point_value IF (occ_temp .LT. 0) THEN ! for safety tlist%DATA(II)%point_value = 0 ELSE diff --git a/Source/Fortran/TripletListModule.F90 b/Source/Fortran/TripletListModule.F90 index c74f7ee3..5c967d2a 100644 --- a/Source/Fortran/TripletListModule.F90 +++ b/Source/Fortran/TripletListModule.F90 @@ -184,7 +184,7 @@ SUBROUTINE CopyTripletList_rc(tripA, tripB) tripB%DATA(II)%index_row = tripA%DATA(II)%index_row tripB%DATA(II)%index_column = tripA%DATA(II)%index_column tripB%DATA(II)%point_value = & - & CMPLX(tripA%DATA(II)%point_value, KIND=NTCOMPLEX) + & CMPLX(tripA%DATA(II)%point_value, KIND=NTCOMPLEX) END DO END SUBROUTINE CopyTripletList_rc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 b/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 index 32492c14..01f7c66e 100644 --- a/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 +++ b/Source/Fortran/distributed_algebra_includes/ScaleDiagonal.f90 @@ -1,24 +1,24 @@ -INTEGER :: II, col + INTEGER :: II, col -!! Merge to the local block -CALL MergeMatrixLocalBlocks(this, lmat) + !! Merge to the local block + CALL MergeMatrixLocalBlocks(this, lmat) -!! Filter out the triplets that aren't stored locally -CALL ConstructTripletList(filtered) -DO II = 1, tlist%CurrentSize - CALL GetTripletAt(tlist, II, trip) - col = trip%index_column - IF (col .GE. this%start_column .AND. col .LT. this%end_column) THEN + !! Filter out the triplets that aren't stored locally + CALL ConstructTripletList(filtered) + DO II = 1, tlist%CurrentSize + CALL GetTripletAt(tlist, II, trip) + col = trip%index_column + IF (col .GE. this%start_column .AND. col .LT. this%end_column) THEN trip%index_column = trip%index_column - this%start_column + 1 trip%index_row = trip%index_column CALL AppendToTripletList(filtered, trip) - END IF -END DO + END IF + END DO -!! Scale -CALL MatrixDiagonalScale(lmat, filtered) + !! Scale + CALL MatrixDiagonalScale(lmat, filtered) -!! Split -CALL SplitMatrixToLocalBlocks(this, lmat) -CALL DestructMatrix(lmat) -CALL DestructTripletList(filtered) + !! Split + CALL SplitMatrixToLocalBlocks(this, lmat) + CALL DestructMatrix(lmat) + CALL DestructTripletList(filtered) diff --git a/Source/Fortran/distributed_includes/GatherMatrixTripletList.f90 b/Source/Fortran/distributed_includes/GatherMatrixTripletList.f90 index f4fa0949..ecb280ff 100644 --- a/Source/Fortran/distributed_includes/GatherMatrixTripletList.f90 +++ b/Source/Fortran/distributed_includes/GatherMatrixTripletList.f90 @@ -1,3 +1,3 @@ -CALL GatherMatrixToProcess(this, lmat) -CALL MatrixToTripletList(lmat, tlist) -CALL DestructMatrix(lmat) \ No newline at end of file + CALL GatherMatrixToProcess(this, lmat) + CALL MatrixToTripletList(lmat, tlist) + CALL DestructMatrix(lmat) diff --git a/Source/Fortran/sparse_includes/DiagonalScale.f90 b/Source/Fortran/sparse_includes/DiagonalScale.f90 index c1997b5b..2bb20fae 100644 --- a/Source/Fortran/sparse_includes/DiagonalScale.f90 +++ b/Source/Fortran/sparse_includes/DiagonalScale.f90 @@ -1,9 +1,9 @@ + + INTEGER :: col, II -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 + 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/Fortran/triplet_includes/GatherTripletList.f90 b/Source/Fortran/triplet_includes/GatherTripletList.f90 index 2c19518f..9b1587df 100644 --- a/Source/Fortran/triplet_includes/GatherTripletList.f90 +++ b/Source/Fortran/triplet_includes/GatherTripletList.f90 @@ -1,66 +1,66 @@ -!! Local Data - Send/Recv Buffers -INTEGER, DIMENSION(:), ALLOCATABLE :: send_buffer_row -INTEGER, DIMENSION(:), ALLOCATABLE :: send_buffer_col -INTEGER, DIMENSION(:), ALLOCATABLE :: recv_buffer_row -INTEGER, DIMENSION(:), ALLOCATABLE :: recv_buffer_col + !! Local Data - Send/Recv Buffers + INTEGER, DIMENSION(:), ALLOCATABLE :: send_buffer_row + INTEGER, DIMENSION(:), ALLOCATABLE :: send_buffer_col + INTEGER, DIMENSION(:), ALLOCATABLE :: recv_buffer_row + INTEGER, DIMENSION(:), ALLOCATABLE :: recv_buffer_col -!! Sizes help -INTEGER, DIMENSION(:), ALLOCATABLE :: recvcounts -INTEGER, DIMENSION(:), ALLOCATABLE :: displ -INTEGER :: gather_size + !! Sizes help + INTEGER, DIMENSION(:), ALLOCATABLE :: recvcounts + INTEGER, DIMENSION(:), ALLOCATABLE :: displ + INTEGER :: gather_size -!! Temporary variables -INTEGER :: num_processes, II, ierr + !! Temporary variables + INTEGER :: num_processes, II, ierr -!! Figure out the comm size -CALL MPI_COMM_SIZE(comm, num_processes, ierr) + !! Figure out the comm size + CALL MPI_COMM_SIZE(comm, num_processes, ierr) -!! Get the count -ALLOCATE(recvcounts(num_processes)) -CALL MPI_Allgather(triplet_in%CurrentSize, 1, MPI_INTEGER, recvcounts, & - & 1, MPI_INTEGER, comm, ierr) + !! Get the count + ALLOCATE(recvcounts(num_processes)) + CALL MPI_Allgather(triplet_in%CurrentSize, 1, MPI_INTEGER, recvcounts, & + & 1, MPI_INTEGER, comm, ierr) -!! Get the displacements -gather_size = SUM(recvcounts) -ALLOCATE(displ(num_processes)) -displ(1) = 0 -DO II = 2, num_processes - displ(II) = displ(II - 1) + recvcounts(II - 1) -END DO + !! Get the displacements + gather_size = SUM(recvcounts) + ALLOCATE(displ(num_processes)) + displ(1) = 0 + DO II = 2, num_processes + displ(II) = displ(II - 1) + recvcounts(II - 1) + END DO -!! Prepare the send buffers -ALLOCATE(send_buffer_row(triplet_in%CurrentSize)) -ALLOCATE(send_buffer_col(triplet_in%CurrentSize)) -ALLOCATE(send_buffer_val(triplet_in%CurrentSize)) -DO II = 1, triplet_in%CurrentSize - CALL GetTripletAt(triplet_in, II, temp_triplet) - send_buffer_row(II) = temp_triplet%index_row - send_buffer_col(II) = temp_triplet%index_column - send_buffer_val(II) = temp_triplet%point_value -END DO + !! Prepare the send buffers + ALLOCATE(send_buffer_row(triplet_in%CurrentSize)) + ALLOCATE(send_buffer_col(triplet_in%CurrentSize)) + ALLOCATE(send_buffer_val(triplet_in%CurrentSize)) + DO II = 1, triplet_in%CurrentSize + CALL GetTripletAt(triplet_in, II, temp_triplet) + send_buffer_row(II) = temp_triplet%index_row + send_buffer_col(II) = temp_triplet%index_column + send_buffer_val(II) = temp_triplet%point_value + END DO -!! Gather Call -ALLOCATE(recv_buffer_row(gather_size)) -ALLOCATE(recv_buffer_col(gather_size)) -ALLOCATE(recv_buffer_val(gather_size)) -CALL MPI_Allgatherv(send_buffer_row, triplet_in%CurrentSize, MPI_INTEGER, & - & recv_buffer_row, recvcounts, displ, MPI_INTEGER, comm, ierr) -CALL MPI_Allgatherv(send_buffer_col, triplet_in%CurrentSize, MPI_INTEGER, & - & recv_buffer_col, recvcounts, displ, MPI_INTEGER, comm, ierr) -CALL MPI_Allgatherv(send_buffer_val, triplet_in%CurrentSize, MPIDATATYPE, & - & recv_buffer_val, recvcounts, displ, MPIDATATYPE, comm, ierr) + !! Gather Call + ALLOCATE(recv_buffer_row(gather_size)) + ALLOCATE(recv_buffer_col(gather_size)) + ALLOCATE(recv_buffer_val(gather_size)) + CALL MPI_Allgatherv(send_buffer_row, triplet_in%CurrentSize, MPI_INTEGER, & + & recv_buffer_row, recvcounts, displ, MPI_INTEGER, comm, ierr) + CALL MPI_Allgatherv(send_buffer_col, triplet_in%CurrentSize, MPI_INTEGER, & + & recv_buffer_col, recvcounts, displ, MPI_INTEGER, comm, ierr) + CALL MPI_Allgatherv(send_buffer_val, triplet_in%CurrentSize, MPIDATATYPE, & + & recv_buffer_val, recvcounts, displ, MPIDATATYPE, comm, ierr) -!! Unpack -CALL ConstructTripletList(gathered_out, gather_size) -DO II = 1, gather_size - gathered_out%DATA(II)%index_row = recv_buffer_row(II) - gathered_out%DATA(II)%index_column = recv_buffer_col(II) - gathered_out%DATA(II)%point_value = recv_buffer_val(II) -END DO + !! Unpack + CALL ConstructTripletList(gathered_out, gather_size) + DO II = 1, gather_size + gathered_out%DATA(II)%index_row = recv_buffer_row(II) + gathered_out%DATA(II)%index_column = recv_buffer_col(II) + gathered_out%DATA(II)%point_value = recv_buffer_val(II) + END DO -!! Cleanup -DEALLOCATE(recvcounts) -DEALLOCATE(displ) -DEALLOCATE(send_buffer_row) -DEALLOCATE(send_buffer_col) -DEALLOCATE(send_buffer_val) + !! Cleanup + DEALLOCATE(recvcounts) + DEALLOCATE(displ) + DEALLOCATE(send_buffer_row) + DEALLOCATE(send_buffer_col) + DEALLOCATE(send_buffer_val)