portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

FInal for paper

Browse Source
This commit is contained in:
Peter Boyle 2024-07-22 15:26:45 -04:00
parent 486412635a
commit c9d5674d5b
1 changed files with 33 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h
View File

@ -279,11 +279,11 @@ 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);
// Glog << "#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<<"] ";
// std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl; std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
} }
//---------------------------------------------------------------------- //----------------------------------------------------------------------
@ -298,6 +298,7 @@ public:
unpackHermitBlockTriDiagMatToEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM); unpackHermitBlockTriDiagMatToEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
for(int ip=Nk; ip<Nm; ++ip){ for(int ip=Nk; ip<Nm; ++ip){
Glog << " ip "<<ip<<" / "<<Nm<<std::endl;
shiftedQRDecompEigen(BTDM,Nu,Nm,eval2[ip],Q); shiftedQRDecompEigen(BTDM,Nu,Nm,eval2[ip],Q);
} }
@ -325,7 +326,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);
// Glog << "#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<<"] ";
@ -467,10 +468,10 @@ public:
// set initial vector // set initial vector
for (int i=0; i<Nu; ++i) { for (int i=0; i<Nu; ++i) {
// Glog << "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);
// Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl; Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
} }
// exit(-43); // exit(-43);
@ -506,11 +507,11 @@ public:
Qt = Eigen::MatrixXcd::Identity(Nr,Nr); Qt = Eigen::MatrixXcd::Identity(Nr,Nr);
diagonalize(eval2,lmd2,lme2,Nu,Nr,Nr,Qt,grid); diagonalize(eval2,lmd2,lme2,Nu,Nr,Nr,Qt,grid);
_sort.push(eval2,Nr); _sort.push(eval2,Nr);
// Glog << "#Ritz value: "<< std::endl; Glog << "#Ritz value: "<< std::endl;
for(int i=0; i<Nr; ++i){ for(int i=0; i<Nr; ++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<<"] ";
// std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl; std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
} }
// Convergence test // Convergence test
@ -570,6 +571,7 @@ public:
Glog << fname + " NOT converged ; Summary :\n"; Glog << fname + " NOT converged ; Summary :\n";
} else { } else {
Glog << fname + " CONVERGED ; Summary :\n"; Glog << fname + " CONVERGED ; Summary :\n";
Nstop = Nconv_guess; // Just take them all
// Sort convered eigenpairs. // Sort convered eigenpairs.
std::vector<Field> Btmp(Nstop,grid); // waste of space replicating std::vector<Field> Btmp(Nstop,grid); // waste of space replicating
@ -642,7 +644,7 @@ private:
// for (int u=0; u<mrhs; ++u) Glog << " out["<<u<<"] = "<<norm2(out[u])<<std::endl; // for (int u=0; u<mrhs; ++u) Glog << " out["<<u<<"] = "<<norm2(out[u])<<std::endl;
k_start +=mrhs; k_start +=mrhs;
} }
// Glog << "LinAlg "<< std::endl; Glog << "LinAlg "<< std::endl;
if (b>0) { if (b>0) {
for (int u=0; u<Nu; ++u) { for (int u=0; u<Nu; ++u) {
@ -676,7 +678,7 @@ private:
} }
w_copy[u] = w[u]; w_copy[u] = w[u];
} }
// Glog << "LinAlg done"<< std::endl; Glog << "LinAlg done"<< std::endl;
// In block version, the steps 6 and 7 in Lanczos construction is // In block version, the steps 6 and 7 in Lanczos construction is
// replaced by the QR decomposition of new basis block. // replaced by the QR decomposition of new basis block.
@ -689,15 +691,15 @@ private:
} }
// re-orthogonalization for numerical stability // re-orthogonalization for numerical stability
// Glog << "Gram Schmidt"<< std::endl; Glog << "Gram Schmidt"<< std::endl;
orthogonalize(w,Nu,evec,R); orthogonalize(w,Nu,evec,R);
// QR part // QR part
for (int u=1; u<Nu; ++u) { for (int u=1; u<Nu; ++u) {
orthogonalize(w[u],w,u); orthogonalize(w[u],w,u);
} }
// Glog << "Gram Schmidt done "<< std::endl; Glog << "Gram Schmidt done "<< std::endl;
// Glog << "LinAlg "<< std::endl; Glog << "LinAlg "<< std::endl;
for (int u=0; u<Nu; ++u) { for (int u=0; u<Nu; ++u) {
//for (int v=0; v<Nu; ++v) { //for (int v=0; v<Nu; ++v) {
for (int v=u; v<Nu; ++v) { for (int v=u; v<Nu; ++v) {
@ -714,7 +716,7 @@ private:
// Glog <<" In block "<< b << "," <<" beta[" << u << "," << k-L << "] = " << lme[u][k] << std::endl; // Glog <<" In block "<< b << "," <<" beta[" << u << "," << k-L << "] = " << lme[u][k] << std::endl;
} }
} }
// Glog << "LinAlg done "<< std::endl; Glog << "LinAlg done "<< std::endl;
if (b < Nm/Nu-1) { if (b < Nm/Nu-1) {
for (int u=0; u<Nu; ++u) { for (int u=0; u<Nu; ++u) {
@ -933,7 +935,7 @@ if (1){
int Nu, int Nb, int Nk, int Nm, int Nu, int Nb, int Nk, int Nm,
Eigen::MatrixXcd& M) Eigen::MatrixXcd& M)
{ {
//Glog << "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);
@ -951,7 +953,7 @@ if (1){
M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu]; M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
} }
} }
//Glog << "unpackHermitBlockTriDiagMatToEigen() end" << endl; Glog << "unpackHermitBlockTriDiagMatToEigen() end" << std::endl;
} }
@ -961,7 +963,7 @@ if (1){
int Nu, int Nb, int Nk, int Nm, int Nu, int Nb, int Nk, int Nm,
Eigen::MatrixXcd& M) Eigen::MatrixXcd& M)
{ {
//Glog << "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 );
@ -977,7 +979,7 @@ if (1){
lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu); lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
} }
} }
//Glog << "packHermitBlockTriDiagMatfromEigen() end" << endl; Glog << "packHermitBlockTriDiagMatfromEigen() end" <<std::endl;
} }
@ -986,7 +988,7 @@ if (1){
RealD Dsh, RealD Dsh,
Eigen::MatrixXcd& Qprod) Eigen::MatrixXcd& Qprod)
{ {
//Glog << "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);
@ -1002,6 +1004,7 @@ if (1){
// lower triangular part used to represent series // lower triangular part used to represent series
// of Q sequence. // of Q sequence.
Glog << "shiftedQRDecompEigen() Housholder & QR" << '\n';
// equivalent operation of Qprod *= Q // equivalent operation of Qprod *= Q
//M = Eigen::MatrixXcd::Zero(Nm,Nm); //M = Eigen::MatrixXcd::Zero(Nm,Nm);
@ -1022,6 +1025,7 @@ if (1){
Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm); Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
Glog << "shiftedQRDecompEigen() Mtmp create" << '\n';
for (int i=0; i<Nm; ++i) { for (int i=0; i<Nm; ++i) {
for (int j=0; j<Nm-(Nu+1); ++j) { for (int j=0; j<Nm-(Nu+1); ++j) {
for (int k=0; k<Nu+1+j; ++k) { for (int k=0; k<Nu+1+j; ++k) {
@ -1029,6 +1033,7 @@ if (1){
} }
} }
} }
Glog << "shiftedQRDecompEigen() Mtmp loop1" << '\n';
for (int i=0; i<Nm; ++i) { for (int i=0; i<Nm; ++i) {
for (int j=Nm-(Nu+1); j<Nm; ++j) { for (int j=Nm-(Nu+1); j<Nm; ++j) {
for (int k=0; k<Nm; ++k) { for (int k=0; k<Nm; ++k) {
@ -1036,6 +1041,7 @@ if (1){
} }
} }
} }
Glog << "shiftedQRDecompEigen() Mtmp loop2" << '\n';
//static int ntimes = 2; //static int ntimes = 2;
//for (int j=0; j<Nm-(ntimes*Nu); ++j) { //for (int j=0; j<Nm-(ntimes*Nu); ++j) {
@ -1061,11 +1067,13 @@ if (1){
Mtmp(j,i) = conj(Mtmp(i,j)); Mtmp(j,i) = conj(Mtmp(i,j));
} }
} }
Glog << "shiftedQRDecompEigen() Mtmp loop3" << '\n';
for (int i=0; i<Nm; ++i) { for (int i=0; i<Nm; ++i) {
Mtmp(i,i) = real(Mtmp(i,i)) + Dsh; Mtmp(i,i) = real(Mtmp(i,i)) + Dsh;
} }
Glog << "shiftedQRDecompEigen() Mtmp loop4" << '\n';
M = Mtmp; M = Mtmp;
//M = Q.adjoint()*(M*Q); //M = Q.adjoint()*(M*Q);
@ -1077,7 +1085,7 @@ if (1){
// } // }
//} //}
//Glog << "shiftedQRDecompEigen() end" << endl; Glog << "shiftedQRDecompEigen() end" <<std::endl;
} }
void exampleQRDecompEigen(void) void exampleQRDecompEigen(void)