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