Feature for Scaling the Matrix by a Diagonal matrix

Source/C/PSMatrix_c.h

@@ -62,5 +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_psr_wrp(int *ih_mat, const int *ih_tlist);
+void MatrixDiagonalScale_psc_wrp(int *ih_mat, const int *ih_tlist);
 
 #endif

Source/C/SMatrix_c.h

@@ -34,6 +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 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);
@@ -70,5 +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 MatrixDiagonalScale_lsc_wrp(int *ih_mat, const int *ih_tlist);
 
 #endif

Source/CPlusPlus/PSMatrix.cc

@@ -217,6 +217,16 @@ void Matrix_ps::Scale(double constant) {
   ScaleMatrix_ps_wrp(ih_this, &constant);
 }
 
+//////////////////////////////////////////////////////////////////////////////
+void Matrix_ps::DiagonalScale(const TripletList_r &tlist) {
+  MatrixDiagonalScale_psr_wrp(ih_this, tlist.ih_this);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+void Matrix_ps::DiagonalScale(const TripletList_c &tlist) {
+  MatrixDiagonalScale_psc_wrp(ih_this, tlist.ih_this);
+}
+
 //////////////////////////////////////////////////////////////////////////////
 double Matrix_ps::Norm() const { return MatrixNorm_ps_wrp(ih_this); }
 

Source/CPlusPlus/PSMatrix.h

@@ -160,13 +160,19 @@ 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);
+  //!\param tlist the triplet list.
+  //!\param threshold for flushing small values.
+  void DiagonalScale(const NTPoly::TripletList_c &tlist);
 
 public:
   //! Destructor.

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) {
+  MatrixDiagonalScale_lsr_wrp(ih_this, tlist.ih_this);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Matrix_lsc::DiagonalScale(const TripletList_c &tlist) {
+  MatrixDiagonalScale_lsc_wrp(ih_this, tlist.ih_this);
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 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);
 
 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);
 
 public:
   //! Transpose a sparse matrix.

Source/Fortran/EigenSolversModule.F90

@@ -13,18 +13,18 @@ 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, &
        & CopySolverParameters
   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

Source/Fortran/FermiOperatorModule.F90

@@ -1,23 +1,26 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !> 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, &
        & 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, &
+       & GatherMatrixTripletList, GetMatrixSize
   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, TripletList_c, &
+       & DestructTripletList, CopyTripletList
   USE NTMPIModule
   IMPLICIT NONE
   PRIVATE
@@ -55,13 +58,13 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, &
     TYPE(Matrix_ps) :: vecs, vecsT, vals, Temp
     TYPE(MatrixMemoryPool_p) :: pool
     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
@@ -100,23 +103,23 @@ 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
        ALLOCATE(occ(num_eigs))
        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))
@@ -148,7 +151,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,30 +162,37 @@ 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
        ELSE
           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
-          tlist%DATA(II)%point_value = occ_temp
+          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
+             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)
 
-    !! Fill
-    CALL ConstructEmptyMatrix(vals, H)
-    CALL FillMatrixFromTripletList(vals, tlist, preduplicated_in = .TRUE.)
+    !! Scale the eigenvectors
+    IF (vecs%is_complex) THEN
+       CALL CopyTripletList(tlist, tlist_c)
+       CALL MatrixDiagonalScale(vecs, tlist_c)
+    ELSE
+       CALL MatrixDiagonalScale(vecs, tlist)
+    END IF
+    CALL FilterMatrix(vecs, params%threshold)
 
     !! 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 +218,7 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, &
     CALL DestructMatrix(vals)
     CALL DestructMatrix(temp)
     CALL DestructTripletList(tlist)
+    CALL DestructTripletList(tlist_c)
     CALL DestructMatrixMemoryPool(pool)
     IF (ALLOCATED(occ)) THEN
        DEALLOCATE(occ)

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)