portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-04 14:04:32 +00:00
Code Releases Activity

first complete version of IRBL; requires practical test and clean up

Browse Source
This commit is contained in:
Yong-Chull Jang
2017-12-21 23:13:39 -05:00
parent fe406e230d
commit 89c4e9b168
3 changed files with 24 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
View File

@@ -198,18 +198,7 @@ until convergence
// clog << "ckpt A2: lmd[" << k << "] = " << lmd[0][k] << '\n'; // clog << "ckpt A2: lmd[" << k << "] = " << lmd[0][k] << '\n';
//} //}
// residual vector
#if 1 // ypj[fixme] temporary to check a case when block has one vector
for ( int i=0; i<Nu; ++i) f_copy[i] = f[i];
for ( int i=0; i<Nu; ++i) {
f[i] = f_copy[0]*lme[0][Nm-Nu+i];
for ( int j=1; j<Nu; ++j) {
f[i] += f_copy[j]*lme[j][Nm-Nu+i];
}
//clog << "ckpt C (i= " << i << ")" << '\n';
//clog << "norm2(f) = " << norm2(f[i]) << std::endl;
}
#endif
// getting eigenvalues // getting eigenvalues
for(int u=0; u<Nu; ++u){ for(int u=0; u<Nu; ++u){
@@ -293,7 +282,19 @@ until convergence
//clog << "ckpt F: norm2_evec[= " << i << "]" << norm2(evec[i]) << std::endl; //clog << "ckpt F: norm2_evec[= " << i << "]" << norm2(evec[i]) << std::endl;
} }
#if 1 // ypj[fixme] temporary to check a case when block has one vector // residual vector
#if 0 // ypj[fixme] temporary to check a case when block has one vector
for ( int i=0; i<Nu; ++i) f_copy[i] = f[i];
for ( int i=0; i<Nu; ++i) {
f[i] = f_copy[0]*lme[0][Nm-Nu+i];
for ( int j=1; j<Nu; ++j) {
f[i] += f_copy[j]*lme[j][Nm-Nu+i];
}
//clog << "ckpt C (i= " << i << ")" << '\n';
//clog << "norm2(f) = " << norm2(f[i]) << std::endl;
}
// ypj[fixme] temporary to check a case when block has one vector
// Compressed vector f and beta(k2) // Compressed vector f and beta(k2)
f[0] *= Q(Nm-1,Nk-1); f[0] *= Q(Nm-1,Nk-1);
f[0] += lme[0][Nk-1] * evec[Nk]; // was commented out f[0] += lme[0][Nk-1] * evec[Nk]; // was commented out
@@ -305,6 +306,8 @@ until convergence
RealD betar = 1.0/beta_k; RealD betar = 1.0/beta_k;
evec[Nk] = betar * f[0]; evec[Nk] = betar * f[0];
lme[0][Nk-1] = beta_k; lme[0][Nk-1] = beta_k;
#else
blockwiseStep(lmd,lme,evec,f,f_copy,Nblock_k-1);
#endif #endif
// Convergence test // Convergence test
@@ -393,7 +396,7 @@ until convergence
clog << fname + " CONVERGED ; Summary :\n"; clog << fname + " CONVERGED ; Summary :\n";
clog <<"**************************************************************************"<< std::endl; clog <<"**************************************************************************"<< std::endl;
clog << " -- Iterations = "<< iter << "\n"; clog << " -- Iterations = "<< iter << "\n";
clog << " -- beta(k) = "<< beta_k << "\n"; //clog << " -- beta(k) = "<< beta_k << "\n";
clog << " -- Nconv = "<< Nconv << "\n"; clog << " -- Nconv = "<< Nconv << "\n";
clog <<"**************************************************************************"<< std::endl; clog <<"**************************************************************************"<< std::endl;
} }

6
tests/lanczos/Test_dwf_block_lanczos.cc
View File

@@ -77,10 +77,10 @@ int main (int argc, char ** argv)
SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf); SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
// SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf); // SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
const int Nstop = 50; const int Nstop = 120;
const int Nu = 1; const int Nu = 1;
const int Nk = 200; const int Nk = 240;
const int Np = 200; const int Np = 240;
const int Nm = Nk+Np; const int Nm = Nk+Np;
const int MaxIt= 10; const int MaxIt= 10;
RealD resid = 1.0e-8; RealD resid = 1.0e-8;

8
tests/lanczos/Test_dwf_lanczos.cc
View File

@@ -75,11 +75,11 @@ int main (int argc, char ** argv)
SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf); SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
// SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf); // SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
const int Nstop = 50; const int Nstop = 120;
const int Nk = 200; const int Nk = 240;
const int Np = 200; const int Np = 240;
const int Nm = Nk+Np; const int Nm = Nk+Np;
const int MaxIt= 100; const int MaxIt= 10;
RealD resid = 1.0e-8; RealD resid = 1.0e-8;
std::vector<double> Coeffs { 0.,-1.}; std::vector<double> Coeffs { 0.,-1.};