WIP: Refactor

CMakeLists.txt

@@ -20,7 +20,11 @@ endif(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CXX_FLAGS)
 ###############################################################################
 # User input options
 ###############################################################################
-option(ENABLE_TESTING "Build and enable tests" OFF)
+option(ENABLE_TESTING "Build and enable tests" ON)
+if(ENABLE_TESTING)
+  enable_testing()
+  find_package(Boost REQUIRED COMPONENTS unit_test_framework)
+endif(ENABLE_TESTING)
 
 ###############################################################################
 # Dependencies

clang.sh

@@ -1,2 +1,3 @@
-clang-format -i src/*.cpp
+clang-format -i src/*.cc
+clang-format -i src/tests/*.cc
 clang-format -i include/*.hpp

include/eigenValueSolver.hpp

@@ -21,39 +21,43 @@
 
 namespace eigenValueSolver {
 
-template <typename real> class eigenSolver {
+template <typename real>
+class eigenSolver {
 
-public:
+ public:
   /**
    * Groups all the options for the solvers
    */
   struct options {
-    std::string solver = "davidson"; // either davidson or jacobi
-    real tolerence = 1e-8;           /** tolerence for convergence*/
+    std::string solver = "davidson";  // either davidson or jacobi
+    real tolerence = 1e-8;            /** tolerence for convergence*/
     int numberOfEigenValues = 1; /** number of eigen values to be calculated*/
     int sizeOfTheMatrix = 100;   /**Size of the input matrix*/
   };
 
+  El::DistMatrix<real> eigenVectors; /**s stores the current eigenvector
+                                        matrix*/
+  El::DistMatrix<real> eigenVectorsFull; /**s stores the current eigenvector
+                                                                              matrix*/
+  El::DistMatrix<real, El::VR, El::STAR> eigenValues; /**vector of eigen
+                                                         values*/
+
   options solverOptions;
   /**
    * This function computes the eigen value-vector pairs for the input matrix
    */
-  void solve(El::DistMatrix<real> &, El::Grid &);
+  void solve(const El::DistMatrix<real> &, El::Grid &);
 
-private:
+ private:
   int columnsOfSearchSpace = solverOptions.numberOfEigenValues *
                              2; /**The columns of the search space*/
   /*************Initialise all the matrices necessary**************/
   El::DistMatrix<real> searchSpace;      /**Matrix that holds the search space*/
   El::DistMatrix<real> searchSpacesub;   /**Matrix that holds the search space*/
   El::DistMatrix<real> correctionVector; /**a matrix that holds the k guess
                                             eigen vectors each of size nx1*/
-  El::DistMatrix<real>
-      eigenVectors; /**s stores the current eigenvector matrix*/
-  El::DistMatrix<real, El::VR, El::STAR>
-      eigenValues; /**vector of eigen values*/
-  El::DistMatrix<real, El::VR, El::STAR>
-      eigenValues_old; /**vector of eigen values*/
+  El::DistMatrix<real, El::VR, El::STAR> eigenValues_old; /**vector of eigen
+                                                             values*/
 
   // Range to loop over just the required number of eigenvalues in theta
   El::Range<int> begTheta;
@@ -79,13 +83,13 @@ template <typename real> class eigenSolver {
   /**
    *Expands the search space with the correction vector
    */
-  void expandSearchSpace(int, El::DistMatrix<real> &, El::Grid &);
+  void expandSearchSpace(int, const El::DistMatrix<real> &, El::Grid &);
 
   /**
    * solve the subspace problem i.e. VTAV and eigenvalue/vectors
    */
-  void subspaceProblem(int, El::DistMatrix<real> &, El::Grid &);
+  void subspaceProblem(int, const El::DistMatrix<real> &, El::Grid &);
 };
-} // namespace eigenValueSolver
+}  // namespace eigenValueSolver
 
 #endif /* eigenValueSolver_h */

src/CMakeLists.txt

