Lanczos.h
Go to the documentation of this file.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
00036 #ifndef MAT_LANCZOS
00037 #define MAT_LANCZOS
00038 #include "MatrixTridiagSymmetric.h"
00039 #include "mat_utils.h"
00040 namespace mat {
00041 namespace arn {
00042
00056 template<typename Treal, typename Tmatrix, typename Tvector>
00057 class Lanczos {
00058 public:
00059 Lanczos(Tmatrix const & AA, Tvector const & startVec,
00060 int maxIt = 100, int cap = 100)
00061 : A(AA), v(new Tvector[cap]), capacity(cap), j(0), maxIter(maxIt),
00062 alpha(0), beta(0) {
00063 assert(cap > 1);
00064 Treal const ONE = 1.0;
00065 v[0] = startVec;
00066 if(v[0].eucl() < template_blas_sqrt(getRelPrecision<Treal>())) {
00067 v[0].rand();
00068 }
00069 v[0] *= (ONE / v[0].eucl());
00070 r = v[0];
00071 }
00072 void restart(Tvector const & startVec) {
00073 j = 0;
00074 alpha = 0;
00075 beta = 0;
00076 delete[] v;
00077 v = new Tvector[capacity];
00078 Treal const ONE = 1.0;
00079 v[0] = startVec;
00080 v[0] *= (ONE / startVec.eucl());
00081 r = startVec;
00082 Tri.clear();
00083 }
00084
00085 virtual void run() {
00086 do {
00087 step();
00088 update();
00089 if (j > maxIter)
00090 throw AcceptableMaxIter("Lanczos::run() did not converge within maxIter");
00091 }
00092 while (!converged());
00093 }
00094
00095 inline void copyTridiag(MatrixTridiagSymmetric<Treal> & Tricopy) {
00096 Tricopy = Tri;
00097 }
00098 virtual ~Lanczos() {
00099 delete[] v;
00100 }
00101 protected:
00102 Tmatrix const & A;
00103 Tvector* v;
00108 Tvector r;
00109 MatrixTridiagSymmetric<Treal> Tri;
00110 int capacity;
00111 int j;
00112 int maxIter;
00113 void increaseCapacity(int const newCapacity);
00114 void step();
00115 void getEigVector(Tvector& eigVec, Treal const * const eVecTri) const;
00116 virtual void update() = 0;
00117 virtual bool converged() const = 0;
00118 private:
00119 Treal alpha;
00120 Treal beta;
00121 };
00122
00123 template<typename Treal, typename Tmatrix, typename Tvector>
00124 void Lanczos<Treal, Tmatrix, Tvector>::step() {
00125 if (j + 1 >= capacity)
00126 increaseCapacity(capacity * 2);
00127 Treal const ONE = 1.0;
00128 A.matVecProd(r, v[j]);
00129 alpha = transpose(v[j]) * r;
00130 r += (-alpha) * v[j];
00131 if (j) {
00132 r += (-beta) * v[j-1];
00133 v[j-1].writeToFile();
00134 }
00135 beta = r.eucl();
00136 v[j+1] = r;
00137 v[j+1] *= ONE / beta;
00138 Tri.increase(alpha, beta);
00139 ++j;
00140 }
00141
00142
00143
00144
00145
00146 template<typename Treal, typename Tmatrix, typename Tvector>
00147 void Lanczos<Treal, Tmatrix, Tvector>::
00148 getEigVector(Tvector& eigVec, Treal const * const eVecTri) const {
00149 if (j <= 1) {
00150 eigVec = v[0];
00151 }
00152 else {
00153 v[0].readFromFile();
00154 eigVec = v[0];
00155 v[0].writeToFile();
00156 }
00157 eigVec *= eVecTri[0];
00158 for (int ind = 1; ind <= j - 2; ++ind) {
00159 v[ind].readFromFile();
00160 eigVec += eVecTri[ind] * v[ind];
00161 v[ind].writeToFile();
00162 }
00163 eigVec += eVecTri[j-1] * v[j-1];
00164 }
00165
00166 template<typename Treal, typename Tmatrix, typename Tvector>
00167 void Lanczos<Treal, Tmatrix, Tvector>::
00168 increaseCapacity(int const newCapacity) {
00169 assert(newCapacity > capacity);
00170 assert(j > 0);
00171 capacity = newCapacity;
00172 Tvector* new_v = new Tvector[capacity];
00173
00174
00175
00176 for (int ind = 0; ind <= j - 2; ind++){
00177 v[ind].readFromFile();
00178 new_v[ind] = v[ind];
00179 new_v[ind].writeToFile();
00180 }
00181 for (int ind = j - 1; ind <= j; ind++){
00182 new_v[ind] = v[ind];
00183 }
00184 delete[] v;
00185 v = new_v;
00186 }
00187
00188
00189 }
00190 }
00191 #endif