portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-03-26 22:16:08 +00:00
Code Releases Activity

added IO for evecs / evals

Browse Source
This commit is contained in:
Patrick Oare
2025-09-08 12:59:48 -04:00
parent 6fd71aea9d
commit 9dcd7ca761
2 changed files with 123 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -32,6 +32,12 @@ See the full license in the file "LICENSE" in the top level distribution directo
#ifndef GRID_KRYLOVSCHUR_H
#define GRID_KRYLOVSCHUR_H
// #include <Grid/Grid.h>
// #include <Grid/parallelIO/IldgIOtypes.h>
// #include <Grid/parallelIO/IldgIO.h>
// #include <Grid/serialisation/emptyUserRecord.h>
NAMESPACE_BEGIN(Grid);
/**
@@ -167,6 +173,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.
*/
@@ -212,6 +227,7 @@ class KrylovSchur {
/* Getters */
int getNk() { return Nk; }
Eigen::MatrixXcd getRayleighQuotient() { return Rayleigh; }
Field getU() { return u; }
std::vector<Field> getBasis() { return basis; }
@@ -321,6 +337,9 @@ class KrylovSchur {
<< i << "." << std::endl;
basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm);
std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl;
// writeEigensystem(path);
return;
}
}
@@ -413,8 +432,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
@@ -653,7 +670,49 @@ class KrylovSchur {
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

67
examples/Example_spec_kryschur.cc
View File

@@ -40,9 +40,57 @@ 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;
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();
for (int i = 0; i < Nk; i++) {
// write Eigenvalues
fEval << i << " " << evals(i);
if (i < Nk - 1) { fEval << "\n"; }
}
fEval.close();
// Write evecs
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
}
}
// Hermitize a DWF operator by squaring it
template<class Matrix,class Field>
class SquaredLinearOperator : public LinearOperatorBase<Field> {
@@ -293,12 +341,14 @@ int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
// Usage : $ ./Example_spec_kryschur <Nm> <Nk> <maxiter> <Nstop>
// Usage : $ ./Example_spec_kryschur <Nm> <Nk> <maxiter> <Nstop> <inFile> <outDir>
// assert (argc == 5);
std::string NmStr = argv[1];
std::string NkStr = argv[2];
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];
//const int Ls=16;
const int Ls = 8;
@@ -342,10 +392,11 @@ int main (int argc, char ** argv)
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");
// 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);
@@ -397,6 +448,10 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "Krylov Schur eigenvalues: " << std::endl << KrySchur.getEvals() << std::endl;
//std::cout << GridLogMessage << "Lanczos eigenvalues: " << std::endl << levals << std::endl;
// KrySchur.writeEigensystem(outDir);
writeEigensystem(KrySchur, outDir);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;