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());
167 void getEigVector(Tvector& eigVec, Treal
const *
const eVecTri)
const;
197 template<
typename Treal,
typename Tmatrix,
typename Tvector>
203 Treal coeff = 0, res;
206 for(
int i = 0; i <= j_curr; ++i)
209 Treal epsilon = mat::getMachineEpsilon<Treal>();
214 for(
int i = j_curr; i >= 0; --i)
224 getEigVector(eigVec, &eigVectorTri[j_curr * i]);
226 tmp += (-coeff) * (eigVec);
233 beta = v[j_curr+1].eucl();
234 Treal
const ONE = 1.0;
235 v[j_curr+1] *= ONE / beta;
236 Tri.update_beta(beta);
259 template<
typename Treal,
typename Tmatrix,
typename Tvector>
263 if (j + 1 >= capacity)
264 increaseCapacity(capacity * 2);
265 Treal
const ONE = 1.0;
266 A.matVecProd(r, v[j]);
287 Treal gamma =
transpose(*deflVec) * v[j];
288 alpha -= sigma*gamma*gamma;
290 r += (-sigma*gamma) * (*deflVec);
293 r += (-alpha) * v[j];
296 r += (-beta) * v[j-1];
297 v[j-1].writeToFile();
309 for(
int i = 0; i < j; ++i )
313 r += (-gamma_i) * v[i];
322 v[j+1] *= ONE / beta;
323 Tri.increase(alpha, beta);
351 template<
typename Treal,
typename Tmatrix,
typename Tvector>
354 if( eigVectorTri != NULL )
delete[] eigVectorTri;
355 if( accTmp != NULL )
delete[] accTmp;
356 if( eValTmp != NULL )
delete[] eValTmp;
358 int num_compute_eigenvalues;
359 if(use_selective_orth)
360 num_compute_eigenvalues = j;
362 num_compute_eigenvalues = number_of_eigenv;
365 int const max_ind = j-1;
366 int const min_ind =
std::max(j - num_compute_eigenvalues, 0);
368 Treal* eigVectors =
new Treal[j * num_compute_eigenvalues];
369 Treal* eigVals =
new Treal[num_compute_eigenvalues];
370 Treal* accMax =
new Treal[num_compute_eigenvalues];
371 assert(eigVectors != NULL);
372 assert(eigVals != NULL);
373 assert(accMax != NULL);
375 Tri.getEigsByIndex(eigVals, eigVectors, accMax,
381 eigVectorTri = eigVectors;
383 size_accTmp = num_compute_eigenvalues;
398 template<
typename Treal,
typename Tmatrix,
typename Tvector>
409 eigVec *= eVecTri[0];
410 for (
int ind = 1; ind <= j - 2; ++ind) {
411 v[ind].readFromFile();
412 eigVec += eVecTri[ind] * v[ind];
413 v[ind].writeToFile();
415 eigVec += eVecTri[j-1] * v[j-1];
418 Treal norm_eigVec = eigVec.eucl();
419 Treal
const ONE = 1.0;
420 eigVec *= ONE / norm_eigVec;
425 template<
typename Treal,
typename Tmatrix,
typename Tvector>
429 if(j < number_of_eigenv)
return false;
430 bool conv = converged_ith(number_of_eigenv);
436 template<
typename Treal,
typename Tmatrix,
typename Tvector>
439 assert(size_accTmp >= i);
450 template<
typename Treal,
typename Tmatrix,
typename Tvector>
453 assert(newCapacity > capacity);
455 capacity = newCapacity;
456 Tvector* new_v =
new Tvector[capacity];
457 assert(new_v != NULL);
461 for (
int ind = 0; ind <= j - 2; ind++){
462 v[ind].readFromFile();
464 new_v[ind].writeToFile();
466 for (
int ind = j - 1; ind <= j; ind++){
Definition: LanczosSeveralLargestEig.h:50
void setRelTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:100
Treal sigma
Definition: LanczosSeveralLargestEig.h:194
Tvector * deflVec
Definition: LanczosSeveralLargestEig.h:193
void increaseCapacity(int const newCapacity)
Definition: LanczosSeveralLargestEig.h:452
void unset_use_full_orth()
Definition: LanczosSeveralLargestEig.h:106
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
int counter_all
Definition: LanczosSeveralLargestEig.h:189
Tridiagonal symmetric matrix class template.
Definition: MatrixTridiagSymmetric.h:47
virtual bool converged() const
Definition: LanczosSeveralLargestEig.h:427
virtual void step()
Definition: LanczosSeveralLargestEig.h:261
Tvector * v
Definition: LanczosSeveralLargestEig.h:155
virtual void computeEigenPairTri()
Definition: LanczosSeveralLargestEig.h:353
Definition: allocate.cc:39
int j
Definition: LanczosSeveralLargestEig.h:164
Treal template_blas_fabs(Treal x)
void copyTridiag(MatrixTridiagSymmetric< Treal > &Tricopy)
Definition: LanczosSeveralLargestEig.h:148
Tmatrix const & A
Definition: LanczosSeveralLargestEig.h:154
Treal * eigVectorTri
Definition: LanczosSeveralLargestEig.h:162
#define max(a, b)
Definition: integrator.cc:87
void selective_orth()
Definition: LanczosSeveralLargestEig.h:199
Treal * accTmp
Definition: LanczosSeveralLargestEig.h:179
void set_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:103
Treal * eValTmp
Definition: LanczosSeveralLargestEig.h:178
virtual void get_ith_eigenpair(int i, Treal &eigVal, Tvector &eigVec, Treal &acc)
Definition: LanczosSeveralLargestEig.h:124
virtual ~LanczosSeveralLargestEig()
Definition: LanczosSeveralLargestEig.h:136
void unset_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:105
virtual bool converged_ith(int i) const
Definition: LanczosSeveralLargestEig.h:438
bool use_full_orth
Definition: LanczosSeveralLargestEig.h:187
void getEigVector(Tvector &eigVec, Treal const *const eVecTri) const
Definition: LanczosSeveralLargestEig.h:400
Treal absTol
Definition: LanczosSeveralLargestEig.h:168
Treal relTol
Definition: LanczosSeveralLargestEig.h:169
int size_accTmp
Definition: LanczosSeveralLargestEig.h:181
int get_num_iter() const
Definition: LanczosSeveralLargestEig.h:134
void set_use_full_orth()
Definition: LanczosSeveralLargestEig.h:104
bool use_selective_orth
Definition: LanczosSeveralLargestEig.h:186
void setAbsTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:101
MatrixTridiagSymmetric< Treal > Tri
Residual vector.
Definition: LanczosSeveralLargestEig.h:161
Treal alpha
Definition: LanczosSeveralLargestEig.h:183
int maxIter
Current step.
Definition: LanczosSeveralLargestEig.h:165
int number_of_eigenv
Definition: LanczosSeveralLargestEig.h:180
Tvector r
Vectors spanning Krylov subspace.
Definition: LanczosSeveralLargestEig.h:160
virtual void run()
Definition: LanczosSeveralLargestEig.h:108
int counter_orth
Definition: LanczosSeveralLargestEig.h:190
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
Treal template_blas_sqrt(Treal x)
virtual void update()
Definition: LanczosSeveralLargestEig.h:172
int capacity
Definition: LanczosSeveralLargestEig.h:163
Treal beta
Definition: LanczosSeveralLargestEig.h:184