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 b4a9b8b8..a1cb2598 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,7 +32,10 @@ See the full license in the file "LICENSE" in the top level distribution directo NAMESPACE_BEGIN(Grid); +//Moved to KrylovSchur +#if 0 /** +<<<<<<< HEAD * Options for which Ritz values to keep in implicit restart. */ enum RitzFilter { @@ -74,6 +75,10 @@ struct ComplexComparator } }; +======= +>>>>>>> 68af1bba67dd62881ead5ab1e54962a5486a0791 +#endif + /** * Implementation of the Arnoldi algorithm. */ diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 5441af1a..d5bf6a20 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 @@ -34,6 +32,86 @@ 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. 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 + EvalImNormSmall, // Keep evals with smallest |imaginary| part + EvalImNormLarge, // Keep evals with largest |imaginary| part +}; + +/** Selects the RitzFilter corresponding to the input string. */ +inline 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 + if (s == "EvalImNormSmall") { return EvalImNormSmall; } else + if (s == "EvalImNormLarge") { return EvalImNormLarge; } else + { assert(0); } +} + +/** Returns a string saying which RitzFilter it is. */ +inline 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"; + case EvalImNormSmall: + return "EvalImNormSmall"; + case EvalImNormLarge: + return "EvalImNormLarge"; + 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::real(z1) < std::real(z2); // DELETE THE ABS HERE!!! + case EvalReLarge: + return std::real(z1) > std::real(z2); + case EvalImSmall: + return std::imag(z1) < std::imag(z2); + case EvalImLarge: + return std::imag(z1) > std::imag(z2); + case EvalImNormSmall: + return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); + case EvalImNormLarge: + 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 @@ -172,6 +250,15 @@ class ComplexSchurDecomposition { }; +// template +// inline void writeFile(const Field &field, const std::string &fname) { +// emptyUserRecord record; +// ScidacWriter WR(field.Grid()->IsBoss()); +// WR.open(fname); +// WR.writeScidacFieldRecord(field, record, 0); // 0 = Lexico +// WR.close(); +// } + /** * Implementation of the Krylov-Schur algorithm. */ @@ -195,21 +282,22 @@ class KrylovSchur { RealD approxLambdaMax; RealD beta_k; Field u; // Residual vector perpendicular to Krylov space (u_{k+1} in notes) - Eigen::VectorXcd b; // b vector in Schur decomposition (e_{k+1} in Arnoldi). + Eigen::VectorXcd b; // b vector in Schur decomposition (e_{k+1} in Arnoldi). std::vector basis; // orthonormal Krylov basis - Eigen::MatrixXcd Rayleigh; // Rayleigh quotient of size Nbasis (after construction) - Eigen::MatrixXcd Qt; // Transpose of basis rotation which projects out high modes. + Eigen::MatrixXcd Rayleigh; // Rayleigh quotient of size Nbasis (after construction) + Eigen::MatrixXcd Qt; // Transpose of basis rotation which projects out high modes. - Eigen::VectorXcd evals; // evals of Rayleigh quotient - Eigen::MatrixXcd littleEvecs; // Nm x Nm evecs matrix + Eigen::VectorXcd evals; // evals of Rayleigh quotient + std::vector ritzEstimates; // corresponding ritz estimates for evals + Eigen::MatrixXcd littleEvecs; // Nm x Nm evecs matrix - RitzFilter ritzFilter; // how to sort evals + RitzFilter ritzFilter; // how to sort evals public: KrylovSchur(LinearOperatorBase &_Linop, GridBase *_Grid, RealD _Tolerance, RitzFilter filter = EvalReSmall) : Linop(_Linop), Grid(_Grid), Tolerance(_Tolerance), ritzFilter(filter), u(_Grid), MaxIter(-1), Nm(-1), Nk(-1), Nstop (-1), - evals (0), evecs (), ssq (0.0), rtol (0.0), beta_k (0.0), approxLambdaMax (0.0) + evals (0), ritzEstimates (), evecs (), ssq (0.0), rtol (0.0), beta_k (0.0), approxLambdaMax (0.0) { u = Zero(); }; @@ -218,11 +306,13 @@ class KrylovSchur { /* Getters */ - Eigen::MatrixXcd getRayleighQuotient() { return Rayleigh; } - Field getU() { return u; } - std::vector getBasis() { return basis; } - Eigen::VectorXcd getEvals() { return evals; } - std::vector getEvecs() { return evecs; } + int getNk() { return Nk; } + Eigen::MatrixXcd getRayleighQuotient() { return Rayleigh; } + Field getU() { return u; } + std::vector getBasis() { return basis; } + Eigen::VectorXcd getEvals() { return evals; } + std::vector getRitzEstimates() { return ritzEstimates; } + std::vector getEvecs() { return evecs; } /** * Runs the Krylov-Schur loop. @@ -239,11 +329,15 @@ class KrylovSchur { ssq = norm2(v0); RealD approxLambdaMax = approxMaxEval(v0); rtol = Tolerance * approxLambdaMax; + std::cout << GridLogMessage << "Approximate max eigenvalue: " << approxLambdaMax << std::endl; // rtol = Tolerance; b = Eigen::VectorXcd::Zero(Nm); // start as e_{k+1} b(Nm-1) = 1.0; + // basis = new std::vector (Nm, Grid); + // evecs.reserve(); + int start = 0; Field startVec = v0; littleEvecs = Eigen::MatrixXcd::Zero(Nm, Nm); @@ -281,6 +375,10 @@ class KrylovSchur { // basisRotate(evecs, Q, 0, Nm, 0, Nm, Nm); std::vector basis2; + // basis2.reserve(Nm); + // for (int i = start; i < Nm; i++) { + // basis2.emplace_back(Grid); + // } constructUR(basis2, basis, Qt, Nm); basis = basis2; @@ -298,18 +396,12 @@ class KrylovSchur { std::cout << GridLogDebug << "Rayleigh before truncation: " << std::endl << Rayleigh << std::endl; - // Rayleigh = Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); Eigen::MatrixXcd RayTmp = Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); Rayleigh = RayTmp; - // basis = std::vector (basis.begin(), basis.begin() + Nk); std::vector basisTmp = std::vector (basis.begin(), basis.begin() + Nk); basis = basisTmp; - // evecs = std::vector (evecs.begin(), evecs.begin() + Nk); - // littleEvecs = littleEvecs(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); - - // b = b.head(Nk); Eigen::VectorXcd btmp = b.head(Nk); b = btmp; @@ -327,8 +419,11 @@ class KrylovSchur { if (Nconv >= Nstop || i == MaxIter - 1) { std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration " << i << "." << std::endl; - basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm); + // basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm); // Think this might have been the issue std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl; + + // writeEigensystem(path); + return; } } @@ -355,15 +450,27 @@ class KrylovSchur { ComplexD coeff; Field w (Grid); // A acting on last Krylov vector. + // basis.reserve(Nm); + // for (int i = start; i < Nm; i++) { + // basis.emplace_back(Grid); + // } + // basis.assign(Nm, Field(Grid)); + // basis.resize(Nm); + // for (int i = start; i < Nm; i++) { + // basis[i] = Field(Grid); + // } + if (start == 0) { // initialize everything that we need. RealD v0Norm = 1 / std::sqrt(ssq); basis.push_back(v0Norm * v0); // normalized source + // basis[0] = v0Norm * v0; // normalized source Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm); u = Zero(); } else { std::cout << GridLogMessage << "start= " < tri; @@ -447,7 +558,7 @@ class KrylovSchur { std::cout << GridLogDebug << "Little evecs: " << littleEvecs << std::endl; // Convert evecs to lattice fields - for (int k = 0; k < evals.size(); k++) { + for (int k = 0; k < n; k++) { Eigen::VectorXcd vec = littleEvecs.col(k); Field tmp (basis[0].Grid()); tmp = Zero(); @@ -455,6 +566,7 @@ class KrylovSchur { tmp = tmp + vec[j] * basis[j]; } evecs.push_back(tmp); + // evecs[k] = tmp; } } @@ -510,11 +622,16 @@ class KrylovSchur { */ int converged() { int Nconv = 0; + int _Nm = evecs.size(); std::cout << GridLogDebug << "b: " << b << std::endl; Field tmp (Grid); Field fullEvec (Grid); - for (int k = 0; k < evecs.size(); k++) { + ritzEstimates.clear(); + // ritzEstimates.resize(_Nm); + for (int k = 0; k < _Nm; k++) { Eigen::VectorXcd evec_k = littleEvecs.col(k); RealD ritzEstimate = std::abs(b.dot(evec_k)); // b^\dagger s + ritzEstimates.push_back(ritzEstimate); + // ritzEstimates[k] = ritzEstimate; std::cout << GridLogMessage << "Ritz estimate for evec " << k << " = " << ritzEstimate << std::endl; if (ritzEstimate < rtol) { Nconv++; @@ -522,7 +639,7 @@ class KrylovSchur { } // Check that Ritz estimate is explicitly || D (Uy) - lambda (Uy) || - // checkRitzEstimate(); + checkRitzEstimate(); return Nconv; } @@ -539,7 +656,7 @@ class KrylovSchur { // rotate basis by Rayleigh to construct UR // std::vector rotated; - std::cout << GridLogDebug << "Rayleigh in KSDecomposition: " << std::endl << Rayleigh << std::endl; + // std::cout << GridLogDebug << "Rayleigh in KSDecomposition: " << std::endl << Rayleigh << std::endl; std::vector rotated = basis; constructUR(rotated, basis, Rayleigh, k); // manually rotate @@ -558,10 +675,10 @@ class KrylovSchur { delta = delta / norm2(tmp); // relative tolerance deltaSum += delta; - std::cout << GridLogDebug << "Iteration " << i << std::endl; - std::cout << GridLogDebug << "Du = " << norm2(tmp) << std::endl; - std::cout << GridLogDebug << "rotated = " << norm2(rotated[i]) << std::endl; - std::cout << GridLogDebug << "b[i] = " << b(i) << std::endl; + // std::cout << GridLogDebug << "Iteration " << i << std::endl; + // std::cout << GridLogDebug << "Du = " << norm2(tmp) << std::endl; + // std::cout << GridLogDebug << "rotated = " << norm2(rotated[i]) << std::endl; + // std::cout << GridLogDebug << "b[i] = " << b(i) << std::endl; std::cout << GridLogMessage << "Deviation in decomp, column " << i << ": " << delta << std::endl; } std::cout << GridLogMessage << "Squared sum of relative deviations in decomposition: " << deltaSum << std::endl; @@ -629,7 +746,9 @@ class KrylovSchur { */ void constructUR(std::vector& UR, std::vector &U, Eigen::MatrixXcd& R, int N) { Field tmp (Grid); + UR.clear(); + // UR.resize(N); std::cout << GridLogDebug << "R to rotate by (should be Rayleigh): " << R << std::endl; @@ -642,6 +761,7 @@ class KrylovSchur { } std::cout << GridLogDebug << "rotated norm at i = " << i << " is: " << norm2(tmp) << std::endl; UR.push_back(tmp); + // UR[i] = tmp; } return; } @@ -652,17 +772,61 @@ class KrylovSchur { void constructRU(std::vector& RU, std::vector &U, Eigen::MatrixXcd& R, int N) { Field tmp (Grid); RU.clear(); + // RU.resize(N); for (int i = 0; i < N; i++) { tmp = Zero(); for (int j = 0; j < N; j++) { tmp = tmp + R(i, j) * U[j]; } RU.push_back(tmp); + // RU[i] = tmp; } return; } + // void writeEvec(Field& in, std::string const fname){ + // #ifdef HAVE_LIME + // // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111 + // std::cout << GridLogMessage << "Writing evec to: " << fname << std::endl; + // Grid::emptyUserRecord record; + // Grid::ScidacWriter WR(in.Grid()->IsBoss()); + // WR.open(fname); + // WR.writeScidacFieldRecord(in,record,0); // Lexico + // WR.close(); + // #endif + // // What is the appropriate way to throw error? + // } + + // /** + // * Writes the eigensystem of a Krylov Schur object to a directory. + // * + // * Parameters + // * ---------- + // * std::string path + // * Directory to write to. + // */ + // void writeEigensystem(std::string outDir) { + // std::cout << GridLogMessage << "Writing output to directory: " << outDir << std::endl; + // // TODO write a scidac density file so that we can easily integrate with visualization toolkit + // std::string evalPath = outDir + "/evals.txt"; + // std::ofstream fEval; + // fEval.open(evalPath); + // for (int i = 0; i < Nk; i++) { + // // write Eigenvalues + // fEval << i << " " << evals(i); + // if (i < Nk - 1) { fEval << "\n"; } + // } + // fEval.close(); + + // for (int i = 0; i < Nk; i++) { + // std::string fName = outDir + "/evec" + std::to_string(i); + // // writeFile(evecs[i], fName); // using method from Grid/HMC/ComputeWilsonFlow.cc + // // writeEvec(evecs[i], fName); + // } + + // } + }; - + NAMESPACE_END(Grid); #endif diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 71a6cf91..6721eaa6 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -536,7 +536,7 @@ int main (int argc, char ** argv) // KrylovSchur KrySchur (Dsq, FGrid, 1e-8, EvalNormLarge); // KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall); // Hacked, really EvalImagSmall - KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalReSmall); + KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); // BlockKrylovSchur KrySchur (HermOp2, UGrid, Nu, resid,EvalNormSmall); KrySchur(src[0], maxIter, Nm, Nk, Nstop); std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl; diff --git a/examples/Example_spec_kryschur.cc b/examples/Example_spec_kryschur.cc index 460881cb..7e70a180 100644 --- a/examples/Example_spec_kryschur.cc +++ b/examples/Example_spec_kryschur.cc @@ -1,12 +1,35 @@ /************************************************************************************* + Runs the Krylov-Schur algorithm on a (pre-conditioned) domain-wall fermion operator + to determine part of its spectrum. + + Usage : + $ ./Example_spec_kryschur + + Nm = Maximum size of approximation subspace. + Nk = Size of truncation subspace + maxiter = Maximum number of iterations. + Nstop = Stop when Nstop eigenvalues have converged. + inFile = Gauge configuration to read in. + outDir = Directory to write output to. + rf = (Optional) RitzFilter to sort with. Takes in any string in + {EvalNormSmall, EvalNormLarge, EvalReSmall, EvalReLarge, EvalImSmall, EvalImLarge} + + Output: + ${outDir}/evals.txt = Contains all eigenvalues. Each line is formatted as `$idx $eval $ritz`, where: + - $idx is the index of the eigenvalue. + - $eval is the eigenvalue, formated as "(re,im)". + - $ritz is the Ritz estimate of the eigenvalue (deviation from being a true eigenvalue) + ${outDir}/evec${idx} = Eigenvector $idx written out in SCIDAC format (if LIME is enabled). + Grid physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_padded_cell.cc Copyright (C) 2023 -Author: Peter Boyle + Author: Peter Boyle + Author: Patrick Oare This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -26,10 +49,6 @@ Author: Peter Boyle *************************************************************************************/ /* END LEGAL */ -// copied here from Test_general_coarse_pvdagm.cc - -// copied here from Test_general_coarse_pvdagm.cc - #include #include @@ -40,9 +59,13 @@ Author: Peter Boyle #include #include +#include +#include + using namespace std; using namespace Grid; +<<<<<<< HEAD namespace Grid { struct LanczosParameters: Serializable { @@ -87,6 +110,54 @@ struct LanczosParameters: Serializable { } #if 0 +======= +template void writeFile(T& in, std::string const fname){ + #ifdef HAVE_LIME + // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111 + std::cout << Grid::GridLogMessage << "Writes to: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter WR(in.Grid()->IsBoss()); + WR.open(fname); + WR.writeScidacFieldRecord(in,record,0); // Lexico + WR.close(); + #endif +} + +/** + * Writes the eigensystem of a Krylov Schur object to a directory. + * + * Parameters + * ---------- + * std::string path + * Directory to write to. + */ +template +void writeEigensystem(KrylovSchur KS, std::string outDir) { + int Nk = KS.getNk(); + std::cout << GridLogMessage << "Writing output to directory: " << outDir << std::endl; + + // Write evals + std::string evalPath = outDir + "/evals.txt"; + std::ofstream fEval; + fEval.open(evalPath); + Eigen::VectorXcd evals = KS.getEvals(); + std::vector ritz = KS.getRitzEstimates(); + for (int i = 0; i < Nk; i++) { + // write eigenvalues and Ritz estimates + fEval << i << " " << evals(i) << " " << ritz[i]; + if (i < Nk - 1) { fEval << "\n"; } + } + fEval.close(); + + // Write evecs (TODO: very heavy on storage costs! Don't write them all out) + // std::vector evecs = KS.getEvecs(); + // for (int i = 0; i < Nk; i++) { + // std::string fName = outDir + "/evec" + std::to_string(i); + // writeFile(evecs[i], fName); // using method from Grid/HMC/ComputeWilsonFlow.cc + // } +} + +>>>>>>> 68af1bba67dd62881ead5ab1e54962a5486a0791 // Hermitize a DWF operator by squaring it template class SquaredLinearOperator : public LinearOperatorBase { @@ -232,6 +303,7 @@ ShiftedComplexPVdagMLinearOperator(ComplexD _shift,Matrix &Mat,Matrix &PV): shif } }; +<<<<<<< HEAD template class MGPreconditioner : public LinearFunction< Lattice > { public: @@ -334,16 +406,28 @@ public: }; #endif +======= +>>>>>>> 68af1bba67dd62881ead5ab1e54962a5486a0791 int main (int argc, char ** argv) { Grid_init(&argc,&argv); - // Usage : $ ./Example_spec_kryschur - // assert (argc == 5); - std::string NmStr = argv[1]; - std::string NkStr = argv[2]; + // Usage : $ ./Example_spec_kryschur + std::string NmStr = argv[1]; + std::string NkStr = argv[2]; std::string maxIterStr = argv[3]; - std::string NstopStr = argv[4]; + std::string NstopStr = argv[4]; + 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; @@ -360,52 +444,24 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - // Construct a coarsened grid - // poare TODO: replace this with the following line? - Coordinate clatt = lat_size; -// Coordinate clatt = GridDefaultLatt(); // [PO] initial line before I edited it - for(int d=0;d seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); - std::vector cseeds({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); - GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds); LatticeFermion src(FGrid); random(RNG5,src); - LatticeFermion result(FGrid); result=Zero(); - LatticeFermion ref(FGrid); ref=Zero(); - LatticeFermion tmp(FGrid); - LatticeFermion err(FGrid); LatticeGaugeField Umu(UGrid); FieldMetaData header; - // std::string file ("/sdcc/u/poare/PETSc-Grid/ckpoint_EODWF_lat.125"); - std::string file("/Users/patrickoare/libraries/PETSc-Grid/ckpoint_EODWF_lat.125"); NerscIO::readConfiguration(Umu,header,file); - RealD mass=0.01; + // RealD mass=0.01; + RealD mass=0.001; RealD M5=1.8; DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5); - // const int nbasis = 20; // size of approximate basis for low-mode space - const int nbasis = 3; // size of approximate basis for low-mode space - const int cb = 0 ; - LatticeFermion prom(FGrid); - - typedef GeneralCoarsenedMatrix LittleDiracOperator; - typedef LittleDiracOperator::CoarseVector CoarseVector; - - NextToNearestStencilGeometry5D geom(Coarse5d); - std::cout< Dsq (Ddwf); NonHermitianLinearOperator DLinOp (Ddwf); - // int Nm = 200; - // int Nk = 110; - // int maxIter = 2000; - // int Nstop = 100; - int Nm = std::stoi(NmStr); int Nk = std::stoi(NkStr); int maxIter = std::stoi(maxIterStr); @@ -430,9 +481,9 @@ int main (int argc, char ** argv) std::cout << GridLogMessage << "Runnning Krylov Schur. Nm = " << Nm << ", Nk = " << Nk << ", maxIter = " << maxIter << ", 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); // use preconditioned PV^\dag D_{dwf} + // KrylovSchur KrySchur (DLinOp, FGrid, 1e-8, RF); // use D_{dwf} KrySchur(src, maxIter, Nm, Nk, Nstop); std::cout< + + Nm = Maximum size of approximation subspace. + Nk = Size of truncation subspace + maxiter = Maximum number of iterations. + Nstop = Stop when Nstop eigenvalues have converged. + inFile = Gauge configuration to read in. + outDir = Directory to write output to. + rf = (Optional) RitzFilter to sort with. Takes in any string in + {EvalNormSmall, EvalNormLarge, EvalReSmall, EvalReLarge, EvalImSmall, EvalImLarge} + + Output: + ${outDir}/evals.txt = Contains all eigenvalues. Each line is formatted as `$idx $eval $ritz`, where: + - $idx is the index of the eigenvalue. + - $eval is the eigenvalue, formated as "(re,im)". + - $ritz is the Ritz estimate of the eigenvalue (deviation from being a true eigenvalue) + ${outDir}/evec${idx} = Eigenvector $idx written out in SCIDAC format (if LIME is enabled). + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_padded_cell.cc + + Copyright (C) 2023 + + Author: Peter Boyle + Author: Patrick Oare + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +using namespace std; +using namespace Grid; + +template void writeFile(T& in, std::string const fname){ + #ifdef HAVE_LIME + // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111 + std::cout << Grid::GridLogMessage << "Writes to: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter WR(in.Grid()->IsBoss()); + WR.open(fname); + WR.writeScidacFieldRecord(in,record,0); // Lexico + WR.close(); + #endif +} + +template void readFile(T& out, std::string const fname){ + #ifdef HAVE_LIME + // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111 + std::cout << Grid::GridLogMessage << "Reads at: " << fname << std::endl; + Grid::emptyUserRecord record; + // Grid::ScidacReader SR(out.Grid()->IsBoss()); + Grid::ScidacReader SR; + SR.open(fname); + SR.readScidacFieldRecord(out, record); + SR.close(); + #endif +} + +/** + * Writes the eigensystem of a Krylov Schur object to a directory. + * + * Parameters + * ---------- + * std::string path + * Directory to write to. + */ +template +void writeEigensystem(KrylovSchur KS, std::string outDir) { + int Nk = KS.getNk(); + std::cout << GridLogMessage << "Writing output to directory: " << outDir << std::endl; + + // Write evals + std::string evalPath = outDir + "/evals.txt"; + std::ofstream fEval; + fEval.open(evalPath); + Eigen::VectorXcd evals = KS.getEvals(); + std::vector ritz = KS.getRitzEstimates(); + for (int i = 0; i < Nk; i++) { + // write eigenvalues and Ritz estimates + fEval << i << " " << evals(i) << " " << ritz[i]; + if (i < Nk - 1) { fEval << "\n"; } + } + fEval.close(); + + // Write evecs + int Nevecs = Nk; // don't write all of them + std::vector evecs = KS.getEvecs(); + for (int i = 0; i < Nevecs; i++) { + std::string fName = outDir + "/evec" + std::to_string(i); + writeFile(evecs[i], fName); // using method from Grid/HMC/ComputeWilsonFlow.cc + } +} + +// Hermitize a DWF operator by squaring it +template +class SquaredLinearOperator : public LinearOperatorBase { + + public: + Matrix &_Mat; + + public: + SquaredLinearOperator(Matrix &Mat): _Mat(Mat) {}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + // std::cout << "Op is overloaded as HermOp" << std::endl; + HermOp(in, out); + } + void AdjOp (const Field &in, Field &out){ + HermOp(in, out); + } + void _Op (const Field &in, Field &out){ + // std::cout << "Op: M "< +class PVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< +class ShiftedPVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; + RealD shift; +public: + ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< +class ShiftedComplexPVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; + ComplexD shift; +public: +ShiftedComplexPVdagMLinearOperator(ComplexD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< lat_size {16, 16, 16, 32}; + std::vector lat_size {32, 32, 32, 32}; + std::cout << "Lattice size: " << lat_size << std::endl; + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size, + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + // GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + // GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + GridCartesian * FGrid = UGrid; + GridRedBlackCartesian * FrbGrid = UrbGrid; + + std::vector seeds4({1,2,3,4}); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + LatticeFermion src(FGrid); random(RNG4, src); + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + NerscIO::readConfiguration(Umu, header, file); + + std::cout << GridLogMessage << "Loaded configuration" << std::endl; + + // RealD mass = 0.01; + RealD M5 = 1.8; + + // Wilson mass + RealD mass = -1.6; + + std::cout << GridLogMessage << "masses specified" << std::endl; + + std::vector boundary = {1,1,1,-1}; + WilsonFermionD::ImplParams Params(boundary); + + // DomainWallFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5); + // NonHermitianLinearOperator DLinOp (Ddwf); + + // WilsonFermionD Dwilson(Umu, *FGrid, *FrbGrid, mass); + WilsonFermionD Dwilson(Umu, *UGrid, *UrbGrid, mass, Params); + NonHermitianLinearOperator DLinOp (Dwilson); + + std::cout << GridLogMessage << "Dirac operator defined" << std::endl; + + std::string eigenPath = "/home/poare/lqcd/multigrid/spectra/32cube-rho0.124-tau4/U_smr_3.000000/Nm72_Nk24_8111835.aurora-pbs-0001.hostmgmt.cm.aurora.alcf.anl.gov/"; + + std::cout << GridLogMessage << "Loading eigenvalues" << std::endl; + std::ifstream evalFile(eigenPath + "evals.txt"); + std::string str; + std::vector evals; + while (std::getline(evalFile, str)) { + std::cout << GridLogMessage << "Reading line: " << str << std::endl; + int i1 = str.find("(") + 1; + int i2 = str.find(",") + 1; + int i3 = str.find(")"); + std::cout << "i1,i2,i3 = " << i1 << "," << i2 << "," << i3 << std::endl; + std::string reStr = str.substr(i1, i2 - i1); + std::string imStr = str.substr(i2, i3 - i2); + std::cout << GridLogMessage << "Parsed re = " << reStr << " and im = " << imStr << std::endl; + // ComplexD z (std::stof(reStr), std::stof(imStr)); + ComplexD z (std::stod(reStr), std::stod(imStr)); + evals.push_back(z); + } + std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl; + + int Nevecs = 24; + std::vector evecs; + LatticeFermion evec (FGrid); + for (int i = 0; i < Nevecs; i++) { + std::string evecPath = eigenPath + "evec" + std::to_string(i); + readFile(evec, evecPath); + evecs.push_back(evec); + } + std::cout << GridLogMessage << "Evecs loaded" << std::endl; + + // Compute < evec | D - \lambda | evec > + std::cout << GridLogMessage << "Testing eigenvectors" << std::endl; + LatticeFermion Devec (FGrid); + ComplexD ritz; + for (int i = 0; i < Nevecs; i++) { + Devec = Zero(); + DLinOp.Op(evecs[i], Devec); + ritz = std::sqrt(norm2(Devec - evals[i] * evecs[i])); + std::cout << GridLogMessage << "i = " << i << ", || (D - lambda) |vi> || = " << ritz << std::endl; + } + // Eigen::MatrixXcd Dw_evecs; + // Dw_evecs = Eigen::MatrixXcd::Zero(Nevecs, Nevecs); + // for (int i = 0; i < Nevecs; i++) { + // Linop.Op(evecs[i], Devec); + // for (int j = 0; j < Nevecs; j++) { + + // } + // } + + std::cout< + + Nm = Maximum size of approximation subspace. + Nk = Size of truncation subspace + maxiter = Maximum number of iterations. + Nstop = Stop when Nstop eigenvalues have converged. + inFile = Gauge configuration to read in. + outDir = Directory to write output to. + rf = (Optional) RitzFilter to sort with. Takes in any string in + {EvalNormSmall, EvalNormLarge, EvalReSmall, EvalReLarge, EvalImSmall, EvalImLarge} + + Output: + ${outDir}/evals.txt = Contains all eigenvalues. Each line is formatted as `$idx $eval $ritz`, where: + - $idx is the index of the eigenvalue. + - $eval is the eigenvalue, formated as "(re,im)". + - $ritz is the Ritz estimate of the eigenvalue (deviation from being a true eigenvalue) + ${outDir}/evec${idx} = Eigenvector $idx written out in SCIDAC format (if LIME is enabled). + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_padded_cell.cc + + Copyright (C) 2023 + + Author: Peter Boyle + Author: Patrick Oare + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +using namespace std; +using namespace Grid; + +template void writeFile(T& in, std::string const fname){ + #ifdef HAVE_LIME + // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111 + std::cout << Grid::GridLogMessage << "Writes to: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter WR(in.Grid()->IsBoss()); + WR.open(fname); + WR.writeScidacFieldRecord(in,record,0); // Lexico + WR.close(); + #endif +} + +/** + * Writes the eigensystem of a Krylov Schur object to a directory. + * + * Parameters + * ---------- + * std::string path + * Directory to write to. + */ +template +void writeEigensystem(KrylovSchur KS, std::string outDir) { + int Nk = KS.getNk(); + std::cout << GridLogMessage << "Writing output to directory: " << outDir << std::endl; + + // Write evals + std::string evalPath = outDir + "/evals.txt"; + std::ofstream fEval; + fEval.open(evalPath); + Eigen::VectorXcd evals = KS.getEvals(); + std::vector ritz = KS.getRitzEstimates(); + for (int i = 0; i < Nk; i++) { + // write eigenvalues and Ritz estimates + fEval << i << " " << evals(i) << " " << ritz[i]; + if (i < Nk - 1) { fEval << "\n"; } + } + fEval.close(); + + // Write evecs + int Nevecs = Nk; // don't write all of them + std::vector evecs = KS.getEvecs(); + for (int i = 0; i < Nevecs; i++) { + std::string fName = outDir + "/evec" + std::to_string(i); + writeFile(evecs[i], fName); // using method from Grid/HMC/ComputeWilsonFlow.cc + } +} + +// Hermitize a DWF operator by squaring it +template +class SquaredLinearOperator : public LinearOperatorBase { + + public: + Matrix &_Mat; + + public: + SquaredLinearOperator(Matrix &Mat): _Mat(Mat) {}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + // std::cout << "Op is overloaded as HermOp" << std::endl; + HermOp(in, out); + } + void AdjOp (const Field &in, Field &out){ + HermOp(in, out); + } + void _Op (const Field &in, Field &out){ + // std::cout << "Op: M "< +class PVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< +class ShiftedPVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; + RealD shift; +public: + ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< +class ShiftedComplexPVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; + ComplexD shift; +public: +ShiftedComplexPVdagMLinearOperator(ComplexD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< + std::string NmStr = argv[1]; + std::string NkStr = argv[2]; + std::string maxIterStr = argv[3]; + std::string NstopStr = argv[4]; + 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; + // } + // RitzFilter RF; + std::string rf = argv[7]; + RitzFilter RF = selectRitzFilter(rf); + std::cout << "Sorting eigenvalues using " << rfToString(RF) << std::endl; + + const int Ls=16; + +// GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + //std::vector lat_size {16, 16, 16, 32}; + std::vector lat_size {32, 32, 32, 32}; + std::cout << "Lattice size: " << lat_size << std::endl; + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size, + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + // GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + // GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + GridCartesian * FGrid = UGrid; + GridRedBlackCartesian * FrbGrid = UrbGrid; + + std::vector seeds4({1,2,3,4}); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + LatticeFermion src(FGrid); random(RNG4, src); + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + NerscIO::readConfiguration(Umu, header, file); + + std::cout << GridLogMessage << "Loaded configuration" << std::endl; + + // RealD mass = 0.01; + RealD M5 = 1.8; + + // Wilson mass + RealD mass = -1.6; + + std::cout << GridLogMessage << "masses specified" << std::endl; + + std::vector boundary = {1,1,1,-1}; + WilsonFermionD::ImplParams Params(boundary); + + // DomainWallFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5); + // NonHermitianLinearOperator DLinOp (Ddwf); + + // WilsonFermionD Dwilson(Umu, *FGrid, *FrbGrid, mass); + WilsonFermionD Dwilson(Umu, *UGrid, *UrbGrid, mass, Params); + NonHermitianLinearOperator DLinOp (Dwilson); + + std::cout << GridLogMessage << "Dirac operator defined" << std::endl; + + // Define PV^dag D (if we want) + // DomainWallFermionD Dpv(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, 1.0, M5); + // typedef PVdagMLinearOperator PVdagM_t; + // PVdagM_t PVdagM(Ddwf, Dpv); + + std::cout< Dsq (DWilson); + // NonHermitianLinearOperator DLinOp (DWilson); + + int Nm = std::stoi(NmStr); + int Nk = std::stoi(NkStr); + int maxIter = std::stoi(maxIterStr); + int Nstop = std::stoi(NstopStr); + + std::cout << GridLogMessage << "Runnning Krylov Schur. Nm = " << Nm << ", Nk = " << Nk << ", maxIter = " << maxIter + << ", Nstop = " << Nstop << std::endl; + + // KrylovSchur KrySchur (PVdagM, FGrid, 1e-8, RF); // use PV^\dag M + KrylovSchur KrySchur (DLinOp, FGrid, 1e-8, RF); // use Ddwf + KrySchur(src, maxIter, Nm, Nk, Nstop); + + std::cout << GridLogMessage << "Checking eigensystem." << std::endl; + KrySchur.checkRitzEstimate(); + + std::cout<