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