diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index 85c4767f..4c1b8d10 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -81,7 +81,7 @@ NAMESPACE_CHECK(PowerMethod); NAMESPACE_CHECK(multigrid); #include -#include #include +#include #endif diff --git a/Grid/algorithms/iterative/Arnoldi.h b/Grid/algorithms/iterative/Arnoldi.h index 19b5f6c7..6232681d 100644 --- a/Grid/algorithms/iterative/Arnoldi.h +++ b/Grid/algorithms/iterative/Arnoldi.h @@ -8,8 +8,6 @@ Copyright (C) 2015 Author: Peter Boyle Author: paboyle -Author: Chulwoo Jung -Author: Christoph Lehner Author: Patrick Oare This program is free software; you can redistribute it and/or modify @@ -34,37 +32,6 @@ See the full license in the file "LICENSE" in the top level distribution directo NAMESPACE_BEGIN(Grid); -/** - * Options for which Ritz values to keep in implicit restart. - */ -enum RitzFilter { - EvalNormSmall, // Keep evals with smallest norm - EvalNormLarge, // Keep evals with largest norm - EvalReSmall, // Keep evals with smallest real part - EvalReLarge // Keep evals with largest real part -}; - -// Select comparison function from RitzFilter -struct ComplexComparator -{ - RitzFilter f; - ComplexComparator (RitzFilter _f) : f(_f) {} - bool operator()(std::complex z1, std::complex z2) { - switch (f) { - case EvalNormSmall: - return std::abs(z1) < std::abs(z2); - case EvalNormLarge: - return std::abs(z1) > std::abs(z2); - case EvalReSmall: - return std::abs(std::real(z1)) < std::abs(std::real(z2)); - case EvalReLarge: - return std::abs(std::real(z1)) > std::abs(std::real(z2)); - default: - assert(0); - } - } -}; - /** * Implementation of the Arnoldi algorithm. */ diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 9906741b..a83e7014 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -8,8 +8,6 @@ Copyright (C) 2015 Author: Peter Boyle Author: paboyle -Author: Chulwoo Jung -Author: Christoph Lehner Author: Patrick Oare This program is free software; you can redistribute it and/or modify @@ -32,14 +30,76 @@ See the full license in the file "LICENSE" in the top level distribution directo #ifndef GRID_KRYLOVSCHUR_H #define GRID_KRYLOVSCHUR_H -// #include -// #include -// #include - -// #include - NAMESPACE_BEGIN(Grid); +/** + * Options for which Ritz values to keep in implicit restart. TODO move this and utilities into a new file + */ +enum RitzFilter { + EvalNormSmall, // Keep evals with smallest norm + EvalNormLarge, // Keep evals with largest norm + EvalReSmall, // Keep evals with smallest real part + EvalReLarge, // Keep evals with largest real part + EvalImSmall, // Keep evals with smallest imaginary part + EvalImLarge // Keep evals with largest imaginary part +}; + +/** Selects the RitzFilter corresponding to the input string. */ +RitzFilter selectRitzFilter(std::string s) { + if (s == "EvalNormSmall") { return EvalNormSmall; } else + if (s == "EvalNormLarge") { return EvalNormLarge; } else + if (s == "EvalReSmall") { return EvalReSmall; } else + if (s == "EvalReLarge") { return EvalReLarge; } else + if (s == "EvalImSmall") { return EvalImSmall; } else + if (s == "EvalImLarge") { return EvalImLarge; } else + { assert(0); } +} + +/** Returns a string saying which RitzFilter it is. */ +std::string rfToString(RitzFilter RF) { + switch (RF) { + case EvalNormSmall: + return "EvalNormSmall"; + case EvalNormLarge: + return "EvalNormLarge"; + case EvalReSmall: + return "EvalReSmall"; + case EvalReLarge: + return "EvalReLarge"; + case EvalImSmall: + return "EvalImSmall"; + case EvalImLarge: + return "EvalImLarge"; + default: + assert(0); + } +} + +// Select comparison function from RitzFilter +struct ComplexComparator +{ + RitzFilter RF; + ComplexComparator (RitzFilter _rf) : RF(_rf) {} + bool operator()(std::complex z1, std::complex z2) { + switch (RF) { + case EvalNormSmall: + return std::abs(z1) < std::abs(z2); + case EvalNormLarge: + return std::abs(z1) > std::abs(z2); + case EvalReSmall: + return std::abs(std::real(z1)) < std::abs(std::real(z2)); + case EvalReLarge: + return std::abs(std::real(z1)) > std::abs(std::real(z2)); + case EvalImSmall: + return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); + case EvalImLarge: + return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); + default: + assert(0); + } + } +}; + /** * Computes a complex Schur decomposition of a complex matrix A using Eigen's matrix library. The Schur decomposition, * A = Q^\dag S Q diff --git a/examples/Example_spec_kryschur.cc b/examples/Example_spec_kryschur.cc index 61c6ab30..aad64dbd 100644 --- a/examples/Example_spec_kryschur.cc +++ b/examples/Example_spec_kryschur.cc @@ -350,6 +350,15 @@ int main (int argc, char ** argv) std::string file = argv[5]; std::string outDir = argv[6]; + RitzFilter RF; + if (argc == 8) { + std::string rf = argv[7]; + RF = selectRitzFilter(rf); + } else { + RF = EvalReSmall; + } + std::cout << "Sorting eigenvalues using " << rfToString(RF) << std::endl; + //const int Ls=16; const int Ls = 8; @@ -438,7 +447,7 @@ int main (int argc, char ** argv) << ", Nstop = " << Nstop << std::endl; // Arnoldi Arn(PVdagM, FGrid, 1e-8); // Arn(src, maxIter, Nm, Nk, Nstop); - KrylovSchur KrySchur (PVdagM, FGrid, 1e-8); + KrylovSchur KrySchur (PVdagM, FGrid, 1e-8, RF); KrySchur(src, maxIter, Nm, Nk, Nstop); std::cout<