diff --git a/.gitignore b/.gitignore index 067c8e5..6909559 100644 --- a/.gitignore +++ b/.gitignore @@ -50,3 +50,6 @@ lib/ #clang clang.sh + +#visual studio +.vscode diff --git a/.travis.yml b/.travis.yml index 44f8123..c5dc117 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,9 +1,17 @@ -sudo: true +sudo: false language: cpp - -os: - - linux + +addons: + apt: + sources: + - ubuntu-toolchain-r-test + packages: + - cmake + - libblas-dev + - liblapack-dev + - gfortran-7 + - libopenmpi-dev dist: bionic @@ -16,7 +24,14 @@ branches: before_install: - sudo add-apt-repository ppa:ubuntu-toolchain-r/test -y - sudo apt-get update -q - - sudo apt-get install cmake libblas-dev liblapack-dev libopenmpi-dev + # Install elemental + - git clone --branch master git://github.com/elemental/Elemental.git ${HOME}/elemental + - ls /usr/lib/gcc/x86_64-linux-gnu/7 + - | + cd ${HOME}/elemental && + cmake -H. -Bbuild -DCMAKE_INSTALL_PREFIX=${HOME}/local -DCMAKE_Fortran_COMPILER=/usr/bin/gfortran-7 -DGFORTRAN_LIB=/usr/lib/gcc/x86_64-linux-gnu/7 -DEL_C_INTERFACE=OFF -DEL_BUILT_OPENBLAS=OFF && + cmake --build build --target install script: - cmake -H. -Bbuild - cmake --build build + - ctest --output-on-failure \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 77fc6ac..ed3590a 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/README.md b/README.md new file mode 100644 index 0000000..c115ac2 --- /dev/null +++ b/README.md @@ -0,0 +1,60 @@ +Davidson Eigensolver +=================== +This package contains a C++ implementation of the *Davidson diagonalization algorithms*. +Different schemas are available to compute the correction. + +Available correction methods are: + * **davidson**: Diagonal-Preconditioned-Residue + * **jacobi**: Generalized Jacobi Davidson + + ### Note: +The Davidson method is suitable for **diagonal-dominant symmetric matrices**, that are quite common +in certain scientific problems like [electronic structure](https://en.wikipedia.org/wiki/Electronic_structure). The Davidson method could be not practical +for other kind of symmetric matrices.i + +Usage +----- +The following program calls the `solve` *subroutine* from the `davidson` module and computes +the lowest eigenvalue and corresponding eigenvector, using the *jacobi* correction method with a tolerance of `1e-8`. +```C++ +eigenSolver solver; +solver.solverOptions.numberOfEigenValues = numEig; +solver.solverOptions.tolerence = 1e-8; +solver.solverOptions.solver = "jacobi"; +solver.solverOptions.sizeOfTheMatrix = A.Height(); +solver.solve(A, grid); +``` + +The helper `generateDDHermitianMatrix` function generates a diagonal dominant matrix. + +**Variables**: + * `A` (*in*) matrix to diagonalize + * `solver.eigenValues` (*out*) resulting eigenvalues + * `solver.eigenVectorsFull` (*out*) resulting eigenvectors + * `solver.SolverOptions.solver`(*in*) Either "davidson" or "jacobi" + * `solver.solverOptions.tolerance`(*in*) Numerical tolerance for convergence + +### References: + * [Davidson diagonalization method and its applications to electronic structure calculations](https://www.semanticscholar.org/paper/DAVIDSON-DIAGONALIZATION-METHOD-AND-ITS-APPLICATION-Liao/5811eaf768d1a006f505dfe24f329874a679ba59) + * [Numerical Methods for Large Eigenvalue Problem](https://doi.org/10.1137/1.9781611970739) + +Installation and Testing +------------------------ + +To compile execute: +``` +cmake -H. -Bbuild && cmake --build build +``` + +To run the test: +``` +cmake -H. -Bbuild -DENABLE_TEST=ON && cmake --build build +cd build && ctest -V +``` + +Dependencies +------------ +This packages assumes that you have installed the following packages: + * A C++ compiler >= C++ 14 + * [CMake](https://cmake.org/) + * [Elemental] (https://github.com/elemental/Elemental) diff --git a/clang.sh b/clang.sh index 0c40577..71c705e 100644 --- a/clang.sh +++ b/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 diff --git a/include/eigenValueSolver.hpp b/include/eigenValueSolver.hpp index f68aada..165ecea 100644 --- a/include/eigenValueSolver.hpp +++ b/include/eigenValueSolver.hpp @@ -21,71 +21,81 @@ namespace eigenValueSolver { -template class eigenSolver { +template +class eigenSolver { -public: - /** - * Groups all the options for the solvers - */ + 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*/ - int numberOfEigenValues = 1; /** number of eigen values to be calculated*/ - int sizeOfTheMatrix = 100; /**Size of the input matrix*/ + // Type of solver used - either "davidson" or "jacobi" + std::string solver = "davidson"; + + // Tolerence criterion for convergence + real tolerence = 1e-8; + + // The number of eigen values to be calculated + int numberOfEigenValues = 1; + + // Size of the input matrix + int sizeOfTheMatrix = 100; }; + // MEMBER VARIABLES + + // Eigenvector matrix of the reduced problem + El::DistMatrix eigenVectors; + + // ritzVectors i.e. eigen vectors of the full matrix + El::DistMatrix ritzVectors; + + // Eigenvalues + El::DistMatrix eigenValues; + + // The guess eigenspace + El::DistMatrix V; + + // A*V + El::DistMatrix AV; + + // The guess eigenspace per iteration + El::DistMatrix Vsub; + + // The correction vector + El::DistMatrix correctionVector; + + // The residual + El::DistMatrix residual; + + // Instance of the options struct options solverOptions; - /** - * This function computes the eigen value-vector pairs for the input matrix - */ - void solve(El::DistMatrix &, El::Grid &); - -private: - int columnsOfSearchSpace = solverOptions.numberOfEigenValues * - 2; /**The columns of the search space*/ - /*************Initialise all the matrices necessary**************/ - El::DistMatrix searchSpace; /**Matrix that holds the search space*/ - El::DistMatrix searchSpacesub; /**Matrix that holds the search space*/ - El::DistMatrix correctionVector; /**a matrix that holds the k guess - eigen vectors each of size nx1*/ - El::DistMatrix - eigenVectors; /**s stores the current eigenvector matrix*/ - El::DistMatrix - eigenValues; /**vector of eigen values*/ - El::DistMatrix - eigenValues_old; /**vector of eigen values*/ + + // MEMBER FUNCTIONS + // This function computes the eigen value-vector pairs for the input matrix + void solve(const El::DistMatrix &, El::Grid &); + + private: + // The size of the search space + real sizeOfSearchSpace; // Range to loop over just the required number of eigenvalues in theta El::Range begTheta; El::Range endTheta; - // El::Range begTheta{0, solverOptions.numberOfEigenValues}; - // El::Range endTheta{0, 1}; - - /** - * This function initialises all the required matrices - */ - void initialise(El::Grid &); - - /** - *Calculates the residual - */ - void calculateResidual(); - - /** - *Calculates the correction vector - */ - void calculateCorrectionVector(); - - /** - *Expands the search space with the correction vector - */ - void expandSearchSpace(int, El::DistMatrix &, El::Grid &); - - /** - * solve the subspace problem i.e. VTAV and eigenvalue/vectors - */ - void subspaceProblem(int, El::DistMatrix &, El::Grid &); + + // This function initialises all the required matrices + void initialise(El::Grid &, const El::DistMatrix &); + + // Calculates the residual + bool calculateResidual(int, El::Grid &); + + // Calculates the correction vector + void calculateCorrectionVector(int, const El::DistMatrix &, El::Grid &); + + // Expands the search space with the correction vector + void expandSearchSpace(int &, const El::DistMatrix &, El::Grid &); + + // Solve the subspace problem i.e. VTAV and eigenvalue/vectors + void subspaceProblem(int, const El::DistMatrix &, El::Grid &); }; -} // namespace eigenValueSolver +} // namespace eigenValueSolver #endif /* eigenValueSolver_h */ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ae693df..cf565e1 100644 --- a/src/CMakeLists.txt +++ b/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() diff --git a/src/eigenValueSolver.cc b/src/eigenValueSolver.cc new file mode 100644 index 0000000..076dfa0 --- /dev/null +++ b/src/eigenValueSolver.cc @@ -0,0 +1,264 @@ +// +// eigenValueSolver.cpp +// eigenValueSolver +// +// Created by Adithya Vijaykumar on 17/09/2019. +// Copyright © 2019 Adithya Vijaykumar. All rights reserved. +// + +#include "eigenValueSolver.hpp" +#include + +namespace eigenValueSolver { + +template +void eigenSolver::initialise(El::Grid &grid, + const El::DistMatrix &A) { + // Intitialize the size of the search space + sizeOfSearchSpace = solverOptions.numberOfEigenValues * 3; + + // Intialize the search space V as an identity matrix of size n x + // sizeOfSearchSpace + V.SetGrid(grid); + El::Identity(V, solverOptions.sizeOfTheMatrix, sizeOfSearchSpace); + + // Initialize the Vsub which part of the search space relevent to the current + // iteration of size n x numEigVal + Vsub.SetGrid(grid); + El::Identity(Vsub, solverOptions.sizeOfTheMatrix, + solverOptions.numberOfEigenValues); + + // Initialize the eigen vectors + eigenVectors.SetGrid(grid); + El::Identity(eigenVectors, solverOptions.sizeOfTheMatrix, + solverOptions.numberOfEigenValues); + + // Initalize the residuals for all the eigenvalues + residual.SetGrid(grid); + El::Identity(residual, solverOptions.sizeOfTheMatrix, + solverOptions.numberOfEigenValues); + + // Set the eigen value vector to diagonal values of input matrix and calculate + // respective residuals + eigenValues.SetGrid(grid); + El::Identity(eigenValues, solverOptions.numberOfEigenValues, 1); + + ritzVectors.SetGrid(grid); + AV.SetGrid(grid); + +} // namespace eigenValueSolver + +// Calculate the residual and the correction vector +template +void eigenSolver::calculateCorrectionVector(int iterations, + const El::DistMatrix &A, + El::Grid &grid) { + correctionVector.SetGrid(grid); + El::Identity(correctionVector, solverOptions.sizeOfTheMatrix, + solverOptions.numberOfEigenValues); + for (int j = 0; j < solverOptions.numberOfEigenValues; ++j) { + if (solverOptions.solver == "davidson") { + for (int k = 0; k < solverOptions.sizeOfTheMatrix; ++k) { + // Calculate the correction vector = r/(D-theta*I) + real denominator = + 1.0 / (eigenValues.GetLocal(j, 0) - A.GetLocal(k, k)); + correctionVector.SetLocal(k, j, residual.GetLocal(k, j) * denominator); + } + } /*else if (solverOptions.solver == "jacobi") { + El::DistMatrix proj(grid); // projector matrix + El::Zeros(proj, solverOptions.sizeOfTheMatrix, + solverOptions.sizeOfTheMatrix); + El::DistMatrix 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 projProd(grid); // product of the projectors + El::Zeros(projProd, solverOptions.sizeOfTheMatrix, + solverOptions.sizeOfTheMatrix); + + El::DistMatrix projTemp(grid); // temp intermediate step + El::Zeros(projTemp, solverOptions.sizeOfTheMatrix, + solverOptions.sizeOfTheMatrix); + + El::Gemm(El::NORMAL, El::TRANSPOSE, alpha, proj, Atemp, beta, projTemp); + El::Gemm(El::NORMAL, El::TRANSPOSE, alpha, projTemp, proj, beta, + projProd); + + correctionVector = residual; + correctionVector *= -1.0; + + El::LinearSolve(projProd, correctionVector); + }*/ + } +} + +// Calculate the residual +template +bool eigenSolver::calculateResidual(int iterations, El::Grid &grid) { + // Parameters for GEMM + real alpha = 1, beta = 0; + + // Initialize the ritz vectors + El::Zeros(ritzVectors, solverOptions.sizeOfTheMatrix, + iterations * solverOptions.numberOfEigenValues); + + // Set ranges of V + El::Range begV(0, solverOptions.sizeOfTheMatrix); + El::Range endV(0, iterations * solverOptions.numberOfEigenValues); + + // Calculate the ritz vectors + El::Gemm(El::NORMAL, El::NORMAL, alpha, V(begV, endV), eigenVectors, beta, + ritzVectors); + + int residualCheck = 0; + for (int j = 0; j < solverOptions.numberOfEigenValues; ++j) { + + // Set ranges to get jth eigen vectors + El::Range beg(0, iterations * solverOptions.numberOfEigenValues); + El::Range end(j, j + 1); + + // Set ranges to get jth ritz vectors + El::Range begR(0, solverOptions.sizeOfTheMatrix); + El::Range endR(j, j + 1); + + // Matrix to store AV*s_j, where s is the eigenVector of the reduced problem + // and Vs the ritz vectors + El::DistMatrix AVs(grid); + El::Zeros(AVs, solverOptions.sizeOfTheMatrix, 1); + + // Calculate AV*s_j + El::Gemm(El::NORMAL, El::NORMAL, alpha, AV, eigenVectors(beg, end), beta, + AVs); + + // Matrix to store theta*V*s_j + El::DistMatrix thetaVs(ritzVectors(begR, endR)); + + // Calculate theta_j*I*Vs_j + thetaVs *= eigenValues.GetLocal(j, 0); + + // Calculate residual = AVs - thetaVs + AVs -= thetaVs; + double normAVs = El::Nrm2(AVs); + if (normAVs < solverOptions.tolerence) { + ++residualCheck; + } + + // Add the residual vector to the residual matrix + for (int k = 0; k < solverOptions.sizeOfTheMatrix; ++k) { + residual.SetLocal(k, j, AVs.GetLocal(k, 0)); + } + } + + if (residualCheck == solverOptions.numberOfEigenValues) + return true; + else + return false; +} + +// Calculate the eigen pairs of the V^TAV problem +template +void eigenSolver::subspaceProblem(int iterations, + const El::DistMatrix &A, + El::Grid &grid) { + // Set ranges of V + El::Range begV(0, solverOptions.sizeOfTheMatrix); + El::Range endV(0, iterations * solverOptions.numberOfEigenValues); + + // Initialise T = V^TAV + El::DistMatrix T(grid); + El::Zeros(T, iterations * solverOptions.numberOfEigenValues, + iterations * solverOptions.numberOfEigenValues); + + // Initialize AV + El::Identity(AV, solverOptions.sizeOfTheMatrix, + iterations * solverOptions.numberOfEigenValues); + + // Parameters for GEMM + real alpha = 1, beta = 0; + + // Calculate AV + El::Gemm(El::NORMAL, El::NORMAL, alpha, A, V(begV, endV), beta, AV); + + // Calculate V^TAV + El::Gemm(El::TRANSPOSE, El::NORMAL, alpha, V(begV, endV), AV, beta, T); + + // Get the eigen pairs for the reduced problem V^TAV + El::HermitianEig(El::UPPER, T, eigenValues, eigenVectors); +} + +// Expand the search space with the correction vector +template +void eigenSolver::expandSearchSpace(int &iterations, + const El::DistMatrix &A, + El::Grid &grid) { + //Implement Restart part + if ((iterations * solverOptions.numberOfEigenValues) >= sizeOfSearchSpace) { + iterations = 0; + eigenSolver::initialise(grid, A); + for (int j = 0; j < solverOptions.numberOfEigenValues; ++j) { + // Append the correction vector to the search space + for (int k = 0; k < solverOptions.sizeOfTheMatrix; ++k) { + V.SetLocal(k, j, correctionVector.GetLocal(k, j)); + } + } + } + //If no restart add the correction vector to V + else { + for (int j = 0; j < solverOptions.numberOfEigenValues; ++j) { + // Append the correction vector to the search space + for (int k = 0; k < solverOptions.sizeOfTheMatrix; ++k) { + V.SetLocal(k, (iterations)*solverOptions.numberOfEigenValues + j, + correctionVector.GetLocal(k, j)); + } + } + } +} // namespace eigenValueSolver + +template +void eigenSolver::solve(const El::DistMatrix &A, El::Grid &grid) { + + // Set grid sizes of the member matrices and initialize them + eigenSolver::initialise(grid, A); + int count = 0; + + for (int iterations = 1; iterations < solverOptions.sizeOfTheMatrix / 2; + ++iterations) { + + // Solve the subspace problem VTAV + subspaceProblem(iterations, A, grid); + + // Calculate the residual + bool tolerenceReached = calculateResidual(iterations, grid); + if (tolerenceReached == true) return; + + // Calculate the correction vector based on the solver + calculateCorrectionVector(iterations, A, grid); + + // calculate residual, correction vector and expand the search space + expandSearchSpace(iterations, A, grid); + + // Set ranges of V + El::Range begV(0, solverOptions.sizeOfTheMatrix); + El::Range endV(0, + (iterations + 1) * solverOptions.numberOfEigenValues); + + // R matrix for QR factorization + El::DistMatrix R; + El::Zeros(R, (iterations + 1) * solverOptions.numberOfEigenValues, + solverOptions.sizeOfTheMatrix); + + El::DistMatrix Q(V(begV, endV)); + // QR factorization of V + El::qr::Explicit(Q, R); + ++count; + } +} + +// explicit instantiations +template class eigenSolver; +template class eigenSolver; +} // namespace eigenValueSolver diff --git a/src/eigenValueSolver.cpp b/src/eigenValueSolver.cpp deleted file mode 100644 index 4ebc33a..0000000 --- a/src/eigenValueSolver.cpp +++ /dev/null @@ -1,177 +0,0 @@ -// -// eigenValueSolver.cpp -// eigenValueSolver -// -// Created by Adithya Vijaykumar on 17/09/2019. -// Copyright © 2019 Adithya Vijaykumar. All rights reserved. -// - -#include "eigenValueSolver.hpp" - -namespace eigenValueSolver { -template void eigenSolver::initialise(El::Grid &grid) { - - searchSpace.SetGrid(grid); - searchSpacesub.SetGrid(grid); - correctionVector.SetGrid(grid); - eigenVectors.SetGrid(grid); - eigenValues.SetGrid(grid); - eigenValues_old.SetGrid(grid); - - begTheta = El::Range{0, solverOptions.numberOfEigenValues}; - endTheta = El::Range{0, 1}; - - El::Identity(searchSpace, solverOptions.sizeOfTheMatrix, - solverOptions.sizeOfTheMatrix); - El::Identity(searchSpacesub, solverOptions.sizeOfTheMatrix, - solverOptions.sizeOfTheMatrix); - El::Identity(correctionVector, solverOptions.sizeOfTheMatrix, - columnsOfSearchSpace); - El::Identity(eigenVectors, solverOptions.sizeOfTheMatrix, - solverOptions.sizeOfTheMatrix); -} - -template -void eigenSolver::subspaceProblem(int iterations, El::DistMatrix &A, - El::Grid &grid) { - searchSpacesub = searchSpace; - searchSpacesub.Resize(solverOptions.sizeOfTheMatrix, iterations + 1); - - // Ttemp is AV and T is V^TAV i.e V^T.Ttemp - El::DistMatrix Ttemp(grid), T(grid); - El::Zeros(Ttemp, solverOptions.sizeOfTheMatrix, iterations + 1); - El::Zeros(T, iterations + 1, iterations + 1); - - // Parameters for GEMM - real alpha = 1, beta = 0; - - // AV - El::Gemm(El::NORMAL, El::NORMAL, alpha, A, searchSpacesub, beta, Ttemp); - // V^TAV - El::Gemm(El::TRANSPOSE, El::NORMAL, alpha, searchSpacesub, Ttemp, beta, T); - - // Get the eigen pairs for the reduced problem V^TAV - El::HermitianEig(El::UPPER, T, eigenValues, eigenVectors); -} - -template -void eigenSolver::expandSearchSpace(int iterations, - El::DistMatrix &A, - El::Grid &grid) { - - for (int j = 0; j < columnsOfSearchSpace; ++j) { - El::Range beg(0, iterations + 1); - El::Range end(j, j + 1); - - El::DistMatrix residual(grid); // residual Ay-thetay - El::Zeros(residual, solverOptions.sizeOfTheMatrix, 1); - - // calculate the ritz vector Vs - El::DistMatrix Vs(grid); - El::Zeros(Vs, solverOptions.sizeOfTheMatrix, 1); - real alpha = 1, beta = 0; - El::Gemv(El::NORMAL, alpha, searchSpacesub, eigenVectors(beg, end), beta, - Vs); - - El::DistMatrix I(grid); // Identitiy matrix - El::Identity( - I, solverOptions.sizeOfTheMatrix, - solverOptions.sizeOfTheMatrix); // Initialize as identity matrix - I *= eigenValues.GetLocal(j, 0); - El::DistMatrix Atemp(A); - 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 *= den; - } else if (solverOptions.solver == "jacobi") { - El::DistMatrix proj(grid); // projector matrix - El::Zeros(proj, solverOptions.sizeOfTheMatrix, - solverOptions.sizeOfTheMatrix); - El::DistMatrix 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 projProd(grid); // product of the projectors - El::Zeros(projProd, solverOptions.sizeOfTheMatrix, - solverOptions.sizeOfTheMatrix); - - El::DistMatrix projTemp(grid); // temp intermediate step - El::Zeros(projTemp, solverOptions.sizeOfTheMatrix, - solverOptions.sizeOfTheMatrix); - - El::Gemm(El::NORMAL, El::TRANSPOSE, alpha, proj, Atemp, beta, projTemp); - El::Gemm(El::NORMAL, El::TRANSPOSE, alpha, projTemp, proj, beta, - projProd); - - correctionVector = residual; - correctionVector *= -1.0; - - El::LinearSolve(projProd, correctionVector); - } - - for (int k = 0; k < solverOptions.sizeOfTheMatrix; ++k) { - searchSpace.SetLocal(k, iterations + j + 1, - correctionVector.GetLocal(k, 0)); - } - eigenValues_old -= eigenValues(begTheta, endTheta); - - if (El::Nrm2(eigenValues_old) < solverOptions.tolerence) { - break; - } - } -} - -template -void eigenSolver::solve(El::DistMatrix &A, El::Grid &grid) { - - eigenSolver::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 - { - - for (int i = 0; i < solverOptions.sizeOfTheMatrix; ++i) { - for (int j = 0; j < columnsOfSearchSpace; ++j) { - searchSpace.SetLocal(i, j, correctionVector.GetLocal(i, j)); - } - } - 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 - { - eigenValues_old = eigenValues(begTheta, endTheta); - } - - // Orthogonalize the searchSpace matrix using QR - El::DistMatrix R; // R matrix for QR factorization - El::Zeros(R, solverOptions.sizeOfTheMatrix, solverOptions.sizeOfTheMatrix); - - // QR factorization of V - El::qr::Explicit(searchSpace, R, false); - - // solve the subspace problem VTAV - subspaceProblem(iterations, A, grid); - // expand the search space - expandSearchSpace(iterations, A, grid); - } - El::Print(eigenValues(begTheta, endTheta)); -} - -// explicit instantiations -template class eigenSolver; -template class eigenSolver; -} // namespace eigenValueSolver diff --git a/src/main.cpp b/src/main.cc similarity index 66% rename from src/main.cpp rename to src/main.cc index 913348b..7909697 100644 --- a/src/main.cpp +++ b/src/main.cc @@ -8,6 +8,8 @@ #include "eigenValueSolver.hpp" #include "utils.hpp" +#include +#include using namespace eigenValueSolver; using namespace El; @@ -24,48 +26,71 @@ int main(int argc, char *argv[]) { try { - // initalise all the MPI stuff + // Initalise all the MPI variables + // Size of the block const El::Int blocksize = El::Input("--blocksize", "algorithmic blocksize", 128); + + // Height of the grid El::Int gridHeight = El::Input("--gridHeight", "grid height", 0); - // const El::Int numEigVal = - // El::Input("--numeigval", "size of matrix", - // 10); // Number of eigenvalues to be evaluated + + // Number of right hand sides const El::Int numRhs = El::Input("--numRhs", "# of right-hand sides", 1); + + // Error check const bool error = El::Input("--error", "test Elemental error?", true); + + // Print details const bool details = El::Input("--details", "print norm details?", true); - const El::Int matrixSize = El::Input("--size", "size of matrix", 100); + + // Size of the input matrix + const El::Int matrixSize = El::Input("--size", "size of matrix", 10); + + // The number of eigenvalues to be calculated const El::Int numEig = El::Input("--numeig", "number of eigenvalues", 1); + + // Type of the solver used 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 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); El::Grid grid{comm, gridHeight}; if (commRank == 0) El::Output("Grid is: ", grid.Height(), " x ", grid.Width()); - // Set up random A and B, then make the copies X := B - El::Timer timer; - // The matrix A whose eigenvalues have to be evaluated El::DistMatrix A(grid); + + // Initialize the matrix with zeros El::Zeros(A, matrixSize, matrixSize); + // Generate the Diagonally dominant hermitian matrix generateDDHermitianMatrix(A); + // Create an instance of the eigenSolver class eigenSolver solver; + + // Set solver options solver.solverOptions.numberOfEigenValues = numEig; solver.solverOptions.tolerence = 1e-8; solver.solverOptions.solver = solverType; solver.solverOptions.sizeOfTheMatrix = A.Height(); - El::Output(grid.Height()); + + // Solve function which calculates the eigenvalues and eigenvectors for + // matrix A solver.solve(A, grid); + + // Print eigenvalues + for (int i = 0; i < numEig; ++i) { + std::cout << std::setprecision( + std::numeric_limits::digits10 + 1) + << solver.eigenValues.GetLocal(i, 0) << std::endl; + } } catch (std::exception &e) { El::ReportException(e); } diff --git a/src/tests/test_eigen.cc b/src/tests/test_eigen.cc new file mode 100644 index 0000000..f56fec8 --- /dev/null +++ b/src/tests/test_eigen.cc @@ -0,0 +1,58 @@ +#define BOOST_TEST_MODULE eigen_solver +#include "eigenValueSolver.hpp" +#include "utils.hpp" +#include + +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 A(grid); + El::DistMatrix Ax(grid); + El::DistMatrix lambdax(grid); + El::Zeros(A, matrixSize, matrixSize); + El::Zeros(Ax, matrixSize, 1); + El::Zeros(lambdax, matrixSize, 1); + + // Generate the Diagonally dominant hermitian matrix + generateDDHermitianMatrix(A); + //El::Print(A); + eigenSolver 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 beg(0, solver.solverOptions.sizeOfTheMatrix); + El::Range end(0, 1); + + double eVal = solver.eigenValues.GetLocal(0,0); + El::DistMatrix eVec = solver.ritzVectors(beg,end); + + El::Gemm(El::NORMAL, El::NORMAL, 1.0, A, solver.ritzVectors(beg,end), 0.0, Ax); + El::Scale(eVal, eVec); + + + El::DistMatrix r = Ax; + r -= eVec; + double rNorm = El::Nrm2(r); + + BOOST_CHECK(rNorm < 1e-6); +} diff --git a/src/utils.cpp b/src/utils.cc similarity index 92% rename from src/utils.cpp rename to src/utils.cc index 9802db1..be434a1 100644 --- a/src/utils.cpp +++ b/src/utils.cc @@ -17,7 +17,7 @@ void generateDDHermitianMatrix(El::DistMatrix &A) { if (i == j) A.SetLocal(i, j, i + 21); else - A.SetLocal(i, j, matrixSize * 0.0000001); + A.SetLocal(i, j, 0.1); } } }