Eigensolver Fallback Option (#187)

* Dense infastructure * Dense versions of standard functions * Turn lambda into a function of one variable * Dense density matrix solver * Linting for a more consistent feel * Singular value decomposition * Separate singular values into a new module Avoids cyclic dependency * Example directions fixed * My theory is that the test has a race condition * Try reactivating just the diagonal check of the eigenvalues * Test vectors too * Reactivate SVD test * Dense FOE for reference purposes * Forgot to add files * Removed testing line * Apparantely, those default values cause trouble * Fix capitalization * Fix chemical potential calculation * Fix cmake * Fix typo * Merge with master - fix * Merge in partial electrons * Fix a bug where we checked an unallocated variable's properties * Increase optional parameters for eigenexa (#190) * TRS4 seems faster * Nvals option for serial * Set number of values for eigenexa * Fix for optional vectors * Forgot conjugate in test
william-dawson · Apr 7, 2023 · de6067f · de6067f
1 parent 1605c9e
commit de6067f
Show file tree

Hide file tree

Showing 78 changed files with 3,770 additions and 1,225 deletions.
diff --git a/CMake/FindEigenExa.cmake b/CMake/FindEigenExa.cmake
@@ -0,0 +1,15 @@
+# Find the EigenExa Module
+# Variables set:
+# - EigenExa_FOUND - system found EigenExa.
+# - EigenExa_LIBRARIES - the linker line for EigenExa.
+# - EigenExa_INCLUDE_DIRS - the path to EigenExa.
+
+# First we search for the libraries
+find_library(EigenExa_LIBRARIES EigenExa)
+find_path(EigenExa_INCLUDE_DIRS "eigen_libs_mod.mod")
+
+# Now check if that worked
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(EigenExa DEFAULT_MSG
+  EigenExa_LIBRARIES EigenExa_INCLUDE_DIRS)
+mark_as_advanced(EigenExa_LIBRARIES EigenExa_INCLUDE_DIRS)
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -7,15 +7,23 @@ enable_language(Fortran)
 
 ################################################################################
 ## Packages
+set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/CMake)
 if (NOT CMAKE_TOOLCHAIN_FILE)
   message(WARNING "Building without a toolchain file. "
           "If this does not work, please see the example toolchain files in "
           "the Targets directory and set the parameters that match your "
           "system.")
   find_package(MPI REQUIRED)
   find_package(BLAS)
+  find_package(LAPACK)
 endif()
 find_package(SWIG 3.0)
+find_package(EigenExa)
+if (EigenExa_FOUND)
+  add_definitions("-DEIGENEXA")
+  set(EigenSolver_LIBRARIES ${EigenExa_LIBRARIES})
+  set(EigenSolver_INCLUDE_DIRS ${EigenExa_INCLUDE_DIRS})
+endif()
 
 ################################################################################
 ## Options

diff --git a/Examples/CMakeLinkage/CMakeLists.txt b/Examples/CMakeLinkage/CMakeLists.txt
@@ -13,6 +13,7 @@ find_package(NTPolyWrapper REQUIRED)
 find_package(NTPolyCPP REQUIRED)
 find_package(OpenMP)
 find_package(BLAS)
+find_package(LAPACK)
 
 ################################################################################
 ## Output Locations
@@ -22,7 +23,8 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
 ## Fortran Version
 add_executable(example_fortran ${CMAKE_SOURCE_DIR}/../PremadeMatrix/main.f90)
 target_link_libraries(example_fortran NTPoly::NTPoly ${MPI_Fortran_LIBRARIES}
-                      ${OpenMP_Fortran_LIBRARIES} ${BLAS_LIBRARIES})
+                      ${OpenMP_Fortran_LIBRARIES} ${LAPACK_LIBRARIES} 
+                      ${BLAS_LIBRARIES})
 target_include_directories(example_fortran PRIVATE ${MPI_INCLUDE_PATH})
 install(TARGETS example_fortran DESTINATION bin)
 

