diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h index 8cafeb94..03f2bf35 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h @@ -281,14 +281,12 @@ public: // tevals.resize(size); // tevecs.resize(size); LAPACK_INT NN = N1; -// double evals_tmp[NN]; -// double evec_tmp[NN][NN]; - std::vector evals_tmp(NN); - std::vector evec_tmp(NN*NN); - memset(evec_tmp.data(),0,sizeof(double)*NN*NN); - - std::vector DD(NN); - std::vector EE(NN); + 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 ) { @@ -297,15 +295,11 @@ public: 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 lwork = 1+(18*NN) ; + 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]; - std::vector iwork(liwork); - std::vector work(lwork); - std::vector isuppz(2*NN); + LAPACK_INT iwork[liwork]; + double work[lwork]; + LAPACK_INT isuppz[2*NN]; char jobz = 'V'; // calculate evals & evecs char range = 'I'; // calculate all evals // char range = 'A'; // calculate all evals @@ -324,7 +318,7 @@ public: if (iu > NN) iu=NN; double tol = 0.0; if (1) { - memset(evals_tmp.data(),0,sizeof(double)*NN); + 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 @@ -332,32 +326,34 @@ public: #else LAPACK_dstegr(&jobz, &range, &NN, #endif - DD.data(), EE.data(), + (double*)DD, (double*)EE, &vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A' &tol, // tolerance - &evals_found, evals_tmp.data(), evec_tmp.data(), &NN, - isuppz.data(), - work.data(), &lwork, iwork.data(), &liwork, + &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*NN+j] = evec_tmp[(i - (il-1))*NN+j]; - if (il>1) evec_tmp[(i-(il-1))*NN+j]=0.; + evec_tmp[i][j] = evec_tmp[i - (il-1)][j]; + if (il>1) evec_tmp[i-(il-1)][j]=0.; } } } { - grid->GlobalSumVector(evals_tmp.data(),NN); - grid->GlobalSumVector(evec_tmp.data(),NN*NN); +// 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