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