portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-15 02:05:37 +00:00
Code Releases Activity

add explicit restart method rbl

Browse Source
This commit is contained in:
Yong-Chull Jang 2018-03-05 18:12:42 -05:00
parent 3caf0e8b09
commit 11219a8f7a
2 changed files with 309 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -33,10 +33,12 @@ Author: Guido Cossu
#include <string.h> //memset #include <string.h> //memset
#define clog std::cout << GridLogMessage #define Glog std::cout << GridLogMessage
namespace Grid { namespace Grid {
enum class LanczosType { irbl, rbl };
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
// Implicitly restarted block lanczos // Implicitly restarted block lanczos
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
@ -53,6 +55,7 @@ private:
int Nm; // total number of vectors int Nm; // total number of vectors
int Nblock_k; // Nk/Nu int Nblock_k; // Nk/Nu
int Nblock_m; // Nm/Nu int Nblock_m; // Nm/Nu
int Nconv_test_interval; // Number of skipped vectors when checking a convergence
RealD eresid; RealD eresid;
IRLdiagonalisation diagonalisation; IRLdiagonalisation diagonalisation;
//////////////////////////////////// ////////////////////////////////////
@ -69,6 +72,7 @@ public:
ImplicitlyRestartedBlockLanczos(LinearOperatorBase<Field> &Linop, // op ImplicitlyRestartedBlockLanczos(LinearOperatorBase<Field> &Linop, // op
OperatorFunction<Field> & poly, // polynomial OperatorFunction<Field> & poly, // polynomial
int _Nstop, // really sought vecs int _Nstop, // really sought vecs
int _Nconv_test_interval, // conv check interval
int _Nu, // vecs in the unit block int _Nu, // vecs in the unit block
int _Nk, // sought vecs int _Nk, // sought vecs
int _Nm, // total vecs int _Nm, // total vecs
@ -76,7 +80,8 @@ public:
int _MaxIter, // Max iterations int _MaxIter, // Max iterations
IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen) IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen)
: _Linop(Linop), _poly(poly), : _Linop(Linop), _poly(poly),
Nstop(_Nstop), Nu(_Nu), Nk(_Nk), Nm(_Nm), Nstop(_Nstop), Nconv_test_interval(_Nconv_test_interval),
Nu(_Nu), Nk(_Nk), Nm(_Nm),
Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu), Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu),
//eresid(_eresid), MaxIter(10), //eresid(_eresid), MaxIter(10),
eresid(_eresid), MaxIter(_MaxIter), eresid(_eresid), MaxIter(_MaxIter),
@ -117,30 +122,45 @@ public:
} }
normalize(w); normalize(w);
} }
void calc(std::vector<RealD>& eval, void calc(std::vector<RealD>& eval,
std::vector<Field>& evec, std::vector<Field>& evec,
const std::vector<Field>& src, int& Nconv) const std::vector<Field>& src, int& Nconv, LanczosType Impl)
{ {
std::string fname = std::string(cname+"::calc()"); switch (Impl) {
case LanczosType::irbl:
calc_irbl(eval,evec,src,Nconv);
break;
case LanczosType::rbl:
calc_rbl(eval,evec,src,Nconv);
break;
}
}
void calc_irbl(std::vector<RealD>& eval,
std::vector<Field>& evec,
const std::vector<Field>& src, int& Nconv)
{
std::string fname = std::string(cname+"::calc_irbl()");
GridBase *grid = evec[0]._grid; GridBase *grid = evec[0]._grid;
assert(grid == src[0]._grid); assert(grid == src[0]._grid);
assert( Nu = src.size() ); assert( Nu = src.size() );
clog << std::string(74,'*') << std::endl; Glog << std::string(74,'*') << std::endl;
clog << fname + " starting iteration 0 / "<< MaxIter<< std::endl; Glog << fname + " starting iteration 0 / "<< MaxIter<< std::endl;
clog << std::string(74,'*') << std::endl; Glog << std::string(74,'*') << std::endl;
clog <<" -- seek Nk = "<< Nk <<" vectors"<< std::endl; Glog <<" -- seek Nk = "<< Nk <<" vectors"<< std::endl;
clog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl; Glog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl;
clog <<" -- total Nm = "<< Nm <<" vectors"<< std::endl; Glog <<" -- total Nm = "<< Nm <<" vectors"<< std::endl;
clog <<" -- size of eval = "<< eval.size() << std::endl; Glog <<" -- size of eval = "<< eval.size() << std::endl;
clog <<" -- size of evec = "<< evec.size() << std::endl; Glog <<" -- size of evec = "<< evec.size() << std::endl;
if ( diagonalisation == IRLdiagonaliseWithEigen ) { if ( diagonalisation == IRLdiagonaliseWithEigen ) {
clog << "Diagonalisation is Eigen "<< std::endl; Glog << "Diagonalisation is Eigen "<< std::endl;
} else { } else {
abort(); abort();
} }
clog << std::string(74,'*') << std::endl; Glog << std::string(74,'*') << std::endl;
assert(Nm == evec.size() && Nm == eval.size()); assert(Nm == evec.size() && Nm == eval.size());
@ -167,10 +187,10 @@ public:
// set initial vector // set initial vector
for (int i=0; i<Nu; ++i) { for (int i=0; i<Nu; ++i) {
clog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl; Glog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
evec[i] = src[i]; evec[i] = src[i];
orthogonalize(evec[i],evec,i); orthogonalize(evec[i],evec,i);
clog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl; Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
} }
// initial Nblock_k steps // initial Nblock_k steps
@ -180,7 +200,7 @@ public:
int iter; int iter;
for(iter = 0; iter<MaxIter; ++iter){ for(iter = 0; iter<MaxIter; ++iter){
clog <<"#Restart iteration = "<< iter << std::endl; Glog <<"#Restart iteration = "<< iter << std::endl;
// additional (Nblock_m - Nblock_k) steps // additional (Nblock_m - Nblock_k) steps
for(int b=Nblock_k; b<Nblock_m; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b); for(int b=Nblock_k; b<Nblock_m; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
@ -194,7 +214,7 @@ public:
Qt = Eigen::MatrixXcd::Identity(Nm,Nm); Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid); diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid);
_sort.push(eval2,Nm); _sort.push(eval2,Nm);
clog << "#Ritz value before shift: "<< std::endl; Glog << "#Ritz value before shift: "<< std::endl;
for(int i=0; i<Nm; ++i){ for(int i=0; i<Nm; ++i){
std::cout.precision(13); std::cout.precision(13);
std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] "; std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
@ -203,7 +223,7 @@ public:
//---------------------------------------------------------------------- //----------------------------------------------------------------------
if ( Nm>Nk ) { if ( Nm>Nk ) {
clog <<" #Apply shifted QR transformations "<<std::endl; Glog <<" #Apply shifted QR transformations "<<std::endl;
//int k2 = Nk+Nu; //int k2 = Nk+Nu;
int k2 = Nk; int k2 = Nk;
@ -240,7 +260,7 @@ public:
Qt = Eigen::MatrixXcd::Identity(Nm,Nm); Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid); diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
_sort.push(eval2,Nk); _sort.push(eval2,Nk);
clog << "#Ritz value after shift: "<< std::endl; Glog << "#Ritz value after shift: "<< std::endl;
for(int i=0; i<Nk; ++i){ for(int i=0; i<Nk; ++i){
std::cout.precision(13); std::cout.precision(13);
std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] "; std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
@ -250,7 +270,7 @@ public:
//---------------------------------------------------------------------- //----------------------------------------------------------------------
// Convergence test // Convergence test
clog <<" #Convergence test: "<<std::endl; Glog <<" #Convergence test: "<<std::endl;
for(int k = 0; k<Nk; ++k) B[k]=0.0; for(int k = 0; k<Nk; ++k) B[k]=0.0;
for(int j = 0; j<Nk; ++j){ for(int j = 0; j<Nk; ++j){
for(int k = 0; k<Nk; ++k){ for(int k = 0; k<Nk; ++k){
@ -284,7 +304,7 @@ public:
} // i-loop end } // i-loop end
clog <<" #modes converged: "<<Nconv<<std::endl; Glog <<" #modes converged: "<<Nconv<<std::endl;
for(int i=0; i<Nconv; ++i){ for(int i=0; i<Nconv; ++i){
std::cout.precision(13); std::cout.precision(13);
std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<Iconv[i]<<"] "; std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<Iconv[i]<<"] ";
@ -296,11 +316,11 @@ public:
} // end of iter loop } // end of iter loop
clog << std::string(74,'*') << std::endl; Glog << std::string(74,'*') << std::endl;
if ( Nconv<Nstop ) { if ( Nconv<Nstop ) {
clog << fname + " NOT converged ; Summary :\n"; Glog << fname + " NOT converged ; Summary :\n";
} else { } else {
clog << fname + " CONVERGED ; Summary :\n"; Glog << fname + " CONVERGED ; Summary :\n";
// Sort convered eigenpairs. // Sort convered eigenpairs.
eval.resize(Nconv); eval.resize(Nconv);
evec.resize(Nconv,grid); evec.resize(Nconv,grid);
@ -310,11 +330,185 @@ public:
} }
_sort.push(eval,evec,Nconv); _sort.push(eval,evec,Nconv);
} }
clog << std::string(74,'*') << std::endl; Glog << std::string(74,'*') << std::endl;
clog << " -- Iterations = "<< iter << "\n"; Glog << " -- Iterations = "<< iter << "\n";
//clog << " -- beta(k) = "<< beta_k << "\n"; //Glog << " -- beta(k) = "<< beta_k << "\n";
clog << " -- Nconv = "<< Nconv << "\n"; Glog << " -- Nconv = "<< Nconv << "\n";
clog << std::string(74,'*') << std::endl; Glog << std::string(74,'*') << std::endl;
}
void calc_rbl(std::vector<RealD>& eval,
std::vector<Field>& evec,
const std::vector<Field>& src, int& Nconv)
{
std::string fname = std::string(cname+"::calc_rbl()");
GridBase *grid = evec[0]._grid;
assert(grid == src[0]._grid);
assert( Nu = src.size() );
int Np = (Nm-Nk);
if (Np > 0 && MaxIter > 1) Np /= MaxIter;
int Nblock_p = Np/Nu;
Glog << std::string(74,'*') << std::endl;
Glog << fname + " starting iteration 0 / "<< MaxIter<< std::endl;
Glog << std::string(74,'*') << std::endl;
Glog <<" -- seek (min) Nk = "<< Nk <<" vectors"<< std::endl;
Glog <<" -- seek (inc) Np = "<< Np <<" vectors"<< std::endl;
Glog <<" -- seek (max) Nm = "<< Nm <<" vectors"<< std::endl;
Glog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl;
Glog <<" -- size of eval = "<< eval.size() << std::endl;
Glog <<" -- size of evec = "<< evec.size() << std::endl;
if ( diagonalisation == IRLdiagonaliseWithEigen ) {
Glog << "Diagonalisation is Eigen "<< std::endl;
} else {
abort();
}
Glog << std::string(74,'*') << std::endl;
assert(Nm == evec.size() && Nm == eval.size());
std::vector<std::vector<ComplexD>> lmd(Nu,std::vector<ComplexD>(Nm,0.0));
std::vector<std::vector<ComplexD>> lme(Nu,std::vector<ComplexD>(Nm,0.0));
std::vector<std::vector<ComplexD>> lmd2(Nu,std::vector<ComplexD>(Nm,0.0));
std::vector<std::vector<ComplexD>> lme2(Nu,std::vector<ComplexD>(Nm,0.0));
std::vector<RealD> eval2(Nk);
std::vector<RealD> resid(Nm);
Eigen::MatrixXcd Qt = Eigen::MatrixXcd::Zero(Nm,Nm);
Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
std::vector<int> Iconv(Nm);
std::vector<Field> B(Nm,grid); // waste of space replicating
std::vector<Field> f(Nu,grid);
std::vector<Field> f_copy(Nu,grid);
Field v(grid);
Nconv = 0;
RealD beta_k;
// set initial vector
for (int i=0; i<Nu; ++i) {
Glog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
evec[i] = src[i];
orthogonalize(evec[i],evec,i);
Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
}
// initial Nblock_k steps
for(int b=0; b<Nblock_k; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
// restarting loop begins
int iter;
int Nblock_l, Nblock_r;
int Nl, Nr;
int Nconv_guess = 0;
for(iter = 0; iter<MaxIter; ++iter){
Glog <<"#Restart iteration = "<< iter << std::endl;
Nblock_l = Nblock_k + iter*Nblock_p;
Nblock_r = Nblock_l + Nblock_p;
Nl = Nblock_l*Nu;
Nr = Nblock_r*Nu;
eval2.resize(Nr);
// additional Nblock_p steps
for(int b=Nblock_l; b<Nblock_r; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
// getting eigenvalues
for(int u=0; u<Nu; ++u){
for(int k=0; k<Nr; ++k){
lmd2[u][k] = lmd[u][k];
lme2[u][k] = lme[u][k];
}
}
Qt = Eigen::MatrixXcd::Identity(Nr,Nr);
diagonalize(eval2,lmd2,lme2,Nu,Nr,Nr,Qt,grid);
_sort.push(eval2,Nr);
Glog << "#Ritz value: "<< std::endl;
for(int i=0; i<Nr; ++i){
std::cout.precision(13);
std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
}
// Convergence test
Glog <<" #Convergence test: "<<std::endl;
for(int k = 0; k<Nr; ++k) B[k]=0.0;
for(int j = 0; j<Nr; j+=Nconv_test_interval){
Glog <<" #rotation for evec"
<< std::setw(4)<< std::setiosflags(std::ios_base::right)
<< "["<< j <<"]" <<std::endl;
for(int k = 0; k<Nr; ++k){
B[j].checkerboard = evec[k].checkerboard;
B[j] += evec[k]*Qt(k,j);
}
}
Nconv = 0;
for(int i=0; i<Nr; i+=Nconv_test_interval){
_Linop.HermOp(B[i],v);
RealD vnum = real(innerProduct(B[i],v)); // HermOp.
RealD vden = norm2(B[i]);
eval2[i] = vnum/vden;
v -= eval2[i]*B[i];
RealD vv = norm2(v);
resid[i] = vv;
std::cout.precision(13);
std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
std::cout << "eval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i];
std::cout << " resid^2 = "<< std::setw(20)<< std::setiosflags(std::ios_base::right)<< vv<< std::endl;
// change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged
//if( (vv<eresid*eresid) && (i == Nconv) ){
if (vv<eresid*eresid) {
Iconv[Nconv] = i;
++Nconv;
}
} // i-loop end
Glog <<" #modes converged: "<<Nconv<<std::endl;
for(int i=0; i<Nconv; ++i){
std::cout.precision(13);
std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<Iconv[i]<<"] ";
std::cout << "eval_conv = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[Iconv[i]];
std::cout << " resid^2 = "<< std::setw(20)<< std::setiosflags(std::ios_base::right)<< resid[Iconv[i]]<< std::endl;
}
(Nconv > 0 ) ? Nconv_guess = 1 + (Nconv-1)*Nconv_test_interval : Nconv_guess = 0;
if ( Nconv_guess >= Nstop ) break;
} // end of iter loop
Glog << std::string(74,'*') << std::endl;
if ( Nconv_guess < Nstop ) {
Glog << fname + " NOT converged ; Summary :\n";
} else {
Glog << fname + " CONVERGED ; Summary :\n";
// Sort convered eigenpairs.
eval.resize(Nconv);
evec.resize(Nconv,grid);
for(int i=0; i<Nconv; ++i){
eval[i] = eval2[Iconv[i]];
evec[i] = B[Iconv[i]];
}
_sort.push(eval,evec,Nconv);
}
Glog << std::string(74,'*') << std::endl;
Glog << " -- Iterations = "<< iter << "\n";
//Glog << " -- beta(k) = "<< beta_k << "\n";
Glog << " -- Nconv = "<< Nconv << "\n";
Glog << " -- Nconv (guess) = "<< Nconv_guess << "\n";
Glog << std::string(74,'*') << std::endl;
} }
@ -412,9 +606,9 @@ private:
//lme[0][L] = beta; //lme[0][L] = beta;
for (int u=0; u<Nu; ++u) { for (int u=0; u<Nu; ++u) {
clog << "norm2(w[" << u << "])= "<< norm2(w[u]) << std::endl; Glog << "norm2(w[" << u << "])= "<< norm2(w[u]) << std::endl;
for (int k=L+u; k<R; ++k) { for (int k=L+u; k<R; ++k) {
clog <<" In block "<< b << ","; Glog <<" In block "<< b << ",";
std::cout <<" beta[" << u << "," << k-L << "] = "; std::cout <<" beta[" << u << "," << k-L << "] = ";
std::cout << lme[u][k] << std::endl; std::cout << lme[u][k] << std::endl;
} }
@ -508,7 +702,7 @@ private:
int Nu, int Nb, int Nk, int Nm, int Nu, int Nb, int Nk, int Nm,
Eigen::MatrixXcd& M) Eigen::MatrixXcd& M)
{ {
//clog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; //Glog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n';
assert( Nk%Nu == 0 && Nm%Nu == 0 ); assert( Nk%Nu == 0 && Nm%Nu == 0 );
assert( Nk <= Nm ); assert( Nk <= Nm );
M = Eigen::MatrixXcd::Zero(Nk,Nk); M = Eigen::MatrixXcd::Zero(Nk,Nk);
@ -526,7 +720,7 @@ private:
M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu]; M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
} }
} }
//clog << "unpackHermitBlockTriDiagMatToEigen() end" << endl; //Glog << "unpackHermitBlockTriDiagMatToEigen() end" << endl;
} }
@ -536,7 +730,7 @@ private:
int Nu, int Nb, int Nk, int Nm, int Nu, int Nb, int Nk, int Nm,
Eigen::MatrixXcd& M) Eigen::MatrixXcd& M)
{ {
//clog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; //Glog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n';
assert( Nk%Nu == 0 && Nm%Nu == 0 ); assert( Nk%Nu == 0 && Nm%Nu == 0 );
assert( Nk <= Nm ); assert( Nk <= Nm );
@ -552,7 +746,7 @@ private:
lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu); lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
} }
} }
//clog << "packHermitBlockTriDiagMatfromEigen() end" << endl; //Glog << "packHermitBlockTriDiagMatfromEigen() end" << endl;
} }
@ -561,7 +755,7 @@ private:
RealD Dsh, RealD Dsh,
Eigen::MatrixXcd& Qprod) Eigen::MatrixXcd& Qprod)
{ {
//clog << "shiftedQRDecompEigen() begin" << '\n'; //Glog << "shiftedQRDecompEigen() begin" << '\n';
Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm); Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
Eigen::MatrixXcd R = Eigen::MatrixXcd::Zero(Nm,Nm); Eigen::MatrixXcd R = Eigen::MatrixXcd::Zero(Nm,Nm);
Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm); Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
@ -652,7 +846,7 @@ private:
// } // }
//} //}
//clog << "shiftedQRDecompEigen() end" << endl; //Glog << "shiftedQRDecompEigen() end" << endl;
} }
void exampleQRDecompEigen(void) void exampleQRDecompEigen(void)
@ -672,74 +866,74 @@ private:
A(2,1) = 24.0; A(2,1) = 24.0;
A(2,2) = -41.0; A(2,2) = -41.0;
clog << "matrix A before ColPivHouseholder" << std::endl; Glog << "matrix A before ColPivHouseholder" << std::endl;
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n'; Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> QRD(A); Eigen::ColPivHouseholderQR<Eigen::MatrixXd> QRD(A);
clog << "matrix A after ColPivHouseholder" << std::endl; Glog << "matrix A after ColPivHouseholder" << std::endl;
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n'; Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
clog << "HouseholderQ with sequence lenth = nonzeroPiviots" << std::endl; Glog << "HouseholderQ with sequence lenth = nonzeroPiviots" << std::endl;
Q = QRD.householderQ().setLength(QRD.nonzeroPivots()); Q = QRD.householderQ().setLength(QRD.nonzeroPivots());
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n'; Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
clog << "HouseholderQ with sequence lenth = 1" << std::endl; Glog << "HouseholderQ with sequence lenth = 1" << std::endl;
Q = QRD.householderQ().setLength(1); Q = QRD.householderQ().setLength(1);
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n'; Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
clog << "HouseholderQ with sequence lenth = 2" << std::endl; Glog << "HouseholderQ with sequence lenth = 2" << std::endl;
Q = QRD.householderQ().setLength(2); Q = QRD.householderQ().setLength(2);
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n'; Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
clog << "matrixR" << std::endl; Glog << "matrixR" << std::endl;
R = QRD.matrixR(); R = QRD.matrixR();
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "R[" << i << "," << j << "] = " << R(i,j) << '\n'; Glog << "R[" << i << "," << j << "] = " << R(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
clog << "rank = " << QRD.rank() << std::endl; Glog << "rank = " << QRD.rank() << std::endl;
clog << "threshold = " << QRD.threshold() << std::endl; Glog << "threshold = " << QRD.threshold() << std::endl;
clog << "matrixP" << std::endl; Glog << "matrixP" << std::endl;
P = QRD.colsPermutation(); P = QRD.colsPermutation();
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "P[" << i << "," << j << "] = " << P(i,j) << '\n'; Glog << "P[" << i << "," << j << "] = " << P(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
clog << "QR decomposition without column pivoting" << std::endl; Glog << "QR decomposition without column pivoting" << std::endl;
A(0,0) = 12.0; A(0,0) = 12.0;
A(0,1) = -51.0; A(0,1) = -51.0;
@ -751,35 +945,35 @@ private:
A(2,1) = 24.0; A(2,1) = 24.0;
A(2,2) = -41.0; A(2,2) = -41.0;
clog << "matrix A before Householder" << std::endl; Glog << "matrix A before Householder" << std::endl;
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n'; Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
Eigen::HouseholderQR<Eigen::MatrixXd> QRDplain(A); Eigen::HouseholderQR<Eigen::MatrixXd> QRDplain(A);
clog << "HouseholderQ" << std::endl; Glog << "HouseholderQ" << std::endl;
Q = QRDplain.householderQ(); Q = QRDplain.householderQ();
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n'; Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
clog << "matrix A after Householder" << std::endl; Glog << "matrix A after Householder" << std::endl;
for ( int i=0; i<3; i++ ) { for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) { for ( int j=0; j<3; j++ ) {
clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n'; Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
} }
} }
clog << std::endl; Glog << std::endl;
} }
}; };
} }
#undef clog #undef Glog
#endif #endif