diff --git a/Examples/ComplexMatrix/ReadMe.md b/Examples/ComplexMatrix/ReadMe.md
@@ -38,7 +38,7 @@ Fortran Build Instructions:
 ```
 mpif90 main.f90 -o example \
   -I../../Build/include \
-  -L../../Build/lib -lNTPoly -fopenmp -lblas
+  -L../../Build/lib -lNTPoly -fopenmp -llapack -lblas
 
 ```
 
@@ -49,7 +49,7 @@ mpicxx main.cc -c \
 
 mpif90 main.o -o example \
   -L../../Build/lib -lNTPolyCPP -lNTPolyWrapper -lNTPoly -fopenmp -lstdc++ \
-  -lblas -lmpi_cxx
+  -llapack -lblas -lmpi_cxx
 
 ```
 

diff --git a/Examples/GraphTheory/ReadMe.md b/Examples/GraphTheory/ReadMe.md
@@ -22,7 +22,7 @@ Fortran Build Instructions:
 ```
 mpif90 main.f90 -o example \
   -I../../Build/include \
-  -L../../Build/lib -lNTPoly -fopenmp -lblas
+  -L../../Build/lib -lNTPoly -fopenmp -llapack -lblas
 
 ```
 
@@ -33,7 +33,7 @@ mpicxx main.cc -c \
 
 mpif90 main.o -o example \
   -L../../Build/lib -lNTPolyCPP -lNTPolyWrapper -lNTPoly -fopenmp -lstdc++ \
-  -lblas -lmpi_cxx
+  -llapack -lblas -lmpi_cxx
 
 ```
 

diff --git a/Examples/HydrogenAtom/ReadMe.md b/Examples/HydrogenAtom/ReadMe.md
@@ -49,7 +49,7 @@ Fortran Build Instructions:
 ```
 mpif90 main.f90 -o example \
   -I../../Build/include \
-  -L../../Build/lib -lNTPoly -fopenmp -lblas
+  -L../../Build/lib -lNTPoly -fopenmp -llapack -lblas
 
 ```
 C++ Build Instructions:
@@ -59,7 +59,7 @@ mpicxx main.cc -c \
 
 mpif90 main.o -o example \
   -L../../Build/lib -lNTPolyCPP -lNTPolyWrapper -lNTPoly -fopenmp -lstdc++ \
-  -lblas -lmpi_cxx
+  -llapack -lblas -lmpi_cxx
 
 ```
 

diff --git a/Examples/MatrixMaps/ReadMe.md b/Examples/MatrixMaps/ReadMe.md
@@ -28,7 +28,7 @@ Fortran Build Instructions:
 ```
 mpif90 main.f90 -o example \
   -I../../Build/include \
-  -L../../Build/lib -lNTPoly -fopenmp -lblas
+  -L../../Build/lib -lNTPoly -fopenmp -llapack -lblas
 
 ```
 
@@ -39,7 +39,7 @@ mpicxx main.cc -c \
 
 mpif90 main.o -o example \
   -L../../Build/lib -lNTPolyCPP -lNTPolyWrapper -lNTPoly -fopenmp -lstdc++ \
-  -lblas -lmpi_cxx
+  -llapack -lblas -lmpi_cxx
 
 ```
 

diff --git a/Examples/OverlapMatrix/ReadMe.md b/Examples/OverlapMatrix/ReadMe.md
@@ -61,7 +61,7 @@ Build with:
 ```
 mpif90 main.f90 -o example \
   -I../../Build/include \
-  -L../../Build/lib -lNTPoly -fopenmp -lblas
+  -L../../Build/lib -lNTPoly -fopenmp -llapack -lblas
 
 ```
 

diff --git a/Examples/PremadeMatrix/ReadMe.md b/Examples/PremadeMatrix/ReadMe.md
@@ -47,7 +47,7 @@ Fortran Build Instructions:
 ```
 mpif90 main.f90 -o example \
   -I../../Build/include \
-  -L../../Build/lib -lNTPoly -fopenmp -lblas
+  -L../../Build/lib -lNTPoly -fopenmp -llapack -lblas
 
 ```
 
