Go to the documentation of this file.
40 #ifndef MAT_LANCZOSSEVERALLARGESTMAGNITUDEEIG
41 #define MAT_LANCZOSSEVERALLARGESTMAGNITUDEEIG
49 template<
typename Treal,
typename Tmatrix,
typename Tvector>
61 int maxit = 100,
int cap = 100, Tvector * deflVec_ = NULL, Treal sigma_ = 0)
80 Treal
const ONE = 1.0;
85 v[0] *= (ONE /
v[0].eucl());
128 v[i].readFromFile();
v[k].readFromFile();
131 throw std::runtime_error(
"Lanczos::run() : detected loss of orthogonality! Discard results.");
132 v[i].writeToFile();
v[k].writeToFile();
139 throw std::runtime_error(
"Lanczos::run() : detected loss of orthogonality! Discard results.");
192 void getEigVector(Tvector& eigVec, Treal
const *
const eVecTri)
const;
224 template<
typename Treal,
typename Tmatrix,
typename Tvector>
230 Treal coeff = 0, res;
233 for(
int i = 0; i <= j_curr; ++i)
236 Treal epsilon = mat::getMachineEpsilon<Treal>();
241 for(
int i = j_curr; i >= 0; --i)
251 getEigVector(eigVec, &eigVectorTri[j_curr * i]);
253 tmp += (-coeff) * (eigVec);
260 beta = v[j_curr+1].eucl();
261 Treal
const ONE = 1.0;
262 v[j_curr+1] *= ONE / beta;
263 Tri.update_beta(beta);
269 template<
typename Treal,
typename Tmatrix,
typename Tvector>
273 if (j + 1 >= capacity)
274 increaseCapacity(capacity * 2);
275 Treal
const ONE = 1.0;
276 A.matVecProd(r, v[j]);
297 Treal gamma =
transpose(*deflVec) * v[j];
298 alpha -= sigma*gamma*gamma;
300 r += (-sigma*gamma) * (*deflVec);
303 r += (-alpha) * v[j];
306 r += (-beta) * v[j-1];
307 v[j-1].writeToFile();
318 for(
int i = 0; i < j; ++i )
322 r += (-gamma_i) * v[i];
326 r += (-gamma_i) * v[j];
332 v[j+1] *= ONE / beta;
333 Tri.increase(alpha, beta);
341 template<
typename Treal,
typename Tmatrix,
typename Tvector>
344 if( eigVectorTri != NULL )
delete[] eigVectorTri;
345 if( accTmp != NULL )
delete[] accTmp;
346 if( eValTmp != NULL )
delete[] eValTmp;
348 int num_compute_eigenvalues;
349 if(use_selective_orth)
350 num_compute_eigenvalues = j;
352 num_compute_eigenvalues = number_of_eigenv;
355 int const max_ind = j-1;
356 int const min_ind =
std::max(j - num_compute_eigenvalues, 0);
358 Treal* eigVectors =
new Treal[j * num_compute_eigenvalues];
359 Treal* eigVals =
new Treal[num_compute_eigenvalues];
360 Treal* accMax =
new Treal[num_compute_eigenvalues];
361 assert(eigVectors != NULL);
362 assert(eigVals != NULL);
363 assert(accMax != NULL);
365 Tri.getEigsByIndex(eigVals, eigVectors, accMax,
371 eigVectorTri = eigVectors;
373 size_accTmp = num_compute_eigenvalues;
388 template<
typename Treal,
typename Tmatrix,
typename Tvector>
390 getEigVector(Tvector& eigVec, Treal
const *
const eVecTri)
const {
399 eigVec *= eVecTri[0];
400 for (
int ind = 1; ind <= j - 2; ++ind) {
401 v[ind].readFromFile();
402 eigVec += eVecTri[ind] * v[ind];
403 v[ind].writeToFile();
405 eigVec += eVecTri[j-1] * v[j-1];
408 Treal norm_eigVec = eigVec.eucl();
409 Treal
const ONE = 1.0;
410 eigVec *= ONE / norm_eigVec;
415 template<
typename Treal,
typename Tmatrix,
typename Tvector>
419 if(j < number_of_eigenv)
return false;
421 if(number_of_eigenv > 1)
422 conv1 = converged_ith(number_of_eigenv-1);
423 bool conv = converged_ith(number_of_eigenv);
425 return conv && conv1;
429 template<
typename Treal,
typename Tmatrix,
typename Tvector>
432 assert(size_accTmp >= i);
443 template<
typename Treal,
typename Tmatrix,
typename Tvector>
446 assert(newCapacity > capacity);
448 capacity = newCapacity;
449 Tvector* new_v =
new Tvector[capacity];
450 assert(new_v != NULL);
454 for (
int ind = 0; ind <= j - 2; ind++){
455 v[ind].readFromFile();
457 new_v[ind].writeToFile();
459 for (
int ind = j - 1; ind <= j; ind++){
int counter_all
Definition: LanczosSeveralLargestEig.h:216
LanczosSeveralLargestEig(Tmatrix const &AA, Tvector const &startVec, int num_eigs, int maxit=100, int cap=100, Tvector *deflVec_=NULL, Treal sigma_=0)
Definition: LanczosSeveralLargestEig.h:60
Treal template_blas_sqrt(Treal x)
int j
Definition: LanczosSeveralLargestEig.h:189
void copyTridiag(MatrixTridiagSymmetric< Treal > &Tricopy)
Definition: LanczosSeveralLargestEig.h:173
Tvector * deflVec
Definition: LanczosSeveralLargestEig.h:220
void setRelTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:100
int total_num_iter
Definition: LanczosSeveralLargestEig.h:211
Tmatrix const & A
Definition: LanczosSeveralLargestEig.h:179
void set_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:103
bool use_full_orth
Definition: LanczosSeveralLargestEig.h:214
virtual void computeEigenPairTri()
Definition: LanczosSeveralLargestEig.h:343
virtual bool converged() const
Definition: LanczosSeveralLargestEig.h:417
Treal absTol
Definition: LanczosSeveralLargestEig.h:193
Treal beta
Definition: LanczosSeveralLargestEig.h:209
Treal template_blas_fabs(Treal x)
Treal sigma
Definition: LanczosSeveralLargestEig.h:221
int number_of_eigenv
Definition: LanczosSeveralLargestEig.h:205
int maxIter
Current step.
Definition: LanczosSeveralLargestEig.h:190
int capacity
Definition: LanczosSeveralLargestEig.h:188
Tridiagonal symmetric matrix class template.
Definition: MatrixTridiagSymmetric.h:47
Definition: LanczosSeveralLargestEig.h:51
MatrixTridiagSymmetric< Treal > Tri
Residual vector.
Definition: LanczosSeveralLargestEig.h:186
void increaseCapacity(int const newCapacity)
Definition: LanczosSeveralLargestEig.h:445
void unset_use_full_orth()
Definition: LanczosSeveralLargestEig.h:106
void unset_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:105
virtual bool converged_ith(int i) const
Definition: LanczosSeveralLargestEig.h:431
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
Treal * eigVectorTri
Definition: LanczosSeveralLargestEig.h:187
Definition: allocate.cc:39
Treal relTol
Definition: LanczosSeveralLargestEig.h:194
int size_accTmp
Definition: LanczosSeveralLargestEig.h:206
void setAbsTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:101
int counter_orth
Definition: LanczosSeveralLargestEig.h:217
int get_num_iter() const
Definition: LanczosSeveralLargestEig.h:159
#define max(a, b)
Definition: integrator.cc:87
virtual void run()
Definition: LanczosSeveralLargestEig.h:108
virtual ~LanczosSeveralLargestEig()
Definition: LanczosSeveralLargestEig.h:161
bool use_selective_orth
Definition: LanczosSeveralLargestEig.h:213
virtual void step()
Definition: LanczosSeveralLargestEig.h:271
Treal * accTmp
Definition: LanczosSeveralLargestEig.h:204
virtual void update()
Definition: LanczosSeveralLargestEig.h:197
Tvector r
Vectors spanning Krylov subspace.
Definition: LanczosSeveralLargestEig.h:185
void selective_orth()
Definition: LanczosSeveralLargestEig.h:226
void getEigVector(Tvector &eigVec, Treal const *const eVecTri) const
Definition: LanczosSeveralLargestEig.h:390
Treal alpha
Definition: LanczosSeveralLargestEig.h:208
Treal * eValTmp
Definition: LanczosSeveralLargestEig.h:203
Tvector * v
Definition: LanczosSeveralLargestEig.h:180
virtual void get_ith_eigenpair(int i, Treal &eigVal, Tvector &eigVec, Treal &acc)
Definition: LanczosSeveralLargestEig.h:149
void set_use_full_orth()
Definition: LanczosSeveralLargestEig.h:104