57
tests/lanczos/Test_dwf_block_lanczos.cc
View File

@ -48,10 +48,13 @@ class CmdJobParams
std::vector<ComplexD> omega; std::vector<ComplexD> omega;
std::vector<ComplexD> boundary_phase; std::vector<ComplexD> boundary_phase;
LanczosType Impl;
int Nu; int Nu;
int Nk; int Nk;
int Np; int Np;
int Nm;
int Nstop; int Nstop;
int Ntest;
int MaxIter; int MaxIter;
double resid; double resid;
@ -61,10 +64,11 @@ class CmdJobParams
CmdJobParams() CmdJobParams()
: gaugefile("Hot"), : gaugefile("Hot"),
Ls(8), mass(0.01), M5(1.8), mob_b(1.5), Ls(8), mass(0.01), M5(1.8), mob_b(1.5),
Nu(4), Nk(200), Np(200), Nstop(100), MaxIter(10), resid(1.0e-8), Impl(LanczosType::irbl),
Nu(4), Nk(200), Np(200), Nstop(100), Ntest(1), MaxIter(10), resid(1.0e-8),
low(0.2), high(5.5), order(11) low(0.2), high(5.5), order(11)
{}; {Nm=Nk+Np;};
void Parse(char **argv, int argc); void Parse(char **argv, int argc);
}; };
@ -155,6 +159,28 @@ void CmdJobParams::Parse(char **argv,int argc)
Np = vi[2]; Np = vi[2];
Nstop = vi[3]; Nstop = vi[3];
MaxIter = vi[4]; MaxIter = vi[4];
// ypj[fixme] mode overriding message is needed.
Impl = LanczosType::irbl;
Nm = Nk+Np;
}
// block Lanczos with explicit extension of its dimensions
if( GridCmdOptionExists(argv,argv+argc,"--rbl") ){
arg = GridCmdOptionPayload(argv,argv+argc,"--rbl");
GridCmdOptionIntVector(arg,vi);
Nu = vi[0];
Nk = vi[1];
Np = vi[2]; // vector space is enlarged by adding Np vectors
Nstop = vi[3];
MaxIter = vi[4];
// ypj[fixme] mode overriding message is needed.
Impl = LanczosType::rbl;
Nm = Nk+Np*MaxIter;
}
if( GridCmdOptionExists(argv,argv+argc,"--check_int") ){
arg = GridCmdOptionPayload(argv,argv+argc,"--check_int");
GridCmdOptionInt(arg,Ntest);
} }
if( GridCmdOptionExists(argv,argv+argc,"--resid") ){ if( GridCmdOptionExists(argv,argv+argc,"--resid") ){
@ -193,7 +219,9 @@ void CmdJobParams::Parse(char **argv,int argc)
std::cout << GridLogMessage <<" Nu "<< Nu << '\n'; std::cout << GridLogMessage <<" Nu "<< Nu << '\n';
std::cout << GridLogMessage <<" Nk "<< Nk << '\n'; std::cout << GridLogMessage <<" Nk "<< Nk << '\n';
std::cout << GridLogMessage <<" Np "<< Np << '\n'; std::cout << GridLogMessage <<" Np "<< Np << '\n';
std::cout << GridLogMessage <<" Nm "<< Nm << '\n';
std::cout << GridLogMessage <<" Nstop "<< Nstop << '\n'; std::cout << GridLogMessage <<" Nstop "<< Nstop << '\n';
std::cout << GridLogMessage <<" Ntest "<< Ntest << '\n';
std::cout << GridLogMessage <<" MaxIter "<< MaxIter << '\n'; std::cout << GridLogMessage <<" MaxIter "<< MaxIter << '\n';
std::cout << GridLogMessage <<" resid "<< resid << '\n'; std::cout << GridLogMessage <<" resid "<< resid << '\n';
std::cout << GridLogMessage <<" Cheby Poly "<< low << "," << high << "," << order << std::endl; std::cout << GridLogMessage <<" Cheby Poly "<< low << "," << high << "," << order << std::endl;
@ -230,6 +258,7 @@ int main (int argc, char ** argv)
} else { } else {
FieldMetaData header; FieldMetaData header;
NerscIO::readConfiguration(Umu,header,JP.gaugefile); NerscIO::readConfiguration(Umu,header,JP.gaugefile);
// ypj [fixme] additional checks for the loaded configuration?
} }
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
@ -238,6 +267,8 @@ int main (int argc, char ** argv)
RealD mass = JP.mass; RealD mass = JP.mass;
RealD M5 = JP.M5; RealD M5 = JP.M5;
// ypj [fixme] flexible support for a various Fermions
// RealD mob_b = JP.mob_b; // Gparity // RealD mob_b = JP.mob_b; // Gparity
// std::vector<ComplexD> omega; // ZMobius // std::vector<ComplexD> omega; // ZMobius
@ -254,10 +285,6 @@ int main (int argc, char ** argv)
ZMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,JP.omega,1.,0.,params); ZMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,JP.omega,1.,0.,params);
SchurDiagTwoOperator<ZMobiusFermionR,FermionField> HermOp(Ddwf); SchurDiagTwoOperator<ZMobiusFermionR,FermionField> HermOp(Ddwf);
int Nu = JP.Nu;
int Nk = JP.Nk;
int Nm = Nk+JP.Np;
//std::vector<double> Coeffs { 0.,-1.}; //std::vector<double> Coeffs { 0.,-1.};
// ypj [note] this may not be supported by some compilers // ypj [note] this may not be supported by some compilers
std::vector<double> Coeffs({ 0.,-1.}); std::vector<double> Coeffs({ 0.,-1.});
@ -267,23 +294,23 @@ int main (int argc, char ** argv)
// Cheb.csv(std::cout); // Cheb.csv(std::cout);
ImplicitlyRestartedBlockLanczos<FermionField> IRBL(HermOp, ImplicitlyRestartedBlockLanczos<FermionField> IRBL(HermOp,
Cheb, Cheb,
JP.Nstop, JP.Nstop, JP.Ntest,
Nu,Nk,Nm, JP.Nu, JP.Nk, JP.Nm,
JP.resid, JP.resid,
JP.MaxIter); JP.MaxIter);
std::vector<RealD> eval(Nm); std::vector<RealD> eval(JP.Nm);
std::vector<FermionField> src(Nu,FrbGrid); std::vector<FermionField> src(JP.Nu,FrbGrid);
for ( int i=0; i<Nu; ++i ) gaussian(RNG5rb,src[i]); for ( int i=0; i<JP.Nu; ++i ) gaussian(RNG5rb,src[i]);
std::vector<FermionField> evec(Nm,FrbGrid); std::vector<FermionField> evec(JP.Nm,FrbGrid);
for(int i=0;i<1;++i){ for(int i=0;i<1;++i){
std::cout << GridLogMessage << i <<" / "<< Nm <<" grid pointer "<< evec[i]._grid << std::endl; std::cout << GridLogMessage << i <<" / "<< JP.Nm <<" grid pointer "<< evec[i]._grid << std::endl;
}; };
int Nconv; int Nconv;
IRBL.calc(eval,evec,src,Nconv); IRBL.calc(eval,evec,src,Nconv,JP.Impl);
Grid_finalize(); Grid_finalize();