Try to simplify the code

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,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))
@@ -170,33 +170,22 @@ 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
              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)
 
     !! 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)

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

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)

Source/Fortran/distributed_includes/GatherMatrixTripletList.f90

@@ -0,0 +1,3 @@
+CALL GatherMatrixToProcess(this, lmat)
+CALL MatrixToTripletList(lmat, tlist)
+CALL DestructMatrix(lmat)