mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-15 02:05:37 +00:00
Adding SimpleLanczos
This commit is contained in:
Chulwoo Jung
parent
8313367a50
commit
05d04ceff8
776
lib/algorithms/iterative/SimpleLanczos.h
Normal file
776
lib/algorithms/iterative/SimpleLanczos.h
Normal file
@ -0,0 +1,776 @@
|
||||
/*************************************************************************************
|
||||
|
||||
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>
|
||||
Author: Chulwoo Jung <chulwoo@bnl.gov>
|
||||
|
||||
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_LANC_H
|
||||
#define GRID_LANC_H
|
||||
|
||||
#include <string.h> //memset
|
||||
|
||||
#ifdef USE_LAPACK
|
||||
#ifdef USE_MKL
|
||||
#include<mkl_lapack.h>
|
||||
#else
|
||||
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);
|
||||
//#include <lapacke/lapacke.h>
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#include <Grid/algorithms/densematrix/DenseMatrix.h>
|
||||
#include <Grid/algorithms/iterative/EigenSort.h>
|
||||
|
||||
// eliminate temorary vector in calc()
|
||||
#define MEM_SAVE
|
||||
|
||||
namespace Grid {
|
||||
|
||||
struct Bisection {
|
||||
|
||||
#if 0
|
||||
static void get_eig2(int row_num,std::vector<RealD> &ALPHA,std::vector<RealD> &BETA, std::vector<RealD> & eig)
|
||||
{
|
||||
int i,j;
|
||||
std::vector<RealD> evec1(row_num+3);
|
||||
std::vector<RealD> evec2(row_num+3);
|
||||
RealD eps2;
|
||||
ALPHA[1]=0.;
|
||||
BETHA[1]=0.;
|
||||
for(i=0;i<row_num-1;i++) {
|
||||
ALPHA[i+1] = A[i*(row_num+1)].real();
|
||||
BETHA[i+2] = A[i*(row_num+1)+1].real();
|
||||
}
|
||||
ALPHA[row_num] = A[(row_num-1)*(row_num+1)].real();
|
||||
bisec(ALPHA,BETHA,row_num,1,row_num,1e-10,1e-10,evec1,eps2);
|
||||
bisec(ALPHA,BETHA,row_num,1,row_num,1e-16,1e-16,evec2,eps2);
|
||||
|
||||
// Do we really need to sort here?
|
||||
int begin=1;
|
||||
int end = row_num;
|
||||
int swapped=1;
|
||||
while(swapped) {
|
||||
swapped=0;
|
||||
for(i=begin;i<end;i++){
|
||||
if(mag(evec2[i])>mag(evec2[i+1])) {
|
||||
swap(evec2+i,evec2+i+1);
|
||||
swapped=1;
|
||||
}
|
||||
}
|
||||
end--;
|
||||
for(i=end-1;i>=begin;i--){
|
||||
if(mag(evec2[i])>mag(evec2[i+1])) {
|
||||
swap(evec2+i,evec2+i+1);
|
||||
swapped=1;
|
||||
}
|
||||
}
|
||||
begin++;
|
||||
}
|
||||
|
||||
for(i=0;i<row_num;i++){
|
||||
for(j=0;j<row_num;j++) {
|
||||
if(i==j) H[i*row_num+j]=evec2[i+1];
|
||||
else H[i*row_num+j]=0.;
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
static void bisec(std::vector<RealD> &c,
|
||||
std::vector<RealD> &b,
|
||||
int n,
|
||||
int m1,
|
||||
int m2,
|
||||
RealD eps1,
|
||||
RealD relfeh,
|
||||
std::vector<RealD> &x,
|
||||
RealD &eps2)
|
||||
{
|
||||
std::vector<RealD> wu(n+2);
|
||||
|
||||
RealD h,q,x1,xu,x0,xmin,xmax;
|
||||
int i,a,k;
|
||||
|
||||
b[1]=0.0;
|
||||
xmin=c[n]-fabs(b[n]);
|
||||
xmax=c[n]+fabs(b[n]);
|
||||
for(i=1;i<n;i++){
|
||||
h=fabs(b[i])+fabs(b[i+1]);
|
||||
if(c[i]+h>xmax) xmax= c[i]+h;
|
||||
if(c[i]-h<xmin) xmin= c[i]-h;
|
||||
}
|
||||
xmax *=2.;
|
||||
|
||||
eps2=relfeh*((xmin+xmax)>0.0 ? xmax : -xmin);
|
||||
if(eps1<=0.0) eps1=eps2;
|
||||
eps2=0.5*eps1+7.0*(eps2);
|
||||
x0=xmax;
|
||||
for(i=m1;i<=m2;i++){
|
||||
x[i]=xmax;
|
||||
wu[i]=xmin;
|
||||
}
|
||||
|
||||
for(k=m2;k>=m1;k--){
|
||||
xu=xmin;
|
||||
i=k;
|
||||
do{
|
||||
if(xu<wu[i]){
|
||||
xu=wu[i];
|
||||
i=m1-1;
|
||||
}
|
||||
i--;
|
||||
}while(i>=m1);
|
||||
if(x0>x[k]) x0=x[k];
|
||||
while((x0-xu)>2*relfeh*(fabs(xu)+fabs(x0))+eps1){
|
||||
x1=(xu+x0)/2;
|
||||
|
||||
a=0;
|
||||
q=1.0;
|
||||
for(i=1;i<=n;i++){
|
||||
q=c[i]-x1-((q!=0.0)? b[i]*b[i]/q:fabs(b[i])/relfeh);
|
||||
if(q<0) a++;
|
||||
}
|
||||
// printf("x1=%0.14e a=%d\n",x1,a);
|
||||
if(a<k){
|
||||
if(a<m1){
|
||||
xu=x1;
|
||||
wu[m1]=x1;
|
||||
}else {
|
||||
xu=x1;
|
||||
wu[a+1]=x1;
|
||||
if(x[a]>x1) x[a]=x1;
|
||||
}
|
||||
}else x0=x1;
|
||||
}
|
||||
printf("x0=%0.14e xu=%0.14e k=%d\n",x0,xu,k);
|
||||
x[k]=(x0+xu)/2;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Implicitly restarted lanczos
|
||||
/////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
template<class Field>
|
||||
class SimpleLanczos {
|
||||
|
||||
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 OrthoTime;
|
||||
|
||||
RealD eresid;
|
||||
|
||||
SortEigen<Field> _sort;
|
||||
|
||||
LinearOperatorBase<Field> &_Linop;
|
||||
|
||||
OperatorFunction<Field> &_poly;
|
||||
|
||||
/////////////////////////
|
||||
// Constructor
|
||||
/////////////////////////
|
||||
void init(void){};
|
||||
void Abort(int ff, DenseVector<RealD> &evals, DenseVector<DenseVector<RealD> > &evecs);
|
||||
|
||||
SimpleLanczos(
|
||||
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);
|
||||
};
|
||||
|
||||
/////////////////////////
|
||||
// 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<<GridLogMessage << "RitzMatrix "<<std::endl;
|
||||
for(int i=0;i<k;i++){
|
||||
_poly(_Linop,evec[i],w);
|
||||
std::cout<<GridLogMessage << "["<<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<<GridLogMessage<<"oops"<<std::endl;
|
||||
abort();
|
||||
} else
|
||||
std::cout<<GridLogMessage << " 0 ";
|
||||
} else {
|
||||
std::cout<<GridLogMessage << " "<<in<<" ";
|
||||
}
|
||||
}
|
||||
std::cout<<GridLogMessage << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void step(DenseVector<RealD>& lmd,
|
||||
DenseVector<RealD>& lme,
|
||||
Field& last,
|
||||
Field& current,
|
||||
Field & next,
|
||||
uint64_t k)
|
||||
{
|
||||
if (lmd.size()<=k) lmd.resize(k+Nm);
|
||||
if (lme.size()<=k) lme.resize(k+Nm);
|
||||
|
||||
|
||||
_poly(_Linop,current,next ); // 3. wk:=Avk−βkv_{k−1}
|
||||
if(k>0){
|
||||
next -= lme[k-1] * last;
|
||||
}
|
||||
// std::cout<<GridLogMessage << "<last|next>" << innerProduct(last,next) <<std::endl;
|
||||
|
||||
ComplexD zalph = innerProduct(current,next); // 4. αk:=(wk,vk)
|
||||
RealD alph = real(zalph);
|
||||
|
||||
next = next - alph * current ;// 5. wk:=wk−αkvk
|
||||
// std::cout<<GridLogMessage << "<current|next>" << innerProduct(current,next) <<std::endl;
|
||||
|
||||
RealD beta = normalise(next); // 6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop
|
||||
// 7. vk+1 := wk/βk+1
|
||||
// norm=beta;
|
||||
|
||||
int interval = Nm/100+1;
|
||||
if ((k%interval)==0)
|
||||
std::cout<<GridLogMessage << k << " : alpha = " << zalph << " beta "<<beta<<std::endl;
|
||||
const RealD tiny = 1.0e-20;
|
||||
if ( beta < tiny ) {
|
||||
std::cout<<GridLogMessage << " beta is tiny "<<beta<<std::endl;
|
||||
}
|
||||
lmd[k] = alph;
|
||||
lme[k] = beta;
|
||||
|
||||
}
|
||||
|
||||
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
|
||||
#ifdef USE_MKL
|
||||
#define LAPACK_INT MKL_INT
|
||||
#else
|
||||
#define LAPACK_INT long long
|
||||
#endif
|
||||
void diagonalize_lapack(DenseVector<RealD>& lmd,
|
||||
DenseVector<RealD>& lme,
|
||||
int N1, // all
|
||||
int N2, // get
|
||||
GridBase *grid) {
|
||||
const int size = Nm;
|
||||
LAPACK_INT NN = N1;
|
||||
double evals_tmp[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];
|
||||
}
|
||||
LAPACK_INT evals_found;
|
||||
LAPACK_INT lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ;
|
||||
LAPACK_INT liwork = 3+NN*10 ;
|
||||
LAPACK_INT iwork[liwork];
|
||||
double work[lwork];
|
||||
LAPACK_INT isuppz[2*NN];
|
||||
char jobz = 'N'; // calculate evals only
|
||||
char range = 'I'; // calculate il-th to iu-th 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];
|
||||
LAPACK_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;
|
||||
LAPACK_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);
|
||||
#ifdef USE_MKL
|
||||
dstegr(&jobz, &range, &NN,
|
||||
#else
|
||||
LAPACK_dstegr(&jobz, &range, &NN,
|
||||
#endif
|
||||
(double*)DD, (double*)EE,
|
||||
&vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A'
|
||||
&tol, // tolerance
|
||||
&evals_found, evals_tmp, (double*)NULL, &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.;
|
||||
}
|
||||
}
|
||||
{
|
||||
grid->GlobalSumVector(evals_tmp,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.
|
||||
}
|
||||
#undef LAPACK_INT
|
||||
#endif
|
||||
|
||||
|
||||
void diagonalize(DenseVector<RealD>& lmd,
|
||||
DenseVector<RealD>& lme,
|
||||
int N2,
|
||||
int N1,
|
||||
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,grid);
|
||||
|
||||
// diagonalize_lapack(lmd2,lme2,Nm2,Nm,Qt,grid);
|
||||
#endif
|
||||
}
|
||||
|
||||
#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)
|
||||
{
|
||||
double t0=-usecond()/1e6;
|
||||
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);
|
||||
t0+=usecond()/1e6;
|
||||
OrthoTime +=t0;
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
void calc(
|
||||
DenseVector<RealD>& eval,
|
||||
const Field& src,
|
||||
int& Nconv)
|
||||
{
|
||||
|
||||
GridBase *grid = src._grid;
|
||||
// assert(grid == src._grid);
|
||||
|
||||
std::cout<<GridLogMessage << " -- Nk = " << Nk << " Np = "<< Np << std::endl;
|
||||
std::cout<<GridLogMessage << " -- Nm = " << Nm << std::endl;
|
||||
std::cout<<GridLogMessage << " -- size of eval = " << eval.size() << std::endl;
|
||||
|
||||
// assert(c.size() && Nm == eval.size());
|
||||
|
||||
DenseVector<RealD> lme(Nm);
|
||||
DenseVector<RealD> lmd(Nm);
|
||||
|
||||
|
||||
Field current(grid);
|
||||
Field last(grid);
|
||||
Field next(grid);
|
||||
|
||||
Nconv = 0;
|
||||
|
||||
RealD beta_k;
|
||||
|
||||
// Set initial vector
|
||||
// (uniform vector) Why not src??
|
||||
// evec[0] = 1.0;
|
||||
current = src;
|
||||
std:: cout<<GridLogMessage <<"norm2(src)= " << norm2(src)<<std::endl;
|
||||
normalise(current);
|
||||
std:: cout<<GridLogMessage <<"norm2(evec[0])= " << norm2(current) <<std::endl;
|
||||
|
||||
// Initial Nk steps
|
||||
OrthoTime=0.;
|
||||
double t0=usecond()/1e6;
|
||||
RealD norm; // sqrt norm of last vector
|
||||
|
||||
uint64_t iter=0;
|
||||
|
||||
while(1){
|
||||
std::vector <RealD> lme2(Nm);
|
||||
std::vector <RealD> lmd2(Nm);
|
||||
for(uint64_t k=0; k<Nm ; ++k,iter++){
|
||||
step(lmd,lme,last,current,next,iter);
|
||||
last=current;
|
||||
current=next;
|
||||
}
|
||||
double t1=usecond()/1e6;
|
||||
std::cout<<GridLogMessage <<"IRL::Initial steps: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
|
||||
std::cout<<GridLogMessage <<"IRL::Initial steps:OrthoTime "<<OrthoTime<< "seconds"<<std::endl;
|
||||
|
||||
// getting eigenvalues
|
||||
lmd2.resize(iter+2);
|
||||
lme2.resize(iter+2);
|
||||
for(uint64_t k=0; k<iter; ++k){
|
||||
lmd2[k+1] = lmd[k];
|
||||
lme2[k+2] = lme[k];
|
||||
}
|
||||
t1=usecond()/1e6;
|
||||
std::cout<<GridLogMessage <<"IRL:: copy: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
|
||||
{
|
||||
int total = grid->_Nprocessors;
|
||||
int node = grid->_processor;
|
||||
int interval = (Nstop/total)+1;
|
||||
int iu = (iter+1) - (interval*node+1);
|
||||
int il = (iter+1) - (interval*(node+1));
|
||||
RealD eps2;
|
||||
Bisection::bisec(lmd2,lme2,iter,il,iu,1e-16,1e-10, eval,eps2);
|
||||
// diagonalize(eval2,lme2,iter,Nk,grid);
|
||||
for(int i=il;i<=iu;i++)
|
||||
printf("eval[%d]=%0.14e\n",i,eval[i]);
|
||||
t1=usecond()/1e6;
|
||||
std::cout<<GridLogMessage <<"IRL:: diagonalize: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
|
||||
}
|
||||
|
||||
for(uint64_t k=0; k<Nk; ++k){
|
||||
// eval[k] = eval2[k];
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**
|
||||
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 = fabs(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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user