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

Broken Lanczos. Going back to an older verion temporarily.

Browse Source
This commit is contained in:
Chulwoo Jung 2017-06-06 12:14:45 -04:00
parent b1b15f0b70
commit cfed2c1ea0
1 changed files with 24 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
View File

@ -281,14 +281,12 @@ public:
// tevals.resize(size); // tevals.resize(size);
// tevecs.resize(size); // tevecs.resize(size);
LAPACK_INT NN = N1; LAPACK_INT NN = N1;
// double evals_tmp[NN]; double evals_tmp[NN];
// double evec_tmp[NN][NN]; double evec_tmp[NN][NN];
std::vector<double> evals_tmp(NN); memset(evec_tmp[0],0,sizeof(double)*NN*NN);
std::vector<double> evec_tmp(NN*NN); // double AA[NN][NN];
memset(evec_tmp.data(),0,sizeof(double)*NN*NN); double DD[NN];
double EE[NN];
std::vector<double> DD(NN);
std::vector<double> EE(NN);
for (int i = 0; i< NN; i++) for (int i = 0; i< NN; i++)
for (int j = i - 1; j <= i + 1; j++) for (int j = i - 1; j <= i + 1; j++)
if ( j < NN && j >= 0 ) { if ( j < NN && j >= 0 ) {
@ -297,15 +295,11 @@ public:
if (j==(i-1)) EE[j] = lme[j]; if (j==(i-1)) EE[j] = lme[j];
} }
LAPACK_INT evals_found; LAPACK_INT evals_found;
// LAPACK_INT lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ; LAPACK_INT lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ;
LAPACK_INT lwork = 1+(18*NN) ;
LAPACK_INT liwork = 3+NN*10 ; LAPACK_INT liwork = 3+NN*10 ;
// LAPACK_INT iwork[liwork]; LAPACK_INT iwork[liwork];
// double work[lwork]; double work[lwork];
// LAPACK_INT isuppz[2*NN]; LAPACK_INT isuppz[2*NN];
std::vector<LAPACK_INT> iwork(liwork);
std::vector<double> work(lwork);
std::vector<LAPACK_INT> isuppz(2*NN);
char jobz = 'V'; // calculate evals & evecs char jobz = 'V'; // calculate evals & evecs
char range = 'I'; // calculate all evals char range = 'I'; // calculate all evals
// char range = 'A'; // calculate all evals // char range = 'A'; // calculate all evals
@ -324,7 +318,7 @@ public:
if (iu > NN) iu=NN; if (iu > NN) iu=NN;
double tol = 0.0; double tol = 0.0;
if (1) { if (1) {
memset(evals_tmp.data(),0,sizeof(double)*NN); memset(evals_tmp,0,sizeof(double)*NN);
if ( il <= NN){ if ( il <= NN){
printf("total=%d node=%d il=%d iu=%d\n",total,node,il,iu); printf("total=%d node=%d il=%d iu=%d\n",total,node,il,iu);
#ifdef USE_MKL #ifdef USE_MKL
@ -332,32 +326,34 @@ public:
#else #else
LAPACK_dstegr(&jobz, &range, &NN, LAPACK_dstegr(&jobz, &range, &NN,
#endif #endif
DD.data(), EE.data(), (double*)DD, (double*)EE,
&vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A' &vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A'
&tol, // tolerance &tol, // tolerance
&evals_found, evals_tmp.data(), evec_tmp.data(), &NN, &evals_found, evals_tmp, (double*)evec_tmp, &NN,
isuppz.data(), isuppz,
work.data(), &lwork, iwork.data(), &liwork, work, &lwork, iwork, &liwork,
&info); &info);
for (int i = iu-1; i>= il-1; i--){ 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)]); 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)]; evals_tmp[i] = evals_tmp[i - (il-1)];
if (il>1) evals_tmp[i-(il-1)]=0.; if (il>1) evals_tmp[i-(il-1)]=0.;
for (int j = 0; j< NN; j++){ for (int j = 0; j< NN; j++){
evec_tmp[i*NN+j] = evec_tmp[(i - (il-1))*NN+j]; evec_tmp[i][j] = evec_tmp[i - (il-1)][j];
if (il>1) evec_tmp[(i-(il-1))*NN+j]=0.; if (il>1) evec_tmp[i-(il-1)][j]=0.;
} }
} }
} }
{ {
grid->GlobalSumVector(evals_tmp.data(),NN); // QMP_sum_double_array(evals_tmp,NN);
grid->GlobalSumVector(evec_tmp.data(),NN*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. // 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 i=0;i<NN;i++){
for(int j=0;j<NN;j++) for(int j=0;j<NN;j++)
Qt[(NN-1-i)*N2+j]=evec_tmp[i*NN+j]; Qt[(NN-1-i)*N2+j]=evec_tmp[i][j];
lmd [NN-1-i]=evals_tmp[i]; lmd [NN-1-i]=evals_tmp[i];
} }
} }
@ -1022,9 +1018,9 @@ until convergence
} }
#else #else
// ConvRotate0( Nk, Nm, Qt, evec, eval,eval2,Iconv,Nconv); // ConvRotate0( Nk, Nm, Qt, evec, eval,eval2,Iconv,Nconv);
ConvRotate( Nk, Nm, Qt, evec, eval,eval2,Iconv,Nconv); // ConvRotate( Nk, Nm, Qt, evec, eval,eval2,Iconv,Nconv);
// ConvCheck only counts Iconv[j]=j. ignore Iconv[] // ConvCheck only counts Iconv[j]=j. ignore Iconv[]
// Rotate0(Nm,Qt,evec,0,Nk,Nk); Rotate0(Nm,Qt,evec,0,Nk,Nk);
FinalCheck( Nk, Nm, eval,evec); FinalCheck( Nk, Nm, eval,evec);
//exit(-1); //exit(-1);
// ConvRotate2( Nk, Nm, Qt, evec, eval,eval2,Iconv,Nconv); // ConvRotate2( Nk, Nm, Qt, evec, eval,eval2,Iconv,Nconv);