portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 15:55:37 +00:00
Code Releases Activity
b8b5934193
Grid/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h

1213 lines
31 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
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 */
#ifndef GRID_IRL_H
#define GRID_IRL_H
#include <string.h> //memset
#ifdef USE_LAPACK
void LAPACK_dstegr(char *jobz, char *range, int *n, double *d, double *e,
double *vl, double *vu, int *il, int *iu, double *abstol,
int *m, double *w, double *z, int *ldz, int *isuppz,
double *work, int *lwork, int *iwork, int *liwork,
int *info);
#endif
#include <Grid/algorithms/densematrix/DenseMatrix.h>
#include <Grid/algorithms/iterative/EigenSort.h>
namespace Grid {
/////////////////////////////////////////////////////////////
// Implicitly restarted lanczos
/////////////////////////////////////////////////////////////
template<class Field>
class ImplicitlyRestartedLanczos {
const RealD small = 1.0e-16;
public:
int lock;
int get;
int Niter;
int converged;
int Nstop; // Number of evecs checked for convergence
int Nk; // Number of converged sought
int Np; // Np -- Number of spare vecs in kryloc space
int Nm; // Nm -- total number of vectors
RealD eresid;
SortEigen<Field> _sort;
// GridCartesian &_fgrid;
LinearOperatorBase<Field> &_Linop;
OperatorFunction<Field> &_poly;
/////////////////////////
// Constructor
/////////////////////////
void init(void){};
void Abort(int ff, DenseVector<RealD> &evals, DenseVector<DenseVector<RealD> > &evecs);
ImplicitlyRestartedLanczos(
LinearOperatorBase<Field> &Linop, // op
OperatorFunction<Field> & poly, // polynmial
int _Nstop, // sought vecs
int _Nk, // sought vecs
int _Nm, // spare vecs
RealD _eresid, // resid in lmdue deficit
int _Niter) : // Max iterations
_Linop(Linop),
_poly(poly),
Nstop(_Nstop),
Nk(_Nk),
Nm(_Nm),
eresid(_eresid),
Niter(_Niter)
{
Np = Nm-Nk; assert(Np>0);
};
ImplicitlyRestartedLanczos(
LinearOperatorBase<Field> &Linop, // op
OperatorFunction<Field> & poly, // polynmial
int _Nk, // sought vecs
int _Nm, // spare vecs
RealD _eresid, // resid in lmdue deficit
int _Niter) : // Max iterations
_Linop(Linop),
_poly(poly),
Nstop(_Nk),
Nk(_Nk),
Nm(_Nm),
eresid(_eresid),
Niter(_Niter)
{
Np = Nm-Nk; assert(Np>0);
};
/////////////////////////
// Sanity checked this routine (step) against Saad.
/////////////////////////
void RitzMatrix(DenseVector<Field>& evec,int k){
if(1) return;
GridBase *grid = evec[0]._grid;
Field w(grid);
std::cout << "RitzMatrix "<<std::endl;
for(int i=0;i<k;i++){
_poly(_Linop,evec[i],w);
std::cout << "["<<i<<"] ";
for(int j=0;j<k;j++){
ComplexD in = innerProduct(evec[j],w);
if ( fabs((double)i-j)>1 ) {
if (abs(in) >1.0e-9 ) {
std::cout<<"oops"<<std::endl;
abort();
} else
std::cout << " 0 ";
} else {
std::cout << " "<<in<<" ";
}
}
std::cout << std::endl;
}
}
/* Saad PP. 195
1. Choose an initial vector v1 of 2-norm unity. Set β1 ≡ 0, v0 ≡ 0
2. For k = 1,2,...,m Do:
3. wk:=Avk−βkv_{k−1}
4. αk:=(wk,vk) //
5. wk:=wk−αkvk // wk orthog vk
6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop
7. vk+1 := wk/βk+1
8. EndDo
*/
void step(DenseVector<RealD>& lmd,
DenseVector<RealD>& lme,
DenseVector<Field>& evec,
Field& w,int Nm,int k)
{
assert( k< Nm );
_poly(_Linop,evec[k],w); // 3. wk:=Avk−βkv_{k−1}
if(k>0){
w -= lme[k-1] * evec[k-1];
}
ComplexD zalph = innerProduct(evec[k],w); // 4. αk:=(wk,vk)
RealD alph = real(zalph);
w = w - alph * evec[k];// 5. wk:=wk−αkvk
RealD beta = normalise(w); // 6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop
// 7. vk+1 := wk/βk+1
// std::cout << "alpha = " << zalph << " beta "<<beta<<std::endl;
const RealD tiny = 1.0e-20;
if ( beta < tiny ) {
std::cout << " beta is tiny "<<beta<<std::endl;
}
lmd[k] = alph;
lme[k] = beta;
if (k>0) {
orthogonalize(w,evec,k); // orthonormalise
}
if(k < Nm-1) evec[k+1] = w;
}
void qr_decomp(DenseVector<RealD>& lmd,
DenseVector<RealD>& lme,
int Nk,
int Nm,
DenseVector<RealD>& Qt,
RealD Dsh,
int kmin,
int kmax)
{
int k = kmin-1;
RealD x;
RealD Fden = 1.0/hypot(lmd[k]-Dsh,lme[k]);
RealD c = ( lmd[k] -Dsh) *Fden;
RealD s = -lme[k] *Fden;
RealD tmpa1 = lmd[k];
RealD tmpa2 = lmd[k+1];
RealD tmpb = lme[k];
lmd[k] = c*c*tmpa1 +s*s*tmpa2 -2.0*c*s*tmpb;
lmd[k+1] = s*s*tmpa1 +c*c*tmpa2 +2.0*c*s*tmpb;
lme[k] = c*s*(tmpa1-tmpa2) +(c*c-s*s)*tmpb;
x =-s*lme[k+1];
lme[k+1] = c*lme[k+1];
for(int i=0; i<Nk; ++i){
RealD Qtmp1 = Qt[i+Nm*k ];
RealD Qtmp2 = Qt[i+Nm*(k+1)];
Qt[i+Nm*k ] = c*Qtmp1 - s*Qtmp2;
Qt[i+Nm*(k+1)] = s*Qtmp1 + c*Qtmp2;
}
// Givens transformations
for(int k = kmin; k < kmax-1; ++k){
RealD Fden = 1.0/hypot(x,lme[k-1]);
RealD c = lme[k-1]*Fden;
RealD s = - x*Fden;
RealD tmpa1 = lmd[k];
RealD tmpa2 = lmd[k+1];
RealD tmpb = lme[k];
lmd[k] = c*c*tmpa1 +s*s*tmpa2 -2.0*c*s*tmpb;
lmd[k+1] = s*s*tmpa1 +c*c*tmpa2 +2.0*c*s*tmpb;
lme[k] = c*s*(tmpa1-tmpa2) +(c*c-s*s)*tmpb;
lme[k-1] = c*lme[k-1] -s*x;
if(k != kmax-2){
x = -s*lme[k+1];
lme[k+1] = c*lme[k+1];
}
for(int i=0; i<Nk; ++i){
RealD Qtmp1 = Qt[i+Nm*k ];
RealD Qtmp2 = Qt[i+Nm*(k+1)];
Qt[i+Nm*k ] = c*Qtmp1 -s*Qtmp2;
Qt[i+Nm*(k+1)] = s*Qtmp1 +c*Qtmp2;
}
}
}
#ifdef USE_LAPACK
void diagonalize_lapack(DenseVector<RealD>& lmd,
DenseVector<RealD>& lme,
int N1,
int N2,
DenseVector<RealD>& Qt,
GridBase *grid){
const int size = Nm;
// tevals.resize(size);
// tevecs.resize(size);
int NN = N1;
double evals_tmp[NN];
double evec_tmp[NN][NN];
memset(evec_tmp[0],0,sizeof(double)*NN*NN);
// double AA[NN][NN];
double DD[NN];
double EE[NN];
for (int i = 0; i< NN; i++)
for (int j = i - 1; j <= i + 1; j++)
if ( j < NN && j >= 0 ) {
if (i==j) DD[i] = lmd[i];
if (i==j) evals_tmp[i] = lmd[i];
if (j==(i-1)) EE[j] = lme[j];
}
int evals_found;
int lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ;
int liwork = 3+NN*10 ;
int iwork[liwork];
double work[lwork];
int isuppz[2*NN];
char jobz = 'V'; // calculate evals & evecs
char range = 'I'; // calculate all evals
// char range = 'A'; // calculate all evals
char uplo = 'U'; // refer to upper half of original matrix
char compz = 'I'; // Compute eigenvectors of tridiagonal matrix
int ifail[NN];
int info;
// int total = QMP_get_number_of_nodes();
// int node = QMP_get_node_number();
// GridBase *grid = evec[0]._grid;
int total = grid->_Nprocessors;
int node = grid->_processor;
int interval = (NN/total)+1;
double vl = 0.0, vu = 0.0;
int il = interval*node+1 , iu = interval*(node+1);
if (iu > NN) iu=NN;
double tol = 0.0;
if (1) {
memset(evals_tmp,0,sizeof(double)*NN);
if ( il <= NN){
printf("total=%d node=%d il=%d iu=%d\n",total,node,il,iu);
LAPACK_dstegr(&jobz, &range, &NN,
(double*)DD, (double*)EE,
&vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A'
&tol, // tolerance
&evals_found, evals_tmp, (double*)evec_tmp, &NN,
isuppz,
work, &lwork, iwork, &liwork,
&info);
for (int i = iu-1; i>= il-1; i--){
printf("node=%d evals_found=%d evals_tmp[%d] = %g\n",node,evals_found, i - (il-1),evals_tmp[i - (il-1)]);
evals_tmp[i] = evals_tmp[i - (il-1)];
if (il>1) evals_tmp[i-(il-1)]=0.;
for (int j = 0; j< NN; j++){
evec_tmp[i][j] = evec_tmp[i - (il-1)][j];
if (il>1) evec_tmp[i-(il-1)][j]=0.;
}
}
}
{
// QMP_sum_double_array(evals_tmp,NN);
// QMP_sum_double_array((double *)evec_tmp,NN*NN);
grid->GlobalSumVector(evals_tmp,NN);
grid->GlobalSumVector((double*)evec_tmp,NN*NN);
}
}
// cheating a bit. It is better to sort instead of just reversing it, but the document of the routine says evals are sorted in increasing order. qr gives evals in decreasing order.
for(int i=0;i<NN;i++){
for(int j=0;j<NN;j++)
Qt[(NN-1-i)*N2+j]=evec_tmp[i][j];
lmd [NN-1-i]=evals_tmp[i];
}
}
#endif
void diagonalize(DenseVector<RealD>& lmd,
DenseVector<RealD>& lme,
int N2,
int N1,
DenseVector<RealD>& Qt,
GridBase *grid)
{
#ifdef USE_LAPACK
const int check_lapack=0; // just use lapack if 0, check against lapack if 1
if(!check_lapack)
return diagonalize_lapack(lmd,lme,N2,N1,Qt,grid);
DenseVector <RealD> lmd2(N1);
DenseVector <RealD> lme2(N1);
DenseVector<RealD> Qt2(N1*N1);
for(int k=0; k<N1; ++k){
lmd2[k] = lmd[k];
lme2[k] = lme[k];
}
for(int k=0; k<N1*N1; ++k)
Qt2[k] = Qt[k];
// diagonalize_lapack(lmd2,lme2,Nm2,Nm,Qt,grid);
#endif
int Niter = 100*N1;
int kmin = 1;
int kmax = N2;
// (this should be more sophisticated)
for(int iter=0; iter<Niter; ++iter){
// determination of 2x2 leading submatrix
RealD dsub = lmd[kmax-1]-lmd[kmax-2];
RealD dd = sqrt(dsub*dsub + 4.0*lme[kmax-2]*lme[kmax-2]);
RealD Dsh = 0.5*(lmd[kmax-2]+lmd[kmax-1] +dd*(dsub/fabs(dsub)));
// (Dsh: shift)
// transformation
qr_decomp(lmd,lme,N2,N1,Qt,Dsh,kmin,kmax);
// Convergence criterion (redef of kmin and kamx)
for(int j=kmax-1; j>= kmin; --j){
RealD dds = fabs(lmd[j-1])+fabs(lmd[j]);
if(fabs(lme[j-1])+dds > dds){
kmax = j+1;
goto continued;
}
}
Niter = iter;
#ifdef USE_LAPACK
if(check_lapack){
const double SMALL=1e-8;
diagonalize_lapack(lmd2,lme2,N2,N1,Qt2,grid);
DenseVector <RealD> lmd3(N2);
for(int k=0; k<N2; ++k) lmd3[k]=lmd[k];
_sort.push(lmd3,N2);
_sort.push(lmd2,N2);
for(int k=0; k<N2; ++k){
if (fabs(lmd2[k] - lmd3[k]) >SMALL) std::cout <<"lmd(qr) lmd(lapack) "<< k << ": " << lmd2[k] <<" "<< lmd3[k] <<std::endl;
// if (fabs(lme2[k] - lme[k]) >SMALL) std::cout <<"lme(qr)-lme(lapack) "<< k << ": " << lme2[k] - lme[k] <<std::endl;
}
for(int k=0; k<N1*N1; ++k){
// if (fabs(Qt2[k] - Qt[k]) >SMALL) std::cout <<"Qt(qr)-Qt(lapack) "<< k << ": " << Qt2[k] - Qt[k] <<std::endl;
}
}
#endif
return;
continued:
for(int j=0; j<kmax-1; ++j){
RealD dds = fabs(lmd[j])+fabs(lmd[j+1]);
if(fabs(lme[j])+dds > dds){
kmin = j+1;
break;
}
}
}
std::cout << "[QL method] Error - Too many iteration: "<<Niter<<"\n";
abort();
}
#if 1
static RealD normalise(Field& v)
{
RealD nn = norm2(v);
nn = sqrt(nn);
v = v * (1.0/nn);
return nn;
}
void orthogonalize(Field& w,
DenseVector<Field>& evec,
int k)
{
typedef typename Field::scalar_type MyComplex;
MyComplex ip;
if ( 0 ) {
for(int j=0; j<k; ++j){
normalise(evec[j]);
for(int i=0;i<j;i++){
ip = innerProduct(evec[i],evec[j]); // are the evecs normalised? ; this assumes so.
evec[j] = evec[j] - ip *evec[i];
}
}
}
for(int j=0; j<k; ++j){
ip = innerProduct(evec[j],w); // are the evecs normalised? ; this assumes so.
w = w - ip * evec[j];
}
normalise(w);
}
void setUnit_Qt(int Nm, DenseVector<RealD> &Qt) {
for(int i=0; i<Qt.size(); ++i) Qt[i] = 0.0;
for(int k=0; k<Nm; ++k) Qt[k + k*Nm] = 1.0;
}
/* Rudy Arthur's thesis pp.137
------------------------
Require: M > K P = M − K †
Compute the factorization AVM = VM HM + fM eM
repeat
Q=I
for i = 1,...,P do
QiRi =HM −θiI Q = QQi
H M = Q †i H M Q i
end for
βK =HM(K+1,K) σK =Q(M,K)
r=vK+1βK +rσK
VK =VM(1:M)Q(1:M,1:K)
HK =HM(1:K,1:K)
→AVK =VKHK +fKe†K † Extend to an M = K + P step factorization AVM = VMHM + fMeM
until convergence
*/
void calc(DenseVector<RealD>& eval,
DenseVector<Field>& evec,
const Field& src,
int& Nconv)
{
GridBase *grid = evec[0]._grid;
assert(grid == src._grid);
std::cout << " -- Nk = " << Nk << " Np = "<< Np << std::endl;
std::cout << " -- Nm = " << Nm << std::endl;
std::cout << " -- size of eval = " << eval.size() << std::endl;
std::cout << " -- size of evec = " << evec.size() << std::endl;
assert(Nm == evec.size() && Nm == eval.size());
DenseVector<RealD> lme(Nm);
DenseVector<RealD> lme2(Nm);
DenseVector<RealD> eval2(Nm);
DenseVector<RealD> Qt(Nm*Nm);
DenseVector<int> Iconv(Nm);
DenseVector<Field> B(Nm,grid); // waste of space replicating
Field f(grid);
Field v(grid);
int k1 = 1;
int k2 = Nk;
Nconv = 0;
RealD beta_k;
// Set initial vector
// (uniform vector) Why not src??
// evec[0] = 1.0;
evec[0] = src;
std:: cout <<"norm2(src)= " << norm2(src)<<std::endl;
// << src._grid << std::endl;
normalise(evec[0]);
std:: cout <<"norm2(evec[0])= " << norm2(evec[0]) <<std::endl;
// << evec[0]._grid << std::endl;
// Initial Nk steps
for(int k=0; k<Nk; ++k) step(eval,lme,evec,f,Nm,k);
// std:: cout <<"norm2(evec[1])= " << norm2(evec[1]) << std::endl;
// std:: cout <<"norm2(evec[2])= " << norm2(evec[2]) << std::endl;
RitzMatrix(evec,Nk);
for(int k=0; k<Nk; ++k){
// std:: cout <<"eval " << k << " " <<eval[k] << std::endl;
// std:: cout <<"lme " << k << " " << lme[k] << std::endl;
}
// Restarting loop begins
for(int iter = 0; iter<Niter; ++iter){
std::cout<<"\n Restart iteration = "<< iter << std::endl;
//
// Rudy does a sort first which looks very different. Getting fed up with sorting out the algo defs.
// We loop over
//
for(int k=Nk; k<Nm; ++k) step(eval,lme,evec,f,Nm,k);
f *= lme[Nm-1];
RitzMatrix(evec,k2);
// getting eigenvalues
for(int k=0; k<Nm; ++k){
eval2[k] = eval[k+k1-1];
lme2[k] = lme[k+k1-1];
}
setUnit_Qt(Nm,Qt);
diagonalize(eval2,lme2,Nm,Nm,Qt,grid);
// sorting
_sort.push(eval2,Nm);
// Implicitly shifted QR transformations
setUnit_Qt(Nm,Qt);
for(int ip=k2; ip<Nm; ++ip){
std::cout << "qr_decomp "<< ip << " "<< eval2[ip] << std::endl;
qr_decomp(eval,lme,Nm,Nm,Qt,eval2[ip],k1,Nm);
}
for(int i=0; i<(Nk+1); ++i) B[i] = 0.0;
for(int j=k1-1; j<k2+1; ++j){
for(int k=0; k<Nm; ++k){
B[j].checkerboard = evec[k].checkerboard;
B[j] += Qt[k+Nm*j] * evec[k];
}
}
for(int j=k1-1; j<k2+1; ++j) evec[j] = B[j];
// Compressed vector f and beta(k2)
f *= Qt[Nm-1+Nm*(k2-1)];
f += lme[k2-1] * evec[k2];
beta_k = norm2(f);
beta_k = sqrt(beta_k);
std::cout<<" beta(k) = "<<beta_k<<std::endl;
RealD betar = 1.0/beta_k;
evec[k2] = betar * f;
lme[k2-1] = beta_k;
// Convergence test
for(int k=0; k<Nm; ++k){
eval2[k] = eval[k];
lme2[k] = lme[k];
}
setUnit_Qt(Nm,Qt);
diagonalize(eval2,lme2,Nk,Nm,Qt,grid);
for(int k = 0; k<Nk; ++k) B[k]=0.0;
for(int j = 0; j<Nk; ++j){
for(int k = 0; k<Nk; ++k){
B[j].checkerboard = evec[k].checkerboard;
B[j] += Qt[k+j*Nm] * evec[k];
}
// std::cout << "norm(B["<<j<<"])="<<norm2(B[j])<<std::endl;
}
// _sort.push(eval2,B,Nk);
Nconv = 0;
// std::cout << std::setiosflags(std::ios_base::scientific);
for(int i=0; i<Nk; ++i){
// _poly(_Linop,B[i],v);
_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);
std::cout.precision(13);
std::cout << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
std::cout << "eval = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2[i];
std::cout <<" |H B[i] - eval[i]B[i]|^2 "<< std::setw(25)<< 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) ){
Iconv[Nconv] = i;
++Nconv;
}
} // i-loop end
// std::cout << std::resetiosflags(std::ios_base::scientific);
std::cout<<" #modes converged: "<<Nconv<<std::endl;
if( Nconv>=Nstop ){
goto converged;
}
} // end of iter loop
std::cout<<"\n NOT converged.\n";
abort();
converged:
// Sorting
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);
std::cout << "\n Converged\n Summary :\n";
std::cout << " -- Iterations = "<< Nconv << "\n";
std::cout << " -- beta(k) = "<< beta_k << "\n";
std::cout << " -- Nconv = "<< Nconv << "\n";
}
/////////////////////////////////////////////////
// Adapted from Rudy's lanczos factor routine
/////////////////////////////////////////////////
int Lanczos_Factor(int start, int end, int cont,
DenseVector<Field> & bq,
Field &bf,
DenseMatrix<RealD> &H){
GridBase *grid = bq[0]._grid;
RealD beta;
RealD sqbt;
RealD alpha;
for(int i=start;i<Nm;i++){
for(int j=start;j<Nm;j++){
H[i][j]=0.0;
}
}
std::cout<<"Lanczos_Factor start/end " <<start <<"/"<<end<<std::endl;
// Starting from scratch, bq[0] contains a random vector and |bq[0]| = 1
int first;
if(start == 0){
std::cout << "start == 0\n"; //TESTING
_poly(_Linop,bq[0],bf);
alpha = real(innerProduct(bq[0],bf));//alpha = bq[0]^dag A bq[0]
std::cout << "alpha = " << alpha << std::endl;
bf = bf - alpha * bq[0]; //bf = A bq[0] - alpha bq[0]
H[0][0]=alpha;
std::cout << "Set H(0,0) to " << H[0][0] << std::endl;
first = 1;
} else {
first = start;
}
// I think start==0 and cont==zero are the same. Test this
// If so I can drop "cont" parameter?
if( cont ) assert(start!=0);
if( start==0 ) assert(cont!=0);
if( cont){
beta = 0;sqbt = 0;
std::cout << "cont is true so setting beta to zero\n";
} else {
beta = norm2(bf);
sqbt = sqrt(beta);
std::cout << "beta = " << beta << std::endl;
}
for(int j=first;j<end;j++){
std::cout << "Factor j " << j <<std::endl;
if(cont){ // switches to factoring; understand start!=0 and initial bf value is right.
bq[j] = bf; cont = false;
}else{
bq[j] = (1.0/sqbt)*bf ;
H[j][j-1]=H[j-1][j] = sqbt;
}
_poly(_Linop,bq[j],bf);
bf = bf - (1.0/sqbt)*bq[j-1]; //bf = A bq[j] - beta bq[j-1] // PAB this comment was incorrect in beta term??
alpha = real(innerProduct(bq[j],bf)); //alpha = bq[j]^dag A bq[j]
bf = bf - alpha*bq[j]; //bf = A bq[j] - beta bq[j-1] - alpha bq[j]
RealD fnorm = norm2(bf);
RealD bck = sqrt( real( conjugate(alpha)*alpha ) + beta );
beta = fnorm;
sqbt = sqrt(beta);
std::cout << "alpha = " << alpha << " fnorm = " << fnorm << '\n';
///Iterative refinement of orthogonality V = [ bq[0] bq[1] ... bq[M] ]
int re = 0;
// FIXME undefined params; how set in Rudy's code
int ref =0;
Real rho = 1.0e-8;
while( re == ref || (sqbt < rho * bck && re < 5) ){
Field tmp2(grid);
Field tmp1(grid);
//bex = V^dag bf
DenseVector<ComplexD> bex(j+1);
for(int k=0;k<j+1;k++){
bex[k] = innerProduct(bq[k],bf);
}
zero_fermion(tmp2);
//tmp2 = V s
for(int l=0;l<j+1;l++){
RealD nrm = norm2(bq[l]);
axpy(tmp1,0.0,bq[l],bq[l]); scale(tmp1,bex[l]); //tmp1 = V[j] bex[j]
axpy(tmp2,1.0,tmp2,tmp1); //tmp2 += V[j] bex[j]
}
//bf = bf - V V^dag bf. Subtracting off any component in span { V[j] }
RealD btc = axpy_norm(bf,-1.0,tmp2,bf);
alpha = alpha + real(bex[j]); sqbt = sqrt(real(btc));
// FIXME is alpha real in RUDY's code?
RealD nmbex = 0;for(int k=0;k<j+1;k++){nmbex = nmbex + real( conjugate(bex[k])*bex[k] );}
bck = sqrt( nmbex );
re++;
}
std::cout << "Iteratively refined orthogonality, changes alpha\n";
if(re > 1) std::cout << "orthagonality refined " << re << " times" <<std::endl;
H[j][j]=alpha;
}
return end;
}
void EigenSort(DenseVector<double> evals,
DenseVector<Field> evecs){
int N= evals.size();
_sort.push(evals,evecs, evals.size(),N);
}
void ImplicitRestart(int TM, DenseVector<RealD> &evals, DenseVector<DenseVector<RealD> > &evecs, DenseVector<Field> &bq, Field &bf, int cont)
{
std::cout << "ImplicitRestart begin. Eigensort starting\n";
DenseMatrix<RealD> H; Resize(H,Nm,Nm);
EigenSort(evals, evecs);
///Assign shifts
int K=Nk;
int M=Nm;
int P=Np;
int converged=0;
if(K - converged < 4) P = (M - K-1); //one
// DenseVector<RealD> shifts(P + shift_extra.size());
DenseVector<RealD> shifts(P);
for(int k = 0; k < P; ++k)
shifts[k] = evals[k];
/// Shift to form a new H and q
DenseMatrix<RealD> Q; Resize(Q,TM,TM);
Unity(Q);
Shift(Q, shifts); // H is implicitly passed in in Rudy's Shift routine
int ff = K;
/// Shifted H defines a new K step Arnoldi factorization
RealD beta = H[ff][ff-1];
RealD sig = Q[TM - 1][ff - 1];
std::cout << "beta = " << beta << " sig = " << real(sig) <<std::endl;
std::cout << "TM = " << TM << " ";
std::cout << norm2(bq[0]) << " -- before" <<std::endl;
/// q -> q Q
times_real(bq, Q, TM);
std::cout << norm2(bq[0]) << " -- after " << ff <<std::endl;
bf = beta* bq[ff] + sig* bf;
/// Do the rest of the factorization
ff = Lanczos_Factor(ff, M,cont,bq,bf,H);
if(ff < M)
Abort(ff, evals, evecs);
}
///Run the Eigensolver
void Run(int cont, DenseVector<Field> &bq, Field &bf, DenseVector<DenseVector<RealD> > & evecs,DenseVector<RealD> &evals)
{
init();
int M=Nm;
DenseMatrix<RealD> H; Resize(H,Nm,Nm);
Resize(evals,Nm);
Resize(evecs,Nm);
int ff = Lanczos_Factor(0, M, cont, bq,bf,H); // 0--M to begin with
if(ff < M) {
std::cout << "Krylov: aborting ff "<<ff <<" "<<M<<std::endl;
abort(); // Why would this happen?
}
int itcount = 0;
bool stop = false;
for(int it = 0; it < Niter && (converged < Nk); ++it) {
std::cout << "Krylov: Iteration --> " << it << std::endl;
int lock_num = lock ? converged : 0;
DenseVector<RealD> tevals(M - lock_num );
DenseMatrix<RealD> tevecs; Resize(tevecs,M - lock_num,M - lock_num);
//check residual of polynominal
TestConv(H,M, tevals, tevecs);
if(converged >= Nk)
break;
ImplicitRestart(ff, tevals,tevecs,H);
}
Wilkinson<RealD>(H, evals, evecs, small);
// Check();
std::cout << "Done "<<std::endl;
}
///H - shift I = QR; H = Q* H Q
void Shift(DenseMatrix<RealD> & H,DenseMatrix<RealD> &Q, DenseVector<RealD> shifts) {
int P; Size(shifts,P);
int M; SizeSquare(Q,M);
Unity(Q);
int lock_num = lock ? converged : 0;
RealD t_Househoulder_vector(0.0);
RealD t_Househoulder_mult(0.0);
for(int i=0;i<P;i++){
RealD x, y, z;
DenseVector<RealD> ck(3), v(3);
x = H[lock_num+0][lock_num+0]-shifts[i];
y = H[lock_num+1][lock_num+0];
ck[0] = x; ck[1] = y; ck[2] = 0;
normalise(ck); ///Normalization cancels in PHP anyway
RealD beta;
Householder_vector<RealD>(ck, 0, 2, v, beta);
Householder_mult<RealD>(H,v,beta,0,lock_num+0,lock_num+2,0);
Householder_mult<RealD>(H,v,beta,0,lock_num+0,lock_num+2,1);
///Accumulate eigenvector
Householder_mult<RealD>(Q,v,beta,0,lock_num+0,lock_num+2,1);
int sw = 0;
for(int k=lock_num+0;k<M-2;k++){
x = H[k+1][k];
y = H[k+2][k];
z = (RealD)0.0;
if(k+3 <= M-1){
z = H[k+3][k];
}else{
sw = 1; v[2] = 0.0;
}
ck[0] = x; ck[1] = y; ck[2] = z;
normalise(ck);
Householder_vector<RealD>(ck, 0, 2-sw, v, beta);
Householder_mult<RealD>(H,v, beta,0,k+1,k+3-sw,0);
Householder_mult<RealD>(H,v, beta,0,k+1,k+3-sw,1);
///Accumulate eigenvector
Householder_mult<RealD>(Q,v, beta,0,k+1,k+3-sw,1);
}
}
}
void TestConv(DenseMatrix<RealD> & H,int SS,
DenseVector<Field> &bq, Field &bf,
DenseVector<RealD> &tevals, DenseVector<DenseVector<RealD> > &tevecs,
int lock, int converged)
{
std::cout << "Converged " << converged << " so far." << std::endl;
int lock_num = lock ? converged : 0;
int M = Nm;
///Active Factorization
DenseMatrix<RealD> AH; Resize(AH,SS - lock_num,SS - lock_num );
AH = GetSubMtx(H,lock_num, SS, lock_num, SS);
int NN=tevals.size();
int AHsize=SS-lock_num;
RealD small=1.0e-16;
Wilkinson<RealD>(AH, tevals, tevecs, small);
EigenSort(tevals, tevecs);
RealD resid_nrm= norm2(bf);
if(!lock) converged = 0;
#if 0
for(int i = SS - lock_num - 1; i >= SS - Nk && i >= 0; --i){
RealD diff = 0;
diff = abs( tevecs[i][Nm - 1 - lock_num] ) * resid_nrm;
std::cout << "residual estimate " << SS-1-i << " " << diff << " of (" << tevals[i] << ")" << std::endl;
if(diff < converged) {
if(lock) {
DenseMatrix<RealD> Q; Resize(Q,M,M);
bool herm = true;
Lock(H, Q, tevals[i], converged, small, SS, herm);
times_real(bq, Q, bq.size());
bf = Q[M - 1][M - 1]* bf;
lock_num++;
}
converged++;
std::cout << " converged on eval " << converged << " of " << Nk << std::endl;
} else {
break;
}
}
#endif
std::cout << "Got " << converged << " so far " <<std::endl;
}
///Check
void Check(DenseVector<RealD> &evals,
DenseVector<DenseVector<RealD> > &evecs) {
DenseVector<RealD> goodval(this->get);
EigenSort(evals,evecs);
int NM = Nm;
DenseVector< DenseVector<RealD> > V; Size(V,NM);
DenseVector<RealD> QZ(NM*NM);
for(int i = 0; i < NM; i++){
for(int j = 0; j < NM; j++){
// evecs[i][j];
}
}
}
/**
There is some matrix Q such that for any vector y
Q.e_1 = y and Q is unitary.
**/
template<class T>
static T orthQ(DenseMatrix<T> &Q, DenseVector<T> y){
int N = y.size(); //Matrix Size
Fill(Q,0.0);
T tau;
for(int i=0;i<N;i++){
Q[i][0]=y[i];
}
T sig = conj(y[0])*y[0];
T tau0 = abs(sqrt(sig));
for(int j=1;j<N;j++){
sig += conj(y[j])*y[j];
tau = abs(sqrt(sig) );
if(abs(tau0) > 0.0){
T gam = conj( (y[j]/tau)/tau0 );
for(int k=0;k<=j-1;k++){
Q[k][j]=-gam*y[k];
}
Q[j][j]=tau0/tau;
} else {
Q[j-1][j]=1.0;
}
tau0 = tau;
}
return tau;
}
/**
There is some matrix Q such that for any vector y
Q.e_k = y and Q is unitary.
**/
template< class T>
static T orthU(DenseMatrix<T> &Q, DenseVector<T> y){
T tau = orthQ(Q,y);
SL(Q);
return tau;
}
/**
Wind up with a matrix with the first con rows untouched
say con = 2
Q is such that Qdag H Q has {x, x, val, 0, 0, 0, 0, ...} as 1st colum
and the matrix is upper hessenberg
and with f and Q appropriately modidied with Q is the arnoldi factorization
**/
template<class T>
static void Lock(DenseMatrix<T> &H, // Hess mtx
DenseMatrix<T> &Q, // Lock Transform
T val, // value to be locked
int con, // number already locked
RealD small,
int dfg,
bool herm)
{
//ForceTridiagonal(H);
int M = H.dim;
DenseVector<T> vec; Resize(vec,M-con);
DenseMatrix<T> AH; Resize(AH,M-con,M-con);
AH = GetSubMtx(H,con, M, con, M);
DenseMatrix<T> QQ; Resize(QQ,M-con,M-con);
Unity(Q); Unity(QQ);
DenseVector<T> evals; Resize(evals,M-con);
DenseMatrix<T> evecs; Resize(evecs,M-con,M-con);
Wilkinson<T>(AH, evals, evecs, small);
int k=0;
RealD cold = abs( val - evals[k]);
for(int i=1;i<M-con;i++){
RealD cnew = abs( val - evals[i]);
if( cnew < cold ){k = i; cold = cnew;}
}
vec = evecs[k];
ComplexD tau;
orthQ(QQ,vec);
//orthQM(QQ,AH,vec);
AH = Hermitian(QQ)*AH;
AH = AH*QQ;
for(int i=con;i<M;i++){
for(int j=con;j<M;j++){
Q[i][j]=QQ[i-con][j-con];
H[i][j]=AH[i-con][j-con];
}
}
for(int j = M-1; j>con+2; j--){
DenseMatrix<T> U; Resize(U,j-1-con,j-1-con);
DenseVector<T> z; Resize(z,j-1-con);
T nm = norm(z);
for(int k = con+0;k<j-1;k++){
z[k-con] = conj( H(j,k+1) );
}
normalise(z);
RealD tmp = 0;
for(int i=0;i<z.size()-1;i++){tmp = tmp + abs(z[i]);}
if(tmp < small/( (RealD)z.size()-1.0) ){ continue;}
tau = orthU(U,z);
DenseMatrix<T> Hb; Resize(Hb,j-1-con,M);
for(int a = 0;a<M;a++){
for(int b = 0;b<j-1-con;b++){
T sum = 0;
for(int c = 0;c<j-1-con;c++){
sum += H[a][con+1+c]*U[c][b];
}//sum += H(a,con+1+c)*U(c,b);}
Hb[b][a] = sum;
}
}
for(int k=con+1;k<j;k++){
for(int l=0;l<M;l++){
H[l][k] = Hb[k-1-con][l];
}
}//H(Hb[k-1-con][l] , l,k);}}
DenseMatrix<T> Qb; Resize(Qb,M,M);
for(int a = 0;a<M;a++){
for(int b = 0;b<j-1-con;b++){
T sum = 0;
for(int c = 0;c<j-1-con;c++){
sum += Q[a][con+1+c]*U[c][b];
}//sum += Q(a,con+1+c)*U(c,b);}
Qb[b][a] = sum;
}
}
for(int k=con+1;k<j;k++){
for(int l=0;l<M;l++){
Q[l][k] = Qb[k-1-con][l];
}
}//Q(Qb[k-1-con][l] , l,k);}}
DenseMatrix<T> Hc; Resize(Hc,M,M);
for(int a = 0;a<j-1-con;a++){
for(int b = 0;b<M;b++){
T sum = 0;
for(int c = 0;c<j-1-con;c++){
sum += conj( U[c][a] )*H[con+1+c][b];
}//sum += conj( U(c,a) )*H(con+1+c,b);}
Hc[b][a] = sum;
}
}
for(int k=0;k<M;k++){
for(int l=con+1;l<j;l++){
H[l][k] = Hc[k][l-1-con];
}
}//H(Hc[k][l-1-con] , l,k);}}
}
}
#endif
};
}
#endif
View Git Blame Copy Permalink