@@ -58,7 +58,7 @@ mpicxx main.cc -c \
 
 mpif90 main.o -o example \
   -L../../Build/lib -lNTPolyCPP -lNTPolyWrapper -lNTPoly -fopenmp -lstdc++ \
-  -lblas -lmpi_cxx
+  -llapack -lblas -lmpi_cxx
 
 ```
 

diff --git a/ReadMe.md b/ReadMe.md
@@ -21,6 +21,7 @@ Installing NTPoly requires the following software:
 * An MPI Installation (MPI-3 Standard+).
 * CMake (Version 3.2+).
 * BLAS: for multiplying dense matrices, if they emerge in the calculation.
+* LAPACK: for dense solvers.
 
 The following optional software can greatly enhance the NTPoly experience:
 
@@ -32,6 +33,7 @@ The following optional software can greatly enhance the NTPoly experience:
 * MPI4PY: for testing.
 * SciPy: for testing.
 * NumPy: for testing.
+* EigenExa: for dense, parallel calculations.
 
 NTPoly uses CMake as a build system. First, take a look in the Targets
 directory. You'll find a list of `.cmake` files which have example
@@ -132,6 +134,7 @@ The following features and methods have been implemented in NTPoly:
 * Other
     * Matrix Inverse/Moore-Penrose Pseudo Inverse
     * Sign Function/Polar Decomposition
+    * Interface to Dense Eigen/Singular Value Decomposition
     * Load Balancing Matrices
     * File I/O
 

diff --git a/Source/C/CMakeLists.txt b/Source/C/CMakeLists.txt
@@ -3,6 +3,7 @@ set(Csrc
     ChebyshevSolvers_c.h
     DensityMatrixSolvers_c.h
     EigenBounds_c.h
+    EigenSolvers_c.h
     ExponentialSolvers_c.h
     GeometryOptimization_c.h
     HermiteSolvers_c.h

diff --git a/Source/C/DensityMatrixSolvers_c.h b/Source/C/DensityMatrixSolvers_c.h
@@ -26,6 +26,11 @@ void ScaleAndFold_wrp(const int *ih_Hamiltonian,
                       int *ih_Density, const double *homo, const double *lumo,
                       const double *energy_value_out,
                       const int *ih_solver_parameters);
+void DenseDensity_wrp(const int *ih_Hamiltonian,
+                      const int *ih_InverseSquareRoot, const double *trace,
+                      int *ih_Density, const double *energy_value_out,
+                      const double *chemical_potential_out,
+                      const int *ih_solver_parameters);
 void EnergyDensityMatrix_wrp(const int *ih_Hamiltonian, const int *ih_Density,
                              int *ih_EnergyDensity, const double *threshold);
 #endif
diff --git a/Source/C/EigenSolvers_c.h b/Source/C/EigenSolvers_c.h
@@ -0,0 +1,13 @@
+#ifndef EIGENSOLVERS_ch
+#define EIGENSOLVERS_ch
+
+void EigenDecomposition_wrp(const int *ih_this, int *ih_eigenvectors,
+                            const int *nvals, int *ih_eigenvalues,
+                            const int *ih_solver_parameters);
+void EigenDecomposition_novec_wrp(const int *ih_this, int *ih_eigenvectors,
+                                  const int *nvals,
+                                  const int *ih_solver_parameters);
+void SingularValueDecompostion_wrp(const int *ih_this, int *ih_leftvectors,
+                                   int *ih_rightvectors, int *ih_singularvalues,
+                                   const int *ih_solver_parameters);
+#endif
diff --git a/Source/C/ExponentialSolvers_c.h b/Source/C/ExponentialSolvers_c.h
@@ -3,9 +3,13 @@
 
 void ComputeExponential_wrp(const int *ih_Input, int *ih_Output,
                             const int *ih_solver_parameters);
+void ComputeDenseExponential_wrp(const int *ih_Input, int *ih_Output,
+                                 const int *ih_solver_parameters);
 void ComputeExponentialPade_wrp(const int *ih_Input, int *ih_Output,
                                 const int *ih_solver_parameters);
 void ComputeLogarithm_wrp(const int *ih_Input, int *ih_Output,
                           const int *ih_solver_parameters);
+void ComputeDenseLogarithm_wrp(const int *ih_Input, int *ih_Output,
+                               const int *ih_solver_parameters);
 
 #endif
diff --git a/Source/C/FermiOperator_c.h b/Source/C/FermiOperator_c.h
@@ -0,0 +1,10 @@
+#ifndef DENSITYMATRIXSOLVERS_ch
+#define DENSITYMATRIXSOLVERS_ch
+
+void ComputeDenseFOE_wrp(const int *ih_Hamiltonian,
+                         const int *ih_InverseSquareRoot, const double *trace,
+                         int *ih_Density, const double *inv_temp_in,
+                         const double *energy_value_out,
+                         const double *chemical_potential_out,
+                         const int *ih_solver_parameters);
+#endif
diff --git a/Source/C/InverseSolvers_c.h b/Source/C/InverseSolvers_c.h
@@ -3,6 +3,8 @@
 
 void Invert_wrp(const int *ih_Hamiltonian, int *ih_Inverse,
                 const int *ih_solver_parameters);
+void DenseInvert_wrp(const int *ih_Hamiltonian, int *ih_Inverse,
+                     const int *ih_solver_parameters);
 void PseudoInverse_wrp(const int *ih_Hamiltonian, int *ih_Inverse,
                        const int *ih_solver_parameters);
 #endif
diff --git a/Source/C/SignSolvers_c.h b/Source/C/SignSolvers_c.h
@@ -3,6 +3,8 @@
 
 void SignFunction_wrp(const int *ih_mat1, int *ih_signmat,
                       const int *ih_solver_parameters);
+void DenseSignFunction_wrp(const int *ih_mat1, int *ih_signmat,
+                           const int *ih_solver_parameters);
 void PolarDecomposition_wrp(const int *ih_mat1, int *ih_umat, int *ih_hmat,
                             const int *ih_solver_parameters);
 #endif
diff --git a/Source/C/SquareRootSolvers_c.h b/Source/C/SquareRootSolvers_c.h
@@ -3,7 +3,11 @@
 
 void SquareRoot_wrp(const int *ih_Input, int *ih_Output,
                     const int *ih_solver_parameters);
+void DenseSquareRoot_wrp(const int *ih_Input, int *ih_Output,
+                         const int *ih_solver_parameters);
 void InverseSquareRoot_wrp(const int *ih_Input, int *ih_Output,
                            const int *ih_solver_parameters);
+void DenseInverseSquareRoot_wrp(const int *ih_Input, int *ih_Output,
+                                const int *ih_solver_parameters);
 
 #endif
diff --git a/Source/C/TrigonometrySolvers_c.h b/Source/C/TrigonometrySolvers_c.h
@@ -3,7 +3,11 @@
 
 void Sine_wrp(const int *ih_Input, int *ih_Output,
               const int *ih_solver_parameters);
+void DenseSine_wrp(const int *ih_Input, int *ih_Output,
+                   const int *ih_solver_parameters);
 void Cosine_wrp(const int *ih_Input, int *ih_Output,
                 const int *ih_solver_parameters);
+void DenseCosine_wrp(const int *ih_Input, int *ih_Output,
+                     const int *ih_solver_parameters);
 
 #endif
diff --git a/Source/CPlusPlus/CMakeLists.txt b/Source/CPlusPlus/CMakeLists.txt
@@ -4,7 +4,9 @@ set(Csrc
     ChebyshevSolvers.cc
     DensityMatrixSolvers.cc
     EigenBounds.cc
+    EigenSolvers.cc
     ExponentialSolvers.cc
+    FermiOperator.cc
     GeometryOptimization.cc
     HermiteSolvers.cc
     InverseSolvers.cc
@@ -34,7 +36,9 @@ set(Chead
     ChebyshevSolvers.h
     DensityMatrixSolvers.h
     EigenBounds.h
+    EigenSolvers.h
     ExponentialSolvers.h
+    FermiOperator.h
     GeometryOptimization.h
     HermiteSolvers.h
     InverseSolvers.h

diff --git a/Source/CPlusPlus/DensityMatrixSolvers.cc b/Source/CPlusPlus/DensityMatrixSolvers.cc
@@ -52,6 +52,16 @@ void DensityMatrixSolvers::HPCP(const Matrix_ps &Hamiltonian,
            GetIH(solver_parameters));
 }
 
+////////////////////////////////////////////////////////////////////////////////
+void DensityMatrixSolvers::DenseDensity(
+    const Matrix_ps &Hamiltonian, const Matrix_ps &Overlap, double trace,
+    Matrix_ps &Density, double &energy_value_out,
+    double &chemical_potential_out, const SolverParameters &solver_parameters) {
+  DenseDensity_wrp(GetIH(Hamiltonian), GetIH(Overlap), &trace, GetIH(Density),
+                   &energy_value_out, &chemical_potential_out,
+                   GetIH(solver_parameters));
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 void DensityMatrixSolvers::ScaleAndFold(
     const Matrix_ps &Hamiltonian, const Matrix_ps &Overlap, double trace,

diff --git a/Source/CPlusPlus/DensityMatrixSolvers.h b/Source/CPlusPlus/DensityMatrixSolvers.h
@@ -84,6 +84,19 @@ class DensityMatrixSolvers : public SolverBase {
                            Matrix_ps &Density, double homo, double lumo,
                            double &energy_value_out,
                            const SolverParameters &solver_parameters);
+  //! Compute the density matrix using a dense solver.
+  //!\param Hamiltonian the matrix to compute the corresponding density from.
+  //!\param InverseSquareRoot of the overlap matrix.
+  //!\param trace of the density matrix (usually the number of electrons).
+  //!\param Density the density matrix computed by this routine.
+  //!\param energy_value_out the energy of the system (optional).
+  //!\param chemical_potential_out the chemical potential calculated.
+  //!\param solver_parameters parameters for the solver
+  static void DenseDensity(const Matrix_ps &Hamiltonian,
+                           const Matrix_ps &InverseSquareRoot, double trace,
+                           Matrix_ps &Density, double &energy_value_out,
+                           double &chemical_potential_out,
+                           const SolverParameters &solver_parameters);
   //! Compute the energy-weighted density matrix.
   //!\param Hamiltonian the matrix to compute from.
   //!\param Density the density matrix.

diff --git a/Source/CPlusPlus/EigenSolvers.cc b/Source/CPlusPlus/EigenSolvers.cc
@@ -0,0 +1,33 @@
+#include "EigenSolvers.h"
+using namespace NTPoly;
+
+////////////////////////////////////////////////////////////////////////////////
+extern "C" {
+#include "EigenSolvers_c.h"
+}
+
+////////////////////////////////////////////////////////////////////////////////
+namespace NTPoly {
+////////////////////////////////////////////////////////////////////////////////
+void EigenSolvers::EigenDecomposition(
+    const Matrix_ps &matrix, Matrix_ps &eigenvalues, int nvals,
+    Matrix_ps &eigenvectors, const SolverParameters &solver_parameters) {
+  EigenDecomposition_wrp(GetIH(matrix), GetIH(eigenvalues), &nvals,
+                         GetIH(eigenvectors), GetIH(solver_parameters));
+}
+////////////////////////////////////////////////////////////////////////////////
+void EigenSolvers::EigenValues(const Matrix_ps &matrix, Matrix_ps &eigenvalues,
+                               int nvals,
+                               const SolverParameters &solver_parameters) {
+  EigenDecomposition_novec_wrp(GetIH(matrix), GetIH(eigenvalues), &nvals,
+                               GetIH(solver_parameters));
+}
+////////////////////////////////////////////////////////////////////////////////
+void EigenSolvers::SingularValueDecomposition(
+    const Matrix_ps &matrix, Matrix_ps &leftvectors, Matrix_ps &rightvectors,
+    Matrix_ps &singularvalues, const SolverParameters &solver_parameters) {
+  SingularValueDecompostion_wrp(GetIH(matrix), GetIH(leftvectors),
+                                GetIH(rightvectors), GetIH(singularvalues),
+                                GetIH(solver_parameters));
+}
+} // namespace NTPoly