00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00040 #ifndef MAT_LANCZOSSEVERALLARGESTMAGNITUDEEIG
00041 #define MAT_LANCZOSSEVERALLARGESTMAGNITUDEEIG
00042
00043 #include <limits>
00044 #include <vector>
00045
00046 namespace mat {
00047 namespace arn {
00048
00049 template<typename Treal, typename Tmatrix, typename Tvector>
00050 class LanczosSeveralLargestEig
00051 {
00052 public:
00053
00054
00055
00056
00057
00058
00059
00060 LanczosSeveralLargestEig(Tmatrix const & AA, Tvector const & startVec, int num_eigs,
00061 int maxit = 100, int cap = 100, Tvector * deflVec_ = NULL, Treal sigma_ = 0)
00062 : A(AA),
00063 v(new Tvector[cap]),
00064 eigVectorTri(0),
00065 capacity(cap),
00066 j(0),
00067 maxIter(maxit),
00068 eValTmp(0),
00069 accTmp(0),
00070 number_of_eigenv(num_eigs),
00071 alpha(0),
00072 beta(0),
00073 use_selective_orth(false),
00074 use_full_orth(true),
00075 counter_all(0),
00076 counter_orth(0),
00077 deflVec(deflVec_)
00078 {
00079 assert(cap > 1);
00080 Treal const ONE = 1.0;
00081 v[0] = startVec;
00082 if(v[0].eucl() < template_blas_sqrt(getRelPrecision<Treal>())) {
00083 v[0].rand();
00084 }
00085 v[0] *= (ONE / v[0].eucl());
00086 r = v[0];
00087
00088 if(number_of_eigenv == 1)
00089 {unset_use_full_orth(); unset_use_selective_orth();}
00090
00091 absTol = 1e-12;
00092 relTol = 1e-12;
00093 sigma = sigma_;
00094
00095 }
00096
00097
00098
00099
00100 void setRelTol(Treal const newTol) { relTol = newTol; }
00101 void setAbsTol(Treal const newTol) { absTol = newTol; }
00102
00103 void set_use_selective_orth(){ use_selective_orth = true; }
00104 void set_use_full_orth(){ use_full_orth = true; }
00105 void unset_use_selective_orth(){ use_selective_orth = false; }
00106 void unset_use_full_orth(){ use_full_orth = false; }
00107
00108 virtual void run() {
00109 do {
00110 if(j > 1 && use_selective_orth)
00111 selective_orth();
00112 step();
00113 update();
00114 if (j > maxIter)
00115 throw AcceptableMaxIter("Lanczos::run() did not converge within maxIter");
00116 }
00117 while (!converged());
00118 }
00119
00120
00121
00122
00123
00124 virtual void get_ith_eigenpair(int i, Treal& eigVal, Tvector& eigVec, Treal & acc)
00125 {
00126 assert(i > 0);
00127 assert(i <= size_accTmp);
00128 eigVal = eValTmp[size_accTmp - i];
00129 assert(eigVectorTri);
00130 getEigVector(eigVec, &eigVectorTri[j * (size_accTmp - i)]);
00131 acc = accTmp[size_accTmp - i];
00132 }
00133
00134 int get_num_iter() const{ return j;}
00135
00136 virtual ~LanczosSeveralLargestEig() {
00137
00138 if(use_selective_orth)
00139 printf("Orthogonalized %d of total possible %d, this is %lf %%\n", counter_orth, counter_all, (double)counter_orth/counter_all*100);
00140
00141 delete[] eigVectorTri;
00142 delete[] eValTmp;
00143 delete[] accTmp;
00144 delete[] v;
00145 }
00146
00147
00148 inline void copyTridiag(MatrixTridiagSymmetric<Treal> & Tricopy) {
00149 Tricopy = Tri;
00150 }
00151
00152
00153 protected:
00154 Tmatrix const & A;
00155 Tvector* v;
00160 Tvector r;
00161 MatrixTridiagSymmetric<Treal> Tri;
00162 Treal* eigVectorTri;
00163 int capacity;
00164 int j;
00165 int maxIter;
00166 void increaseCapacity(int const newCapacity);
00167 void getEigVector(Tvector& eigVec, Treal const * const eVecTri) const;
00168 Treal absTol;
00169 Treal relTol;
00170 virtual void step();
00171 virtual void computeEigenPairTri();
00172 virtual void update() {
00173 computeEigenPairTri();
00174 }
00175 void selective_orth();
00176 virtual bool converged() const;
00177 virtual bool converged_ith(int i) const;
00178 Treal* eValTmp;
00179 Treal* accTmp;
00180 int number_of_eigenv;
00181 int size_accTmp;
00182 private:
00183 Treal alpha;
00184 Treal beta;
00185
00186 bool use_selective_orth;
00187 bool use_full_orth;
00188
00189 int counter_all;
00190 int counter_orth;
00191
00192
00193 Tvector * deflVec;
00194 Treal sigma;
00195 };
00196
00197 template<typename Treal, typename Tmatrix, typename Tvector>
00198 void LanczosSeveralLargestEig<Treal, Tmatrix, Tvector>::
00199 selective_orth()
00200 {
00201 int j_curr = j-1;
00202
00203 Treal coeff = 0, res;
00204 Treal normT = 0;
00205
00206 for(int i = 0; i <= j_curr; ++i)
00207 if(template_blas_fabs(eValTmp[i]) > normT) normT = template_blas_fabs(eValTmp[i]);
00208
00209 Treal epsilon = mat::getMachineEpsilon<Treal>();
00210 Tvector tmp;
00211 tmp = v[j_curr+1];
00212 tmp *= beta;
00213
00214 for(int i = j_curr; i >= 0; --i)
00215 {
00216 counter_all++;
00217
00218 res = accTmp[i];
00219 Treal tol = template_blas_sqrt(epsilon) * normT;
00220 if(res <= tol)
00221 {
00222 counter_orth++;
00223 Tvector eigVec;
00224 getEigVector(eigVec, &eigVectorTri[j_curr * i]);
00225 coeff = transpose(eigVec) * tmp;
00226 tmp += (-coeff) * (eigVec);
00227 }
00228 }
00229
00230
00231
00232 v[j_curr+1] = tmp;
00233 beta = v[j_curr+1].eucl();
00234 Treal const ONE = 1.0;
00235 v[j_curr+1] *= ONE / beta;
00236 Tri.update_beta(beta);
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 }
00255
00256
00257
00258
00259 template<typename Treal, typename Tmatrix, typename Tvector>
00260 void LanczosSeveralLargestEig<Treal, Tmatrix, Tvector>::
00261 step()
00262 {
00263 if (j + 1 >= capacity)
00264 increaseCapacity(capacity * 2);
00265 Treal const ONE = 1.0;
00266 A.matVecProd(r, v[j]);
00267 alpha = transpose(v[j]) * r;
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 if(deflVec != NULL)
00281 {
00282
00283
00284
00285
00286
00287 Treal gamma = transpose(*deflVec) * v[j];
00288 alpha -= sigma*gamma*gamma;
00289
00290 r += (-sigma*gamma) * (*deflVec);
00291 }
00292
00293 r += (-alpha) * v[j];
00294
00295 if (j) {
00296 r += (-beta) * v[j-1];
00297 v[j-1].writeToFile();
00298 }
00299
00300
00301
00302 if(use_full_orth)
00303 {
00304
00305
00306 Treal gamma_i = 0;
00307 Tvector tmp;
00308 tmp = r;
00309 for(int i = 0; i < j; ++i )
00310 {
00311 v[i].readFromFile();
00312 gamma_i = transpose(tmp) * v[i];
00313 r += (-gamma_i) * v[i];
00314 v[i].writeToFile();
00315 }
00316 tmp.clear();
00317 }
00318
00319
00320 beta = r.eucl();
00321 v[j+1] = r;
00322 v[j+1] *= ONE / beta;
00323 Tri.increase(alpha, beta);
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 j++;
00345 }
00346
00347
00348
00349
00350
00351 template<typename Treal, typename Tmatrix, typename Tvector>
00352 void LanczosSeveralLargestEig<Treal, Tmatrix, Tvector>::
00353 computeEigenPairTri() {
00354 if( eigVectorTri != NULL ) delete[] eigVectorTri;
00355 if( accTmp != NULL ) delete[] accTmp;
00356 if( eValTmp != NULL ) delete[] eValTmp;
00357
00358 int num_compute_eigenvalues;
00359 if(use_selective_orth)
00360 num_compute_eigenvalues = j;
00361 else
00362 num_compute_eigenvalues = number_of_eigenv;
00363
00364
00365 int const max_ind = j-1;
00366 int const min_ind = std::max(j - num_compute_eigenvalues, 0);
00367
00368 Treal* eigVectors = new Treal[j * num_compute_eigenvalues];
00369 Treal* eigVals = new Treal[num_compute_eigenvalues];
00370 Treal* accMax = new Treal[num_compute_eigenvalues];
00371 assert(eigVectors != NULL);
00372 assert(eigVals != NULL);
00373 assert(accMax != NULL);
00374
00375 Tri.getEigsByIndex(eigVals, eigVectors, accMax,
00376 min_ind, max_ind);
00377
00378 eValTmp = eigVals;
00379
00380
00381 eigVectorTri = eigVectors;
00382 accTmp = accMax;
00383 size_accTmp = num_compute_eigenvalues;
00384
00385
00386 eigVectors = NULL;
00387 eigVals = NULL;
00388 accMax = NULL;
00389 }
00390
00391
00392
00393
00394
00395
00396
00397
00398 template<typename Treal, typename Tmatrix, typename Tvector>
00399 void LanczosSeveralLargestEig<Treal, Tmatrix, Tvector>::
00400 getEigVector(Tvector& eigVec, Treal const * const eVecTri) const {
00401 if (j <= 1) {
00402 eigVec = v[0];
00403 }
00404 else {
00405 v[0].readFromFile();
00406 eigVec = v[0];
00407 v[0].writeToFile();
00408 }
00409 eigVec *= eVecTri[0];
00410 for (int ind = 1; ind <= j - 2; ++ind) {
00411 v[ind].readFromFile();
00412 eigVec += eVecTri[ind] * v[ind];
00413 v[ind].writeToFile();
00414 }
00415 eigVec += eVecTri[j-1] * v[j-1];
00416
00417
00418 Treal norm_eigVec = eigVec.eucl();
00419 Treal const ONE = 1.0;
00420 eigVec *= ONE / norm_eigVec;
00421 }
00422
00423
00424
00425 template<typename Treal, typename Tmatrix, typename Tvector>
00426 bool LanczosSeveralLargestEig<Treal, Tmatrix, Tvector>::
00427 converged() const {
00428
00429 if(j < number_of_eigenv) return false;
00430 bool conv = converged_ith(number_of_eigenv);
00431
00432 return conv;
00433 }
00434
00435
00436 template<typename Treal, typename Tmatrix, typename Tvector>
00437 bool LanczosSeveralLargestEig<Treal, Tmatrix, Tvector>::
00438 converged_ith(int i) const {
00439 assert(size_accTmp >= i);
00440
00441 bool conv = true;
00442 if (template_blas_fabs(eValTmp[size_accTmp - i]) > 0) {
00443 conv = conv &&
00444 accTmp[size_accTmp - i] / template_blas_fabs(eValTmp[size_accTmp - i]) < relTol;
00445 }
00446 return conv;
00447 }
00448
00449
00450 template<typename Treal, typename Tmatrix, typename Tvector>
00451 void LanczosSeveralLargestEig<Treal, Tmatrix, Tvector>::
00452 increaseCapacity(int const newCapacity) {
00453 assert(newCapacity > capacity);
00454 assert(j > 0);
00455 capacity = newCapacity;
00456 Tvector* new_v = new Tvector[capacity];
00457 assert(new_v != NULL);
00458
00459
00460
00461 for (int ind = 0; ind <= j - 2; ind++){
00462 v[ind].readFromFile();
00463 new_v[ind] = v[ind];
00464 new_v[ind].writeToFile();
00465 }
00466 for (int ind = j - 1; ind <= j; ind++){
00467 new_v[ind] = v[ind];
00468 }
00469 delete[] v;
00470 v = new_v;
00471 }
00472
00473
00474 }
00475
00476
00477 }
00478 #endif