portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-03-23 12:36:09 +00:00
Code Releases Activity

Merge branch 'develop' of github.com:poare/Grid into develop

Browse Source
This commit is contained in:
Chulwoo Jung
2025-11-07 15:50:22 +00:00
parent caa66418bd 68af1bba67
commit fe0ab5f1a9
7 changed files with 1065 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Grid/algorithms/Algorithms.h
View File

@@ -81,7 +81,7 @@ NAMESPACE_CHECK(PowerMethod);
NAMESPACE_CHECK(multigrid);
#include <Grid/algorithms/FFT.h>
#include <Grid/algorithms/iterative/Arnoldi.h>
#include <Grid/algorithms/iterative/KrylovSchur.h>
#include <Grid/algorithms/iterative/Arnoldi.h>
#endif

9
Grid/algorithms/iterative/Arnoldi.h
View File

@@ -8,8 +8,6 @@ Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
Author: Chulwoo Jung <chulwoo@bnl.gov>
Author: Christoph Lehner <clehner@bnl.gov>
Author: Patrick Oare <poare@bnl.gov>
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.
*/

232
Grid/algorithms/iterative/KrylovSchur.h
View File

@@ -8,8 +8,6 @@ Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
Author: Chulwoo Jung <chulwoo@bnl.gov>
Author: Christoph Lehner <clehner@bnl.gov>
Author: Patrick Oare <poare@bnl.gov>
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<double> z1, std::complex<double> 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<class Field>
// 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<Field> 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<RealD> 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<Field> &_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<Field> getBasis() { return basis; }
Eigen::VectorXcd getEvals() { return evals; }
std::vector<Field> getEvecs() { return evecs; }
int getNk() { return Nk; }
Eigen::MatrixXcd getRayleighQuotient() { return Rayleigh; }
Field getU() { return u; }
std::vector<Field> getBasis() { return basis; }
Eigen::VectorXcd getEvals() { return evals; }
std::vector<RealD> getRitzEstimates() { return ritzEstimates; }
std::vector<Field> 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<Field> (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<Field> 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<Field> (basis.begin(), basis.begin() + Nk);
std::vector<Field> basisTmp = std::vector<Field> (basis.begin(), basis.begin() + Nk);
basis = basisTmp;
// evecs = std::vector<Field> (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= " <<start<< " basis.size= "<<basis.size()<<std::endl;
assert( start <= basis.size() ); // should be starting at the end of basis (start = Nk)
// assert( start <= basis.size() ); // should be starting at the end of basis (start = Nk)
// assert( start == basis.size() ); // should be starting at the end of basis (start = Nk)
std::cout << GridLogMessage << "Resetting Rayleigh and b" << std::endl;
Eigen::MatrixXcd RayleighCp = Rayleigh;
Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm);
@@ -372,12 +479,14 @@ class KrylovSchur {
// append b^\dag to Rayleigh, add u to basis
Rayleigh(Nk, Eigen::seqN(0, Nk)) = b.adjoint();
basis.push_back(u);
// basis[start] = u; // TODO make sure this is correct
b = Eigen::VectorXcd::Zero(Nm);
}
// Construct next Arnoldi vector by normalizing w_i = Dv_i - \sum_j v_j h_{ji}
for (int i = start; i < Nm; i++) {
Linop.Op(basis.back(), w);
// Linop.Op(basis[i], w);
for (int j = 0; j < basis.size(); j++) {
coeff = innerProduct(basis[j], w); // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after.
Rayleigh(j, i) = coeff;
@@ -400,6 +509,7 @@ class KrylovSchur {
basis.push_back(
(1.0/coeff) * w
);
// basis[i+1] = (1.0/coeff) * w;
}
// after iterations, update u and beta_k = ||u|| before norm
@@ -422,8 +532,6 @@ class KrylovSchur {
* the Rayleigh quotient. Assumes that the Rayleigh quotient has already been constructed (by
* calling the operator() function).
*
* TODO implement in parent class eventually.
*
* Parameters
* ----------
* Eigen::MatrixXcd& S
@@ -433,8 +541,11 @@ class KrylovSchur {
{
std::cout << GridLogMessage << "Computing eigenvalues." << std::endl;
evecs.clear();
evals = S.diagonal();
int n = evals.size(); // should be regular Nm
evecs.clear();
// evecs.assign(n, Field(Grid));
// TODO: is there a faster way to get the eigenvectors of a triangular matrix?
// Rayleigh.triangularView<Eigen::Upper> 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<Field> 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<Field> 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<Field>& UR, std::vector<Field> &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<Field>& RU, std::vector<Field> &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

2
examples/Example_krylov_schur.cc
View File

@@ -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;

150
examples/Example_spec_kryschur.cc
View File

@@ -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> <Nk> <maxiter> <Nstop> <inFile> <outDir> <?rf>
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 <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Patrick Oare <poare@bnl.edu>
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 <paboyle@ph.ed.ac.uk>
*************************************************************************************/
/* END LEGAL */
// copied here from Test_general_coarse_pvdagm.cc
// copied here from Test_general_coarse_pvdagm.cc
#include <cstdlib>
#include <Grid/Grid.h>
@@ -40,9 +59,13 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
#include <Grid/parallelIO/IldgIOtypes.h>
#include <Grid/parallelIO/IldgIO.h>
using namespace std;
using namespace Grid;
<<<<<<< HEAD
namespace Grid {
struct LanczosParameters: Serializable {
@@ -87,6 +110,54 @@ struct LanczosParameters: Serializable {
}
#if 0
=======
template <class T> 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 <class Field>
void writeEigensystem(KrylovSchur<Field> 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<RealD> 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<Field> 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 Matrix,class Field>
class SquaredLinearOperator : public LinearOperatorBase<Field> {
@@ -232,6 +303,7 @@ ShiftedComplexPVdagMLinearOperator(ComplexD _shift,Matrix &Mat,Matrix &PV): shif
}
};
<<<<<<< HEAD
template<class Fobj,class CComplex,int nbasis>
class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
@@ -334,16 +406,28 @@ public:
};
#endif
=======
>>>>>>> 68af1bba67dd62881ead5ab1e54962a5486a0791
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
// Usage : $ ./Example_spec_kryschur <Nm> <Nk> <maxiter> <Nstop>
// assert (argc == 5);
std::string NmStr = argv[1];
std::string NkStr = argv[2];
// Usage : $ ./Example_spec_kryschur <Nm> <Nk> <maaxiter> <Nstop> <inFile> <outDir>
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<clatt.size();d++){
clatt[d] = clatt[d]/2;
// clatt[d] = clatt[d]/4;
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> 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<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
NextToNearestStencilGeometry5D geom(Coarse5d);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
@@ -418,11 +474,6 @@ int main (int argc, char ** argv)
SquaredLinearOperator<DomainWallFermionD, LatticeFermionD> Dsq (Ddwf);
NonHermitianLinearOperator<DomainWallFermionD, LatticeFermionD> 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<<GridLogMessage << "*******************************************" << std::endl;
@@ -440,7 +491,8 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "*******************************************" << std::endl;
std::cout << GridLogMessage << "Krylov Schur eigenvalues: " << std::endl << KrySchur.getEvals() << std::endl;
//std::cout << GridLogMessage << "Lanczos eigenvalues: " << std::endl << levals << std::endl;
writeEigensystem(KrySchur, outDir);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;

383
examples/Example_wilson_evecs.cc Normal file
View File

@@ -0,0 +1,383 @@
/*************************************************************************************
Script for studying the Wilson eigenvectors resulting from the Krylov-Schur process.
Usage :
$ ./Example_spec_kryschur <Nm> <Nk> <maxiter> <Nstop> <inFile> <outDir> <?rf>
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 <paboyle@ph.ed.ac.uk>
Author: Patrick Oare <poare@bnl.edu>
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 <cstdlib>
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
#include <Grid/parallelIO/IldgIOtypes.h>
#include <Grid/parallelIO/IldgIO.h>
using namespace std;
using namespace Grid;
template <class T> 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 <class T> 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 <class Field>
void writeEigensystem(KrylovSchur<Field> 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<RealD> 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<Field> 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 Matrix,class Field>
class SquaredLinearOperator : public LinearOperatorBase<Field> {
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<Field> &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 "<<std::endl;
_Mat.M(in, out);
}
void _AdjOp (const Field &in, Field &out){
// std::cout << "AdjOp: Mdag "<<std::endl;
_Mat.Mdag(in, out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
// std::cout << "HermOp: Mdag M Mdag M"<<std::endl;
Field tmp(in.Grid());
_Op(in,tmp);
_AdjOp(tmp,out);
}
};
template<class Matrix,class Field>
class PVdagMLinearOperator : public LinearOperatorBase<Field> {
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<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
}
void AdjOp (const Field &in, Field &out){
std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(in,tmp);
_Mat.Mdag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
// _Mat.M(in,tmp);
// _PV.Mdag(tmp,out);
// _PV.M(out,tmp);
// _Mat.Mdag(tmp,out);
Op(in,tmp);
AdjOp(tmp,out);
// std::cout << "HermOp done "<<norm2(out)<<std::endl;
}
};
template<class Matrix,class Field>
class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
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<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
out = out + shift * in;
}
void AdjOp (const Field &in, Field &out){
std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(tmp,out);
_Mat.Mdag(in,tmp);
out = out + shift * in;
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
}
};
template<class Matrix, class Field>
class ShiftedComplexPVdagMLinearOperator : public LinearOperatorBase<Field> {
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<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
out = out + shift * in;
}
void AdjOp (const Field &in, Field &out){
std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(tmp,out);
_Mat.Mdag(in,tmp);
out = out + shift * in;
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
}
void resetShift(ComplexD newShift) {
shift = newShift;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
// Usage : $ ./Example_wilson_evecs ${inFile}
std::string file = argv[1];
const int Ls=16;
// GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
//std::vector<int> lat_size {16, 16, 16, 32};
std::vector<int> 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<int> 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<Complex> boundary = {1,1,1,-1};
WilsonFermionD::ImplParams Params(boundary);
// DomainWallFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5);
// NonHermitianLinearOperator<DomainWallFermionD, LatticeFermionD> DLinOp (Ddwf);
// WilsonFermionD Dwilson(Umu, *FGrid, *FrbGrid, mass);
WilsonFermionD Dwilson(Umu, *UGrid, *UrbGrid, mass, Params);
NonHermitianLinearOperator<WilsonFermionD, LatticeFermionD> 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<ComplexD> 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<LatticeFermion> 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<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}

374
examples/Example_wilson_spectrum.cc Normal file
View File

@@ -0,0 +1,374 @@
/*************************************************************************************
Runs the Krylov-Schur algorithm on a Wilson fermion operator to determine part of its spectrum.
TODO rename this file: really is running the topology change jobs on Aurora.
Usage :
$ ./Example_spec_kryschur <Nm> <Nk> <maxiter> <Nstop> <inFile> <outDir> <?rf>
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 <paboyle@ph.ed.ac.uk>
Author: Patrick Oare <poare@bnl.edu>
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 <cstdlib>
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
#include <Grid/parallelIO/IldgIOtypes.h>
#include <Grid/parallelIO/IldgIO.h>
using namespace std;
using namespace Grid;
template <class T> 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 <class Field>
void writeEigensystem(KrylovSchur<Field> 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<RealD> 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<Field> 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 Matrix,class Field>
class SquaredLinearOperator : public LinearOperatorBase<Field> {
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<Field> &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 "<<std::endl;
_Mat.M(in, out);
}
void _AdjOp (const Field &in, Field &out){
// std::cout << "AdjOp: Mdag "<<std::endl;
_Mat.Mdag(in, out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
// std::cout << "HermOp: Mdag M Mdag M"<<std::endl;
Field tmp(in.Grid());
_Op(in,tmp);
_AdjOp(tmp,out);
}
};
template<class Matrix,class Field>
class PVdagMLinearOperator : public LinearOperatorBase<Field> {
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<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
}
void AdjOp (const Field &in, Field &out){
std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(in,tmp);
_Mat.Mdag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
// _Mat.M(in,tmp);
// _PV.Mdag(tmp,out);
// _PV.M(out,tmp);
// _Mat.Mdag(tmp,out);
Op(in,tmp);
AdjOp(tmp,out);
// std::cout << "HermOp done "<<norm2(out)<<std::endl;
}
};
template<class Matrix,class Field>
class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
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<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
out = out + shift * in;
}
void AdjOp (const Field &in, Field &out){
std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(tmp,out);
_Mat.Mdag(in,tmp);
out = out + shift * in;
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
}
};
template<class Matrix, class Field>
class ShiftedComplexPVdagMLinearOperator : public LinearOperatorBase<Field> {
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<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
out = out + shift * in;
}
void AdjOp (const Field &in, Field &out){
std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(tmp,out);
_Mat.Mdag(in,tmp);
out = out + shift * in;
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
}
void resetShift(ComplexD newShift) {
shift = newShift;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
// Usage : $ ./Example_spec_kryschur <Nm> <Nk> <maaxiter> <Nstop> <inFile> <outDir>
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<int> lat_size {16, 16, 16, 32};
std::vector<int> 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<int> 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<Complex> boundary = {1,1,1,-1};
WilsonFermionD::ImplParams Params(boundary);
// DomainWallFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5);
// NonHermitianLinearOperator<DomainWallFermionD, LatticeFermionD> DLinOp (Ddwf);
// WilsonFermionD Dwilson(Umu, *FGrid, *FrbGrid, mass);
WilsonFermionD Dwilson(Umu, *UGrid, *UrbGrid, mass, Params);
NonHermitianLinearOperator<WilsonFermionD, LatticeFermionD> 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<DomainWallFermionD,LatticeFermionD> PVdagM_t;
// PVdagM_t PVdagM(Ddwf, Dpv);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
// SquaredLinearOperator<WilsonFermionD, LatticeFermionD> Dsq (DWilson);
// NonHermitianLinearOperator<WilsonFermionD, LatticeFermionD> 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<<GridLogMessage << "*******************************************" << std::endl;
std::cout<<GridLogMessage << "***************** RESULTS *****************" << std::endl;
std::cout<<GridLogMessage << "*******************************************" << std::endl;
std::cout << GridLogMessage << "Krylov Schur eigenvalues: " << std::endl << KrySchur.getEvals() << std::endl;
writeEigensystem(KrySchur, outDir);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}