diff --git a/tests/lanczos/LanParams.xml b/tests/lanczos/LanParams.xml
index 635fb243..5328b46b 100644
--- a/tests/lanczos/LanParams.xml
+++ b/tests/lanczos/LanParams.xml
@@ -1,6 +1,14 @@
- -3.5
+ 0.00107
+ 1.8
+ 48
+ 10
+ 15
+ 85
+ 0.003
+ 60
+ 201
diff --git a/tests/lanczos/Test_dwf_G5R5.cc b/tests/lanczos/Test_dwf_G5R5.cc
new file mode 100644
index 00000000..0fd1bc35
--- /dev/null
+++ b/tests/lanczos/Test_dwf_G5R5.cc
@@ -0,0 +1,346 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_dwf_G5R5.cc
+
+Copyright (C) 2015
+
+Author: Chulwoo Jung
+From Duo and Bob's Chirality study
+
+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
+
+using namespace std;
+using namespace Grid;
+ ;
+
+//typedef WilsonFermionD FermionOp;
+typedef DomainWallFermionD FermionOp;
+typedef typename DomainWallFermionD::FermionField FermionField;
+
+
+RealD AllZero(RealD x) { return 0.; }
+
+namespace Grid {
+
+struct LanczosParameters: Serializable {
+ GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters,
+ RealD, mass ,
+ RealD, M5 ,
+ Integer, Ls,
+ Integer, Nstop,
+ Integer, Nk,
+ Integer, Np,
+ RealD, ChebyLow,
+ RealD, ChebyHigh,
+ Integer, ChebyOrder)
+// Integer, StartTrajectory,
+// Integer, Trajectories, /* @brief Number of sweeps in this run */
+// bool, MetropolisTest,
+// Integer, NoMetropolisUntil,
+// std::string, StartingType,
+// Integer, SW,
+// RealD, Kappa,
+// IntegratorParameters, MD)
+
+ LanczosParameters() {
+ ////////////////////////////// Default values
+ mass = 0;
+// MetropolisTest = true;
+// NoMetropolisUntil = 10;
+// StartTrajectory = 0;
+// SW = 2;
+// Trajectories = 10;
+// StartingType = "HotStart";
+ /////////////////////////////////
+ }
+
+ template
+ LanczosParameters(Reader & TheReader){
+ initialize(TheReader);
+ }
+
+ template < class ReaderClass >
+ void initialize(Reader &TheReader){
+// std::cout << GridLogMessage << "Reading HMC\n";
+ read(TheReader, "HMC", *this);
+ }
+
+
+ void print_parameters() const {
+// std::cout << GridLogMessage << "[HMC parameters] Trajectories : " << Trajectories << "\n";
+// std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n";
+// std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
+// std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n";
+// std::cout << GridLogMessage << "[HMC parameters] Starting type : " << StartingType << "\n";
+// MD.print_parameters();
+ }
+
+};
+
+}
+
+int main(int argc, char** argv) {
+ Grid_init(&argc, &argv);
+
+ LanczosParameters LanParams;
+#if 1
+ {
+ XmlReader HMCrd("LanParams.xml");
+ read(HMCrd,"LanczosParameters",LanParams);
+ }
+#else
+ {
+ LanParams.mass = mass;
+ }
+#endif
+ std::cout << GridLogMessage<< LanParams < seeds4({1, 2, 3, 4});
+ std::vector seeds5({5, 6, 7, 8});
+ GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
+ GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
+ GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5);
+
+ LatticeGaugeField Umu(UGrid);
+
+ FieldMetaData header;
+ std::string file("./config");
+
+ int precision32 = 0;
+ int tworow = 0;
+ NerscIO::readConfiguration(Umu,header,file);
+
+/*
+ std::vector U(4, UGrid);
+ for (int mu = 0; mu < Nd; mu++) {
+ U[mu] = PeekIndex(Umu, mu);
+ }
+*/
+
+ int Nstop = 10;
+ int Nk = 20;
+ int Np = 80;
+ Nstop=LanParams.Nstop;
+ Nk=LanParams.Nk;
+ Np=LanParams.Np;
+
+ int Nm = Nk + Np;
+ int MaxIt = 10000;
+ RealD resid = 1.0e-5;
+
+
+//while ( mass > - 5.0){
+ FermionOp Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+ MdagMLinearOperator HermOp(Ddwf); /// <-----
+// Gamma5HermitianLinearOperator HermOp2(WilsonOperator); /// <-----
+ Gamma5R5HermitianLinearOperator G5R5Herm(Ddwf);
+// Gamma5R5HermitianLinearOperator
+ std::vector Coeffs{0, 1.};
+ Polynomial PolyX(Coeffs);
+
+ Chebyshev Cheby(LanParams.ChebyLow,LanParams.ChebyHigh,LanParams.ChebyOrder);
+
+ FunctionHermOp OpCheby(Cheby,HermOp);
+ PlainHermOp Op (HermOp);
+ PlainHermOp Op2 (G5R5Herm);
+
+ ImplicitlyRestartedLanczos IRL(OpCheby, Op, Nstop, Nk, Nm, resid, MaxIt);
+
+ std::vector eval(Nm);
+ FermionField src(FGrid);
+ gaussian(RNG5, src);
+ std::vector evec(Nm, FGrid);
+ for (int i = 0; i < 1; i++) {
+ std::cout << i << " / " << Nm << " grid pointer " << evec[i].Grid()
+ << std::endl;
+ };
+
+ int Nconv;
+ IRL.calc(eval, evec, src, Nconv);
+
+ std::cout << mass <<" : " << eval << std::endl;
+
+#if 0
+ Gamma g5(Gamma::Algebra::Gamma5) ;
+ ComplexD dot;
+ FermionField tmp(FGrid);
+// RealD eMe,eMMe;
+ for (int i = 0; i < Nstop ; i++) {
+// tmp = g5*evec[i];
+ dot = innerProduct(evec[i],evec[i]);
+// G5R5(tmp,evec[i]);
+ G5R5Herm.HermOpAndNorm(evec[i],tmp,eMe,eMMe);
+ std::cout <<"Norm "< G5R5Mevec(Nconv, FGrid);
+ vector finalevec(Nconv, FGrid);
+ vector eMe(Nconv), eMMe(Nconv);
+ for(int i = 0; i < Nconv; i++){
+ G5R5Herm.HermOpAndNorm(evec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
+ }
+ cout << "Re: " << endl;
+ cout << eMe << endl;
+ cout << "" << endl;
+ cout << eMMe << endl;
+ vector> VevecG5R5Mevec(Nconv);
+ Eigen::MatrixXcd evecG5R5Mevec = Eigen::MatrixXcd::Zero(Nconv, Nconv);
+ for(int i = 0; i < Nconv; i++){
+ VevecG5R5Mevec[i].resize(Nconv);
+ for(int j = 0; j < Nconv; j++){
+ VevecG5R5Mevec[i][j] = innerProduct(evec[i], G5R5Mevec[j]);
+ evecG5R5Mevec(i, j) = VevecG5R5Mevec[i][j];
+ }
+ }
+ //calculate eigenvector
+ cout << "Eigen solver" << endl;
+ Eigen::SelfAdjointEigenSolver eigensolver(evecG5R5Mevec);
+ vector eigeneval(Nconv);
+ vector> eigenevec(Nconv);
+ for(int i = 0; i < Nconv; i++){
+ eigeneval[i] = eigensolver.eigenvalues()[i];
+ eigenevec[i].resize(Nconv);
+ for(int j = 0; j < Nconv; j++){
+ eigenevec[i][j] = eigensolver.eigenvectors()(i, j);
+ }
+ }
+ //rotation
+ cout << "Do rotation" << endl;
+ for(int i = 0; i < Nconv; i++){
+ finalevec[i] = finalevec[i] - finalevec[i];
+ for(int j = 0; j < Nconv; j++){
+ finalevec[i] = eigenevec[j][i]*evec[j] + finalevec[i];
+ }
+ }
+ //normalize again;
+ for(int i = 0; i < Nconv; i++){
+ RealD tmp_RealD = norm2(finalevec[i]);
+ tmp_RealD = 1./pow(tmp_RealD, 0.5);
+ finalevec[i] = finalevec[i]*tmp_RealD;
+ }
+
+ //check
+ for(int i = 0; i < Nconv; i++){
+ G5R5Herm.HermOpAndNorm(finalevec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
+ }
+
+ //**********************************************************************
+ //sort the eigenvectors
+ vector finalevec_copy(Nconv, FGrid);
+ for(int i = 0; i < Nconv; i++){
+ finalevec_copy[i] = finalevec[i];
+ }
+ vector eMe_copy(eMe);
+ for(int i = 0; i < Nconv; i++){
+ eMe[i] = fabs(eMe[i]);
+ eMe_copy[i] = eMe[i];
+ }
+ sort(eMe_copy.begin(), eMe_copy.end());
+ for(int i = 0; i < Nconv; i++){
+ for(int j = 0; j < Nconv; j++){
+ if(eMe[j] == eMe_copy[i]){
+ finalevec[i] = finalevec_copy[j];
+ }
+ }
+ }
+ for(int i = 0; i < Nconv; i++){
+ G5R5Herm.HermOpAndNorm(finalevec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
+ }
+ cout << "Re: " << endl;
+ cout << eMe << endl;
+ cout << "" << endl;
+ cout << eMMe << endl;
+
+
+// vector finalevec(Nconv, FGrid);
+// temporary, until doing rotation
+// for(int i = 0; i < Nconv; i++)
+// finalevec[i]=evec[i];
+ //**********************************************************************
+ //calculate chirality matrix
+ vector G5evec(Nconv, FGrid);
+ vector> chiral_matrix(Nconv);
+ vector> chiral_matrix_real(Nconv);
+ for(int i = 0; i < Nconv; i++){
+// G5evec[i] = G5evec[i] - G5evec[i];
+ G5evec[i] = Zero();
+ for(int j = 0; j < Ls/2; j++){
+ axpby_ssp(G5evec[i], 1., finalevec[i], 0., G5evec[i], j, j);
+ }
+ for(int j = Ls/2; j < Ls; j++){
+ axpby_ssp(G5evec[i], -1., finalevec[i], 0., G5evec[i], j, j);
+ }
+ }
+ for(int i = 0; i < Nconv; i++){
+ chiral_matrix_real[i].resize(Nconv);
+ chiral_matrix[i].resize(Nconv);
+ for(int j = 0; j < Nconv; j++){
+ chiral_matrix[i][j] = innerProduct(finalevec[i], G5evec[j]);
+ chiral_matrix_real[i][j] = abs(chiral_matrix[i][j]);
+ std::cout <<" chiral_matrix_real "<