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
00037 #ifndef MAT_LANCZOSLARGESTMAGNITUDEEIG
00038 #define MAT_LANCZOSLARGESTMAGNITUDEEIG
00039 #include <limits>
00040 #include "Lanczos.h"
00041 namespace mat {
00042 namespace arn {
00043 template<typename Treal, typename Tmatrix, typename Tvector>
00044 class LanczosLargestMagnitudeEig
00045 : public Lanczos<Treal, Tmatrix, Tvector> {
00046 public:
00047 LanczosLargestMagnitudeEig(Tmatrix const & AA, Tvector const & startVec,
00048 int maxIter = 100, int cap = 100)
00049 : Lanczos<Treal, Tmatrix, Tvector>(AA, startVec, maxIter, cap),
00050 eVal(0), acc(0), eigVectorTri(0), absTol( std::numeric_limits<Treal>::max() ),
00051 relTol(template_blas_sqrt(template_blas_sqrt(getRelPrecision<Treal>()))),
00052 eValTmp(0), accTmp(0) {}
00053 void setRelTol(Treal const newTol) { relTol = newTol; }
00054 void setAbsTol(Treal const newTol) { absTol = newTol; }
00055
00056 inline void getLargestMagnitudeEig(Treal& ev, Treal& accuracy) {
00057 ev = eVal;
00058 accuracy = acc;
00059 }
00060 void getLargestMagnitudeEigPair(Treal& eValue,
00061 Tvector& eVector,
00062 Treal& accuracy);
00063
00064 virtual void run() {
00065 Lanczos<Treal, Tmatrix, Tvector>::run();
00066 computeEigVec();
00067 eVal = eValTmp;
00068 acc = accTmp;
00069 rerun();
00070
00071
00072
00073
00074
00075 if(this->A.get_nrows() == 5) {
00076 rerun();
00077 rerun();
00078 rerun();
00079 }
00080 }
00081
00082 void rerun() {
00083 #if 1
00084
00085
00086
00087
00088 Tvector newResidual(eVec);
00089 newResidual.rand();
00090 Treal sP = transpose(eVec) * newResidual;
00091 newResidual += -sP * eVec;
00092 this->restart(newResidual);
00093 Lanczos<Treal, Tmatrix, Tvector>::run();
00094
00095
00096
00097 if (template_blas_fabs(eValTmp) > template_blas_fabs(eVal)) {
00098 eVal = eValTmp;
00099 acc = accTmp;
00100 computeEigVec();
00101 }
00102
00103 acc = acc / template_blas_fabs(eVal) > 2 * std::numeric_limits<Treal>::epsilon() ? acc : 2 * template_blas_fabs(eVal) * std::numeric_limits<Treal>::epsilon();
00104 #endif
00105 }
00106
00107 virtual ~LanczosLargestMagnitudeEig() {
00108 delete[] eigVectorTri;
00109 }
00110 protected:
00111 Treal eVal;
00112 Tvector eVec;
00113 Treal acc;
00114 Treal* eigVectorTri;
00117 Treal absTol;
00118 Treal relTol;
00119 void computeEigenPairTri();
00120 void computeEigVec();
00121 virtual void update() {
00122 computeEigenPairTri();
00123 }
00124 virtual bool converged() const;
00125 Treal eValTmp;
00126 Treal accTmp;
00127 private:
00128 };
00129
00130 template<typename Treal, typename Tmatrix, typename Tvector>
00131 void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::
00132 getLargestMagnitudeEigPair(Treal& eValue,
00133 Tvector& eVector,
00134 Treal& accuracy) {
00135 eValue = eVal;
00136 accuracy = acc;
00137 eVector = eVec;
00138 }
00139
00140
00141 template<typename Treal, typename Tmatrix, typename Tvector>
00142 void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::
00143 computeEigenPairTri() {
00144 delete[] eigVectorTri;
00145 Treal* eigVectorMax = new Treal[this->j];
00146 Treal* eigVectorMin = new Treal[this->j];
00147
00148
00149 int const lMax = this->j - 1;
00150 Treal eValMax;
00151 Treal accMax;
00152 this->Tri.getEigsByIndex(&eValMax, eigVectorMax, &accMax,
00153 lMax, lMax);
00154
00155 int const lMin = 0;
00156 Treal eValMin;
00157 Treal accMin;
00158 this->Tri.getEigsByIndex(&eValMin, eigVectorMin, &accMin,
00159 lMin, lMin);
00160 if (template_blas_fabs(eValMin) > template_blas_fabs(eValMax)) {
00161 eValTmp = eValMin;
00162 accTmp = accMin;
00163 delete[] eigVectorMax;
00164 eigVectorTri = eigVectorMin;
00165 }
00166 else {
00167 eValTmp = eValMax;
00168 accTmp = accMax;
00169 delete[] eigVectorMin;
00170 eigVectorTri = eigVectorMax;
00171 }
00172 }
00173
00174
00175 template<typename Treal, typename Tmatrix, typename Tvector>
00176 void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::
00177 computeEigVec() {
00178
00179 assert(eigVectorTri);
00180 this->getEigVector(eVec, eigVectorTri);
00181 }
00182
00183
00184 template<typename Treal, typename Tmatrix, typename Tvector>
00185 bool LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::
00186 converged() const {
00187 bool conv = accTmp < absTol;
00188 if (template_blas_fabs(eValTmp) > 0) {
00189 conv = conv &&
00190 accTmp / template_blas_fabs(eValTmp) < relTol;
00191 }
00192 return conv;
00193 }
00194
00195
00196
00197
00198 #if 1
00199 template<typename Treal, typename Tmatrix, typename Tvector>
00200 class LanczosLargestMagnitudeEigIfSmall
00201 : public LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector> {
00202 public:
00203 LanczosLargestMagnitudeEigIfSmall
00204 (Tmatrix const & AA, Tvector const & startVec,
00205 Treal const maxAbsVal,
00206 int maxIter = 100, int cap = 100)
00207 : LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>
00208 (AA, startVec, maxIter, cap), maxAbsValue(maxAbsVal) {
00209 }
00210 bool largestMagEigIsSmall() {return eigIsSmall;}
00211
00212 virtual void run() {
00213 Lanczos<Treal, Tmatrix, Tvector>::run();
00214 this->computeEigVec();
00215 this->eVal = this->eValTmp;
00216 this->acc = this->accTmp;
00217 if (eigIsSmall)
00218 this->rerun();
00219 }
00220
00221 protected:
00222 Treal const maxAbsValue;
00223 bool eigIsSmall;
00224 virtual void update() {
00225 LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::update();
00226 eigIsSmall = template_blas_fabs(this->eValTmp) < maxAbsValue;
00227 }
00228 virtual bool converged() const;
00229 private:
00230 };
00231
00232 template<typename Treal, typename Tmatrix, typename Tvector>
00233 bool LanczosLargestMagnitudeEigIfSmall<Treal, Tmatrix, Tvector>::
00234 converged() const {
00235 bool convAccuracy =
00236 LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::converged();
00237 return convAccuracy || (!eigIsSmall);
00238 }
00239
00240 #endif
00241
00242 }
00243
00245
00246
00247
00248
00249
00256 template<typename Tmatrix, typename Treal>
00257 Interval<Treal> euclIfSmall(Tmatrix const & M,
00258 Treal const requestedAbsAccuracy,
00259 Treal const requestedRelAccuracy,
00260 Treal const maxAbsVal,
00261 typename Tmatrix::VectorType * const eVecPtr = 0) {
00262 assert(requestedAbsAccuracy >= 0);
00263
00264
00265 Treal euclLowerBound;
00266 Treal euclUpperBound;
00267 M.quickEuclBounds(euclLowerBound, euclUpperBound);
00268 if ( euclUpperBound < requestedAbsAccuracy )
00269
00270 return Interval<Treal>( euclLowerBound, euclUpperBound );
00271 if ( euclLowerBound > maxAbsVal )
00272
00273 return Interval<Treal>( euclLowerBound, euclUpperBound );
00274 int maxIter = M.get_nrows() * 100;
00275 typename Tmatrix::VectorType guess;
00276 SizesAndBlocks cols;
00277 M.getCols(cols);
00278 guess.resetSizesAndBlocks(cols);
00279 guess.rand();
00280 mat::arn::LanczosLargestMagnitudeEigIfSmall<Treal, Tmatrix, typename Tmatrix::VectorType>
00281 lan(M, guess, maxAbsVal, maxIter);
00282 lan.setAbsTol( requestedAbsAccuracy );
00283 lan.setRelTol( requestedRelAccuracy );
00284 lan.run();
00285 Treal eVal = 0;
00286 Treal acc = 0;
00287 Treal eValMin = 0;
00288 if (lan.largestMagEigIsSmall()) {
00289 if (eVecPtr)
00290 lan.getLargestMagnitudeEigPair(eVal, *eVecPtr, acc);
00291 else
00292 lan.getLargestMagnitudeEig(eVal, acc);
00293 eVal = template_blas_fabs(eVal);
00294 eValMin = eVal - acc;
00295 eValMin = eValMin >= 0 ? eValMin : (Treal)0;
00296 return Interval<Treal>(eValMin, eVal + acc);
00297 }
00298 else {
00299 eValMin = euclLowerBound;
00300 eValMin = eValMin >= maxAbsVal ? eValMin : maxAbsVal;
00301 return Interval<Treal>(eValMin, euclUpperBound);
00302 }
00303 }
00304
00305 }
00306 #endif