portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-05-19 00:24:32 +01:00
Code Releases Activity

Starting Harmonic (shift and inverse)

Browse Source
This commit is contained in:
Chulwoo Jung
2025-11-24 17:05:35 -05:00
parent 0b457b9d52
commit 3538faf449
2 changed files with 22 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Grid/algorithms/iterative/KrylovSchur.h
+16 -2
View File
@@ -290,10 +290,11 @@ class KrylovSchur {
RitzFilter ritzFilter; // how to sort evals RitzFilter ritzFilter; // how to sort evals
public: public:
RealD *shift; // for Harmonic (shift and invert)
KrylovSchur(LinearOperatorBase<Field> &_Linop, GridBase *_Grid, RealD _Tolerance, RitzFilter filter = EvalReSmall) 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), : Linop(_Linop), Grid(_Grid), Tolerance(_Tolerance), ritzFilter(filter), u(_Grid), MaxIter(-1), Nm(-1), Nk(-1), Nstop (-1),
evals (0), ritzEstimates (), 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),shift(NULL)
{ {
u = Zero(); u = Zero();
}; };
@@ -316,6 +317,10 @@ class KrylovSchur {
* - Truncate the Krylov-Schur expansion. * - Truncate the Krylov-Schur expansion.
*/ */
void operator()(const Field& v0, int _maxIter, int _Nm, int _Nk, int _Nstop, bool doubleOrthog = true) { void operator()(const Field& v0, int _maxIter, int _Nm, int _Nk, int _Nstop, bool doubleOrthog = true) {
RealD shift_=1.;
shift = &shift_;
if (shift)
std::cout << GridLogMessage << "Shift " << approxLambdaMax << std::endl;
MaxIter = _maxIter; MaxIter = _maxIter;
Nm = _Nm; Nk = _Nk; Nm = _Nm; Nk = _Nk;
Nstop = _Nstop; Nstop = _Nstop;
@@ -347,7 +352,16 @@ class KrylovSchur {
// Perform a Schur decomposition on Rayleigh // Perform a Schur decomposition on Rayleigh
// ComplexSchurDecomposition schur (Rayleigh, false); // ComplexSchurDecomposition schur (Rayleigh, false);
Eigen::MatrixXcd temp = Rayleigh;
for (int m=0;m<Nm;m++) temp(m,m) -= *shift;
Eigen::MatrixXcd RayleighS = temp.inverse();
Eigen::MatrixXcd temp2 = RayleighS*temp;
std::cout << GridLogMessage << "Shift inverse check: shift= "<<*shift <<std::endl;
ComplexSchurDecomposition schur (Rayleigh, false, ritzFilter); ComplexSchurDecomposition schur (Rayleigh, false, ritzFilter);
ComplexSchurDecomposition schurS (RayleighS, false, ritzFilter);
std::cout << GridLogDebug << "Schur decomp holds? " << schur.checkDecomposition() << std::endl; std::cout << GridLogDebug << "Schur decomp holds? " << schur.checkDecomposition() << std::endl;
// Rearrange Schur matrix so wanted evals are on top left (like MATLAB's ordschur) // Rearrange Schur matrix so wanted evals are on top left (like MATLAB's ordschur)
@@ -457,7 +471,7 @@ class KrylovSchur {
basis.push_back(v0Norm * v0); // normalized source basis.push_back(v0Norm * v0); // normalized source
// basis[0] = v0Norm * v0; // normalized source // basis[0] = v0Norm * v0; // normalized source
Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm); Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm); // CJ: B in SLEPc
u = Zero(); u = Zero();
} else { } else {
// 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)
examples/Example_krylov_schur.cc
+5 -5
View File
@@ -368,8 +368,8 @@ int main (int argc, char ** argv)
LatticeGaugeField Umu(UGrid); LatticeGaugeField Umu(UGrid);
FieldMetaData header; FieldMetaData header;
// std::string file("ckpoint_lat.4000"); std::string file("config");
std::string file("/Users/patrickoare/libraries/PETSc-Grid/ckpoint_lat.4000"); // std::string file("Users/patrickoare/libraries/PETSc-Grid/ckpoint_lat.4000");
NerscIO::readConfiguration(Umu,header,file); NerscIO::readConfiguration(Umu,header,file);
RealD mass=0.01; RealD mass=0.01;
@@ -406,8 +406,8 @@ int main (int argc, char ** argv)
// int Nk = Nm+1; // if just running once // int Nk = Nm+1; // if just running once
// int maxIter = 5; // int maxIter = 5;
// int maxIter = 1; // int maxIter = 1;
int maxIter = 5; // int maxIter = 5;
// int maxIter = 100; int maxIter = 100;
int Nstop = 4; int Nstop = 4;
Coordinate origin ({0,0,0,0}); Coordinate origin ({0,0,0,0});
@@ -418,7 +418,7 @@ int main (int argc, char ** argv)
// Run KrylovSchur and Arnoldi on a Hermitian matrix // Run KrylovSchur and Arnoldi on a Hermitian matrix
std::cout << GridLogMessage << "Runnning Krylov Schur" << std::endl; std::cout << GridLogMessage << "Runnning Krylov Schur" << std::endl;
// KrylovSchur KrySchur (Dsq, FGrid, 1e-8, EvalNormLarge); // KrylovSchur KrySchur (Dsq, FGrid, 1e-8, EvalNormLarge);
KrylovSchur KrySchur (Dsq, FGrid, 1e-8); KrylovSchur KrySchur (Dsq, FGrid, 1e-8,EvalImNormSmall);
KrySchur(src, maxIter, Nm, Nk, Nstop); KrySchur(src, maxIter, Nm, Nk, Nstop);
/* /*