Local version of scale with a diagonal

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

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);

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.

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

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

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

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.'''