@@ -1,5 +1,5 @@
-add_library(eigenValueSolver eigenValueSolver.cpp)
-add_executable(main main.cpp utils.cpp )
+add_library(eigenValueSolver eigenValueSolver.cc utils.cc)
+add_executable(main main.cc)
 
 target_include_directories(eigenValueSolver
   PUBLIC
@@ -13,3 +13,7 @@ target_link_libraries(eigenValueSolver
 target_link_libraries(main
   PUBLIC
   eigenValueSolver)
+
+if (ENABLE_TESTING)
+  add_subdirectory(tests)
+endif()

src/eigenValueSolver.cpp → src/eigenValueSolver.cc

@@ -9,12 +9,14 @@
 #include "eigenValueSolver.hpp"
 
 namespace eigenValueSolver {
-template <typename real> void eigenSolver<real>::initialise(El::Grid &grid) {
+template <typename real>
+void eigenSolver<real>::initialise(El::Grid &grid) {
 
   searchSpace.SetGrid(grid);
   searchSpacesub.SetGrid(grid);
   correctionVector.SetGrid(grid);
   eigenVectors.SetGrid(grid);
+  eigenVectorsFull.SetGrid(grid);
   eigenValues.SetGrid(grid);
   eigenValues_old.SetGrid(grid);
 
@@ -29,10 +31,13 @@ template <typename real> void eigenSolver<real>::initialise(El::Grid &grid) {
                columnsOfSearchSpace);
   El::Identity(eigenVectors, solverOptions.sizeOfTheMatrix,
                solverOptions.sizeOfTheMatrix);
+  El::Identity(eigenVectorsFull, solverOptions.sizeOfTheMatrix,
+                            solverOptions.sizeOfTheMatrix);
+    El::Identity(eigenValues, solverOptions.sizeOfTheMatrix, 1);
 }
 
 template <typename real>
-void eigenSolver<real>::subspaceProblem(int iterations, El::DistMatrix<real> &A,
+void eigenSolver<real>::subspaceProblem(int iterations, const El::DistMatrix<real> &A,
                                         El::Grid &grid) {
   searchSpacesub = searchSpace;
   searchSpacesub.Resize(solverOptions.sizeOfTheMatrix, iterations + 1);
@@ -41,6 +46,7 @@ void eigenSolver<real>::subspaceProblem(int iterations, El::DistMatrix<real> &A,
   El::DistMatrix<real> Ttemp(grid), T(grid);
   El::Zeros(Ttemp, solverOptions.sizeOfTheMatrix, iterations + 1);
   El::Zeros(T, iterations + 1, iterations + 1);
+  El::Zeros(eigenVectorsFull, solverOptions.sizeOfTheMatrix, iterations + 1);
 
   // Parameters for GEMM
   real alpha = 1, beta = 0;
@@ -52,18 +58,21 @@ void eigenSolver<real>::subspaceProblem(int iterations, El::DistMatrix<real> &A,
 
   // Get the eigen pairs for the reduced problem V^TAV
   El::HermitianEig(El::UPPER, T, eigenValues, eigenVectors);
+
+  //Calculate the eigenvectors of the full matrix
+  El::Gemm(El::NORMAL, El::NORMAL, alpha, searchSpacesub, eigenVectors, beta, eigenVectorsFull);
 }
 
 template <typename real>
 void eigenSolver<real>::expandSearchSpace(int iterations,
-                                          El::DistMatrix<real> &A,
+                                          const El::DistMatrix<real> &A,
                                           El::Grid &grid) {
 
   for (int j = 0; j < columnsOfSearchSpace; ++j) {
     El::Range<int> beg(0, iterations + 1);
     El::Range<int> end(j, j + 1);
 
-    El::DistMatrix<real> residual(grid); // residual Ay-thetay
+    El::DistMatrix<real> residual(grid);  // residual Ay-thetay
     El::Zeros(residual, solverOptions.sizeOfTheMatrix, 1);
 
     // calculate the ritz vector Vs
@@ -73,39 +82,39 @@ void eigenSolver<real>::expandSearchSpace(int iterations,
     El::Gemv(El::NORMAL, alpha, searchSpacesub, eigenVectors(beg, end), beta,
              Vs);
 
-    El::DistMatrix<real> I(grid); // Identitiy matrix
+    El::DistMatrix<real> I(grid);  // Identitiy matrix
     El::Identity(
         I, solverOptions.sizeOfTheMatrix,
-        solverOptions.sizeOfTheMatrix); // Initialize as identity matrix
+        solverOptions.sizeOfTheMatrix);  // Initialize as identity matrix
     I *= eigenValues.GetLocal(j, 0);
     El::DistMatrix<real> Atemp(A);
-    Atemp -= I; // A-theta*I
+    Atemp -= I;  // A-theta*I
 
     // Calculate the residual r=(A-theta*I)*Vs
     El::Gemv(El::NORMAL, alpha, Atemp, Vs, beta, residual);
 
     if (solverOptions.solver == "davidson") {
       real den = 1.0 / eigenValues.GetLocal(j, 0) - A.GetLocal(j, j);
 
-      correctionVector = residual; // new search direction
+      correctionVector = residual;  // new search direction
       correctionVector *= den;
     } else if (solverOptions.solver == "jacobi") {
-      El::DistMatrix<real> proj(grid); // projector matrix
+      El::DistMatrix<real> proj(grid);  // projector matrix
       El::Zeros(proj, solverOptions.sizeOfTheMatrix,
                 solverOptions.sizeOfTheMatrix);
-      El::DistMatrix<real> Ip(grid); // Identitiy matrix
+      El::DistMatrix<real> Ip(grid);  // Identitiy matrix
       El::Identity(Ip, solverOptions.sizeOfTheMatrix,
                    solverOptions.sizeOfTheMatrix);
 
       El::Gemm(El::NORMAL, El::TRANSPOSE, alpha, Vs, Vs, beta, proj);
       proj -= I;
       proj *= -1.0;
 
-      El::DistMatrix<real> projProd(grid); // product of the projectors
+      El::DistMatrix<real> projProd(grid);  // product of the projectors
       El::Zeros(projProd, solverOptions.sizeOfTheMatrix,
                 solverOptions.sizeOfTheMatrix);
 
-      El::DistMatrix<real> projTemp(grid); // temp intermediate step
+      El::DistMatrix<real> projTemp(grid);  // temp intermediate step
       El::Zeros(projTemp, solverOptions.sizeOfTheMatrix,
                 solverOptions.sizeOfTheMatrix);
 
@@ -132,16 +141,16 @@ void eigenSolver<real>::expandSearchSpace(int iterations,
 }
 
 template <typename real>
-void eigenSolver<real>::solve(El::DistMatrix<real> &A, El::Grid &grid) {
+void eigenSolver<real>::solve(const El::DistMatrix<real> &A, El::Grid &grid) {
 
   eigenSolver<real>::initialise(grid);
 
   int maximumIterations = solverOptions.sizeOfTheMatrix / 2;
 
   for (int iterations = columnsOfSearchSpace; iterations < maximumIterations;
        iterations = iterations + columnsOfSearchSpace) {
-    if (iterations <=
-        columnsOfSearchSpace) // If it is the first iteration copy t to V
+    if (iterations <= columnsOfSearchSpace)  // If it is the first iteration
+                                             // copy t to V
     {
 
       for (int i = 0; i < solverOptions.sizeOfTheMatrix; ++i) {
@@ -150,14 +159,14 @@ void eigenSolver<real>::solve(El::DistMatrix<real> &A, El::Grid &grid) {
         }
       }
       El::Ones(eigenValues_old, solverOptions.sizeOfTheMatrix,
-               1); // so this not to converge immediately
-    } else // if its not the first iteration then set old theta to the new one
+               1);  // so this not to converge immediately
+    } else  // if its not the first iteration then set old theta to the new one
     {
       eigenValues_old = eigenValues(begTheta, endTheta);
     }
 
     // Orthogonalize the searchSpace matrix using QR
-    El::DistMatrix<real> R; // R matrix for QR factorization
+    El::DistMatrix<real> R;  // R matrix for QR factorization
     El::Zeros(R, solverOptions.sizeOfTheMatrix, solverOptions.sizeOfTheMatrix);
 
     // QR factorization of V
@@ -167,11 +176,12 @@ void eigenSolver<real>::solve(El::DistMatrix<real> &A, El::Grid &grid) {
     subspaceProblem(iterations, A, grid);
     // expand the search space
     expandSearchSpace(iterations, A, grid);
+
   }
-  El::Print(eigenValues(begTheta, endTheta));
+    //El::Print(eigenVectors);
 }
 
 // explicit instantiations
 template class eigenSolver<float>;
 template class eigenSolver<double>;
-} // namespace eigenValueSolver
+}  // namespace eigenValueSolver

src/main.cpp → src/main.cc

@@ -37,15 +37,14 @@ int main(int argc, char *argv[]) {
     const El::Int matrixSize = El::Input("--size", "size of matrix", 100);
     const El::Int numEig = El::Input("--numeig", "number of eigenvalues", 1);
     const std::string solverType =
-        El::Input("--solver", "solver used", "Davidson");
+        El::Input("--solver", "solver used", "davidson");
 
     // Set block size
     El::SetBlocksize(blocksize);
 
     // If the grid height wasn't specified, then we should attempt to build
     // a nearly-square process grid
-    if (gridHeight == 0)
-      gridHeight = El::Grid::DefaultHeight(commSize);
+    if (gridHeight == 0) gridHeight = El::Grid::DefaultHeight(commSize);
     El::Grid grid{comm, gridHeight};
     if (commRank == 0)
       El::Output("Grid is: ", grid.Height(), " x ", grid.Width());

src/tests/test_eigen.cc

@@ -0,0 +1,58 @@
+#define BOOST_TEST_MODULE eigen_solver
+#include "eigenValueSolver.hpp"
+#include "utils.hpp"
+#include <boost/test/included/unit_test.hpp>
+
+using namespace eigenValueSolver;
+using namespace El;
+
+/**test the davidson implementation*/
+BOOST_AUTO_TEST_CASE(davidson_solver) {
+  // initialise elemental mpi
+  El::Environment env;
+  El::mpi::Comm comm = El::mpi::COMM_WORLD;
+  const int commRank = El::mpi::Rank(comm);
+  const int commSize = El::mpi::Size(comm);
+  typedef double real;
+  El::Int gridHeight = 0;
+  gridHeight = El::Grid::DefaultHeight(commSize);
+  El::Grid grid{comm, gridHeight};
+
+  // The matrix A whose eigenvalues have to be evaluated
+  real matrixSize = 100;
+  El::DistMatrix<real> A(grid);
+  El::DistMatrix<real> Ax(grid);
+  El::DistMatrix<real> lambdax(grid);
+  El::Zeros(A, matrixSize, matrixSize);
+  El::Zeros(Ax, matrixSize, 1);
+  El::Zeros(lambdax, matrixSize, 1);
+
+  // Generate the Diagonally dominant hermitian matrix
+  generateDDHermitianMatrix<real>(A);
+  //El::Print(A);
+  eigenSolver<real> solver;
+  solver.solverOptions.numberOfEigenValues = 4;
+  solver.solverOptions.tolerence = 1e-8;
+  solver.solverOptions.solver = "davidson";
+  solver.solverOptions.sizeOfTheMatrix = A.Height();
+
+  solver.solve(A, grid);
+
+  /**Check if Ax-lambda.x < 1e-6 holds!!
+  */
+  El::Range<int> beg(0, solver.solverOptions.sizeOfTheMatrix);
+  El::Range<int> end(0, 1);
+
+  double eVal = solver.eigenValues.GetLocal(0,0);
+  El::DistMatrix<real> eVec = solver.eigenVectorsFull(beg,end);
+
+  El::Gemm(El::NORMAL, El::NORMAL, 1.0, A, solver.eigenVectorsFull(beg,end), 0.0, Ax);
+  El::Scale(eVal, eVec);
+
+
+  El::DistMatrix<real> r = Ax;
+  r -= eVec;
+  double rNorm = El::Nrm2(r);
+
+  BOOST_CHECK(rNorm < 1e-6);
+}