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_PURISTEPINFO
00037 #define MAT_PURISTEPINFO
00038 #include <math.h>
00039 #include <iomanip>
00040 #include "PuriStepInfoDebug.h"
00041 #define __ID__ "$Id: PuriStepInfo.h 4437 2012-07-05 09:01:18Z elias $"
00042 namespace mat {
00043
00044
00045
00046
00054 template<typename Treal, typename Tvector, typename TdebugPolicy>
00055 class PuriStepInfo : public PuriStepInfoDebug<Treal, TdebugPolicy> {
00056 public:
00057 explicit PuriStepInfo(int nn = -1, int noc = -1, Treal eigvalConvCrit = 0.0)
00058 : n(nn),nocc(noc), traceX(-1.0), traceX2(-1.0), poly(-1),
00059 chosenThresh(0.0), actualThresh(0.0),
00060 estimatedStepsLeft(-1),
00061 eigInterval(-1.0,2.0), delta(-1.0), correctOccupation(0),
00062 XmX2EuclNorm(0.0, 10.0), eigVecPtr(0),
00063 lumoWasComputed(0), homoWasComputed(0),
00064 n0(0.0), n1(0.0),
00065 homo(0.0, 1.0), lumo(0.0, 1.0), eigConvCrit(eigvalConvCrit),
00066 nnzX(0), nnzX2(0), eigAccLoss(0),
00067 timeThresh(0), timeSquare(0), timeXmX2Norm(0), timeTotal(0),
00068 timeXX2Write(0), timeXX2Read(0) {}
00069
00070 ~PuriStepInfo() {
00071 delete eigVecPtr;
00072 }
00073 bool converged() const {
00074 bool homoLumo_converged = (1 - homo.low() < eigConvCrit &&
00075 lumo.upp() < eigConvCrit);
00076 bool XmX2norm_converged = getXmX2EuclNorm().upp() < eigConvCrit;
00077
00078
00079
00080 return homoLumo_converged || XmX2norm_converged;
00081 }
00082
00083 inline void setChosenThresh(Treal const thr) {chosenThresh = thr;}
00084 inline Treal getChosenThresh() const { return chosenThresh;}
00085 inline void setActualThresh(Treal const thr) {actualThresh = thr;}
00086 inline Treal getActualThresh() const { return actualThresh;}
00087 inline void setEstimatedStepsLeft(int const stepsleft) {
00088 estimatedStepsLeft = stepsleft;
00089 }
00090 inline int getEstimatedStepsLeft() const { return estimatedStepsLeft;}
00091 inline void setTraceX(Treal const trX) {traceX = trX;}
00092 inline void setTraceX2(Treal const trX2) {traceX2 = trX2;}
00093 inline Treal getTraceX() const {return traceX;}
00094 inline Treal getTraceX2() const {return traceX2;}
00095 void setPoly();
00096 inline int getPoly() const {return poly;}
00098 inline void setXmX2EuclNorm(Interval<Treal> const XmX2eucl) {
00099 XmX2EuclNorm = XmX2eucl;
00100 }
00101
00102 inline Interval<Treal> getXmX2EuclNorm() const {return XmX2EuclNorm;}
00107 void setEigVecPtr(Tvector * eigVecPtr_) {
00108 delete eigVecPtr;
00109 eigVecPtr = eigVecPtr_;
00110 }
00111 Tvector const * const getEigVecPtr() const {return eigVecPtr;}
00112
00113 void improveHomoLumo(Interval<Treal> const homoInt,
00114 Interval<Treal> const lumoInt);
00115 inline Interval<Treal> const & getHomo() const {
00116 return homo;
00117 }
00118 inline Interval<Treal> const & getLumo() const {
00119 return lumo;
00120 }
00121 inline Interval<Treal> const & getEigInterval() const {
00122 return eigInterval;
00123 }
00124 void exchangeInfoWithNext(PuriStepInfo<Treal, Tvector, TdebugPolicy> & next);
00127 inline void setCorrectOccupation() {correctOccupation = 1;}
00128 inline int getCorrectOccupation() const {return correctOccupation;}
00129 Treal subspaceError() const;
00133 void improveEigInterval(Interval<Treal> const eInt);
00134
00135 inline void setNnzX(size_t const nzX) {nnzX = nzX;}
00136 inline size_t getNnzX() const {return nnzX;}
00137 inline void setNnzX2(size_t const nzX2) {nnzX2 = nzX2;}
00138 inline size_t getNnzX2() const {return nnzX2;}
00139 void computeEigAccLoss();
00140 inline Treal getEigAccLoss() const {return eigAccLoss;}
00141
00142 inline Treal getN0() const {return n0;}
00143 inline Treal getN1() const {return n1;}
00144 inline Treal getDelta() const {return delta;}
00145
00146
00147 inline void checkIntervals(const char* descriptionString) const {
00148 PuriStepInfoDebug<Treal, TdebugPolicy>::
00149 checkIntervals(eigInterval, homo, lumo, XmX2EuclNorm, descriptionString);
00150 }
00151 template<typename Tmatrix>
00152 inline void computeExactValues(Tmatrix const & X, Tmatrix const & X2) {
00153 PuriStepInfoDebug<Treal, TdebugPolicy>::
00154 computeExactValues(X, X2, n, nocc);
00155 }
00156
00157
00158 void setMemUsageBeforeTrunc() {
00159 MemUsage::getMemUsage(memUsageBeforeTrunc);
00160 }
00161 void setMemUsageInXmX2Diff(MemUsage::Values & memUsage) {
00162 memUsageInXmX2Diff = memUsage;
00163 }
00164 void setTimeThresh(float t) {timeThresh = t;}
00165 void setTimeSquare(float t) {timeSquare = t;}
00166 void setTimeXmX2Norm(float t) {timeXmX2Norm = t;}
00167 void setTimeTotal(float t) {timeTotal = t;}
00168 void setTimeXX2Write(float t) {timeXX2Write = t;}
00169 void setTimeXX2Read(float t) {timeXX2Read = t;}
00170
00171 MemUsage::Values getMemUsageBeforeTrunc() { return memUsageBeforeTrunc; }
00172 MemUsage::Values getMemUsageInXmX2Diff() { return memUsageInXmX2Diff; }
00173 float getTimeThresh() {return timeThresh;}
00174 float getTimeSquare() {return timeSquare;}
00175 float getTimeXmX2Norm() {return timeXmX2Norm;}
00176 float getTimeTotal() {return timeTotal;}
00177 float getTimeXX2Write() {return timeXX2Write;}
00178 float getTimeXX2Read() {return timeXX2Read;}
00179
00180 inline bool homoIsAccuratelyKnown
00181 (Treal accuracyLimit
00184 ) const {
00185 return homo.length() < accuracyLimit;
00186 }
00187 inline bool lumoIsAccuratelyKnown
00188 (Treal accuracyLimit
00191 ) const {
00192 return lumo.length() < accuracyLimit;
00193 }
00194
00195 inline bool getLumoWasComputed() {return lumoWasComputed;}
00196 inline bool getHomoWasComputed() {return homoWasComputed;}
00197
00198 protected:
00202 void computen0n1();
00203
00204 int n;
00205 int nocc;
00206
00207 Treal traceX;
00208 Treal traceX2;
00209 int poly;
00210 Treal chosenThresh;
00211 Treal actualThresh;
00212 int estimatedStepsLeft;
00216 Interval<Treal> eigInterval;
00218 Treal delta;
00219 int correctOccupation;
00226 Interval<Treal> XmX2EuclNorm;
00230 Tvector * eigVecPtr;
00233 bool lumoWasComputed;
00234 bool homoWasComputed;
00235 Treal n0;
00238 Treal n1;
00242 Interval<Treal> homo;
00243 Interval<Treal> lumo;
00244 Treal eigConvCrit;
00247 size_t nnzX;
00248 size_t nnzX2;
00254 Treal eigAccLoss;
00255
00256
00257 MemUsage::Values memUsageBeforeTrunc;
00258 MemUsage::Values memUsageInXmX2Diff;
00259 float timeThresh;
00260 float timeSquare;
00261 float timeXmX2Norm;
00262 float timeTotal;
00263 float timeXX2Write;
00264 float timeXX2Read;
00265 private:
00266 };
00267
00268 #if 1
00269 template<typename Treal, typename Tvector, typename TdebugPolicy>
00270 std::ostream& operator<<(std::ostream& s,
00271 PuriStepInfo<Treal, Tvector, TdebugPolicy> const & psi) {
00272 s<<" Trace X: "<<psi.getTraceX()
00273 <<" Trace X2: "<<psi.getTraceX2()
00274 <<" poly: "<<psi.getPoly()
00275 <<" ||E||_2: "<<psi.getActualThresh()
00276 <<" Eiginterval: "<<psi.getEigInterval()
00277 <<" correctOccupation: "<<psi.getCorrectOccupation()
00278 <<std::endl
00279 <<" XmX2EuclNorm: "<<psi.getXmX2EuclNorm()
00280 <<" n0: "<<psi.getN0()
00281 <<" n1: "<<psi.getN1()
00282 <<" homo: "<<psi.getHomo()
00283 <<" lumo: "<<psi.getLumo()
00284 <<" nnzX: "<<psi.getNnzX()
00285 <<" nnzX2: "<<psi.getNnzX2();
00286 return s;
00287 }
00288
00289 #endif
00290
00291 template<typename Treal, typename Tvector, typename TdebugPolicy>
00292 void PuriStepInfo<Treal, Tvector, TdebugPolicy>::setPoly() {
00293 if (Interval<Treal>::intersect(homo,lumo).empty()) {
00294 ASSERTALWAYS(homo.low() > lumo.upp());
00295 if (homo.low() + lumo.upp() > 1)
00296 poly = 0;
00297 else
00298 poly = 1;
00299 }
00300 else {
00301 if (template_blas_fabs(traceX2 - nocc) < template_blas_fabs(2 * traceX - traceX2 - nocc))
00302 poly = 0;
00303 else
00304 poly = 1;
00305 }
00306 }
00307
00308 template<typename Treal, typename Tvector, typename TdebugPolicy>
00309 void PuriStepInfo<Treal, Tvector, TdebugPolicy>::
00310 improveHomoLumo(Interval<Treal> const homoInt,
00311 Interval<Treal> const lumoInt) {
00312 checkIntervals("PuriStepInfo::improveHomoLumo A0");
00313 homo.intersect(homoInt);
00314 checkIntervals("PuriStepInfo::improveHomoLumo A1");
00315 lumo.intersect(lumoInt);
00316 checkIntervals("PuriStepInfo::improveHomoLumo A2");
00317 ASSERTALWAYS(!homo.empty());
00318 ASSERTALWAYS(!lumo.empty());
00319 if (homo.low() > 0.5 && lumo.upp() < 0.5)
00320 this->setCorrectOccupation();
00321
00322 if (correctOccupation && 1.0 - XmX2EuclNorm.upp() * 4.0 > 0) {
00323 Interval<Treal> tmp =
00324 sqrtInt((XmX2EuclNorm * (Treal)(-4.0)) + (Treal)1.0);
00325 ASSERTALWAYS(tmp.length() > 0);
00326
00327 ASSERTALWAYS(!homo.empty());
00328 ASSERTALWAYS(!lumo.empty());
00329
00330 homo.intersect(Interval<Treal>((1.0 + tmp.low()) / 2.0, homo.upp()));
00331
00332 ASSERTALWAYS(!homo.empty());
00333 ASSERTALWAYS(!lumo.empty());
00334
00335 lumo.intersect(Interval<Treal>(lumo.low(), (1.0 - tmp.low()) / 2.0));
00336 checkIntervals("PuriStepInfo::improveHomoLumo B");
00337 ASSERTALWAYS(!homo.empty());
00338 ASSERTALWAYS(!lumo.empty());
00339 if (!eigInterval.cover((1 + template_blas_sqrt(1.0 + 4.0 * XmX2EuclNorm.low())) / 2) &&
00340 !eigInterval.cover((1 - template_blas_sqrt(1.0 + 4.0 * XmX2EuclNorm.low())) / 2)){
00341
00342 Interval<Treal> homoTmp =
00343 (tmp + (Treal)1.0) / (Treal)2.0;
00344 Interval<Treal> lumoTmp =
00345 ((tmp * (Treal)(-1.0)) + (Treal)1.0) / (Treal)2.0;
00346 ASSERTALWAYS(!(Interval<Treal>::intersect(homo, homoTmp).empty() &&
00347 Interval<Treal>::intersect(lumo, lumoTmp).empty()));
00348 if (Interval<Treal>::intersect(homo, homoTmp).empty()) {
00349
00350 if (eigVecPtr)
00351 lumoWasComputed = 1;
00352 lumo.intersect(lumoTmp);
00353 }
00354 if (Interval<Treal>::intersect(lumo, lumoTmp).empty()) {
00355
00356 if (eigVecPtr)
00357 homoWasComputed = 1;
00358 homo.intersect(homoTmp);
00359 }
00360 checkIntervals("PuriStepInfo::improveHomoLumo C");
00361 }
00362 }
00363 }
00364
00365 template<typename Treal, typename Tvector, typename TdebugPolicy>
00366 void PuriStepInfo<Treal, Tvector, TdebugPolicy>::
00367 exchangeInfoWithNext(PuriStepInfo<Treal, Tvector, TdebugPolicy> & next) {
00368 Interval<Treal> zeroOneInt(0.0,1.0);
00369
00370 next.checkIntervals("PuriStepInfo::exchangeInfoWithNext A");
00371
00372
00373 Interval<Treal> homoForNext(homo);
00374 Interval<Treal> lumoForNext(lumo);
00375 Interval<Treal> eigIntForNext(eigInterval);
00376 homoForNext.puriStep(poly);
00377 lumoForNext.puriStep(poly);
00378 eigIntForNext.puriStep(poly);
00379 ASSERTALWAYS(!homoForNext.empty());
00380 ASSERTALWAYS(!lumoForNext.empty());
00381 ASSERTALWAYS(!eigIntForNext.empty());
00382
00383
00384
00385 homoForNext.increase(eigAccLoss);
00386 lumoForNext.increase(eigAccLoss);
00387 eigIntForNext.increase(eigAccLoss);
00388
00389 homoForNext.increase(next.actualThresh);
00390 lumoForNext.increase(next.actualThresh);
00391 homoForNext.intersect(zeroOneInt);
00392 lumoForNext.intersect(zeroOneInt);
00393 eigIntForNext.increase(next.actualThresh);
00394 next.improveEigInterval(eigIntForNext);
00395 next.improveHomoLumo(homoForNext, lumoForNext);
00396
00397
00398
00399
00400
00401
00402
00403 Interval<Treal> homoTmp(next.homo);
00404 Interval<Treal> lumoTmp(next.lumo);
00405 ASSERTALWAYS(!homoTmp.empty());
00406 ASSERTALWAYS(!lumoTmp.empty());
00407
00408 homoTmp.increase(next.actualThresh);
00409 lumoTmp.increase(next.actualThresh);
00410 homoTmp.intersect(zeroOneInt);
00411 lumoTmp.intersect(zeroOneInt);
00412 homoTmp.invPuriStep(poly);
00413 lumoTmp.invPuriStep(poly);
00414
00415
00416
00417 homoTmp.increase(eigAccLoss);
00418 lumoTmp.increase(eigAccLoss);
00419
00420 this->improveHomoLumo(homoTmp, lumoTmp);
00421 }
00422
00423 template<typename Treal, typename Tvector, typename TdebugPolicy>
00424 Treal PuriStepInfo<Treal, Tvector, TdebugPolicy>::subspaceError() const {
00425 Interval<Treal> gap = Interval<Treal>(lumo.upp(),homo.low());
00426 if (actualThresh >= gap.length())
00427 return 1.0;
00428 else {
00429 Treal error = actualThresh / (gap.length() - actualThresh);
00430 return error < 1.0 ? error : (Treal)1.0;
00431 }
00432 }
00433
00434 template<typename Treal, typename Tvector, typename TdebugPolicy>
00435 void PuriStepInfo<Treal, Tvector, TdebugPolicy>::
00436 improveEigInterval(Interval<Treal> const eInt) {
00437 eigInterval.intersect(eInt);
00438 Treal delta1 = eigInterval.upp() - 1;
00439 Treal delta0 = -eigInterval.low();
00440 delta = delta1 > delta0 ? delta1 : delta0;
00441 this->computen0n1();
00442 }
00443
00444
00445 template<typename Treal, typename Tvector, typename TdebugPolicy>
00446 void PuriStepInfo<Treal, Tvector, TdebugPolicy>::computen0n1() {
00447 Treal beta = 0.5;
00448
00449 if (XmX2EuclNorm.upp() < 1 / (Treal)4)
00450 beta = (1 + template_blas_sqrt(1 - 4 * XmX2EuclNorm.upp())) / 2;
00451 n1 = (traceX2 - delta * (1 - beta) * n -
00452 (1 - delta - beta) * traceX) /
00453 ((1 + 2 * delta) * (delta + beta));
00454 n0 = (traceX2 + beta * (1 + delta) * n -
00455 (1 + delta + beta) * traceX) /
00456 ((1 + 2 * delta) * (delta + beta));
00457 if (n1 > nocc -1 &&
00458 n0 > n - nocc - 1)
00459 correctOccupation = 1;
00460 }
00461
00466 template<typename Treal, typename Tvector, typename TdebugPolicy>
00467 void PuriStepInfo<Treal, Tvector, TdebugPolicy>::computeEigAccLoss() {
00468 Treal nnzPerRowX = nnzX / (Treal)n;
00469 Treal maxAbsErrPerElX2 = getRelPrecision<Treal>() * nnzPerRowX;
00470
00471
00472
00473
00474
00475
00476 eigAccLoss = maxAbsErrPerElX2 * template_blas_sqrt((Treal)nnzX2);
00477 ASSERTALWAYS(eigAccLoss >= 0);
00478 }
00479
00480
00481
00482 }
00483 #undef __ID__
00484 #endif