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_PURIFICATION
00037 #define MAT_PURIFICATION
00038 #include <math.h>
00039 #include <iomanip>
00040 #define __ID__ "$Id: Purification.h 4437 2012-07-05 09:01:18Z elias $"
00041 namespace mat {
00042 template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00043 class Purification :public TdebugPolicy {
00044 public:
00045 typedef typename Tmatrix::VectorType VectorType;
00046 Purification(Tmatrix& M,
00049 normType const normXmX2,
00050 PuriInfo<Treal, VectorType, TdebugPolicy>& info
00057 );
00058 void step();
00059 void purify();
00060 protected:
00061 Tmatrix & X;
00062 Tmatrix X2;
00063 normType const normTruncation;
00064 normType const normXmX2;
00065 PuriInfo<Treal, VectorType, TdebugPolicy>& info;
00066 int niter;
00067
00068 void stepComputeInfo(PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep);
00069
00070 private:
00071 };
00072
00073 template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00074 Purification<Treal, Tmatrix, TdebugPolicy>::
00075 Purification(Tmatrix& M,
00076 normType const normXmX2_inp,
00077 PuriInfo<Treal, VectorType, TdebugPolicy>& infop)
00078 :X(M), X2(M), normTruncation( infop.getNormForTruncation() ), normXmX2(normXmX2_inp), info(infop),niter(0) {
00079 Time tStepTotal;
00080 tStepTotal.tic();
00081
00082 Treal lmin(0);
00083 Treal lmax(0);
00084 PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep = info.getNext();
00085 currentStep.improveEigInterval(Interval<Treal>(0.0,1.0));
00086 Interval<Treal> eigFInt = info.getEigFInterval();
00087 lmin = eigFInt.low();
00088 lmax = eigFInt.upp();
00089 X.add_identity(-lmax);
00090 X *= ((Treal)1.0 / (lmin - lmax));
00091
00092 Interval<Treal> homo = info.getHomoF();
00093 Interval<Treal> lumo = info.getLumoF();
00094 homo = (homo - lmax) / (lmin - lmax);
00095 lumo = (lumo - lmax) / (lmin - lmax);
00096
00097 Treal chosenThresh = info.getOptimalThresh();
00098 currentStep.setChosenThresh(chosenThresh);
00099
00100 currentStep.setMemUsageBeforeTrunc();
00101 Time t;
00102 t.tic();
00103 Treal actualThresh = X.thresh(chosenThresh, normTruncation);
00104 currentStep.setTimeThresh(t.toc());
00105 currentStep.setActualThresh(actualThresh);
00106 if (homo.empty() || lumo.empty()) {
00107 throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>::"
00108 "Purification(Tmatrix&, normType const, "
00109 "PuriInfo<Treal, VectorType, TdebugPolicy>&): "
00110 "HOMO or LUMO empty in purification constructor. ");
00111 }
00112
00113 homo.increase(actualThresh);
00114 lumo.increase(actualThresh);
00115 stepComputeInfo(currentStep);
00116 currentStep.improveHomoLumo(homo,lumo);
00117 currentStep.setTimeTotal(tStepTotal.toc());
00118 }
00120 template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00121 void Purification<Treal, Tmatrix, TdebugPolicy>::step() {
00122 Time tStepTotal;
00123 tStepTotal.tic();
00124 PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep = info.getNext();
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 if (info(niter).getPoly()) {
00139
00140 X2 *= (Treal)-1.0;
00141 X2 += (Treal)2.0 * X;
00142
00143 }
00144
00145
00146
00147
00148
00149 X2.transfer(X);
00150
00151
00152
00153
00154
00155 Treal chosenThresh = info.getOptimalThresh();
00156 currentStep.setChosenThresh(chosenThresh);
00157
00158 currentStep.setMemUsageBeforeTrunc();
00159 Time t;
00160 t.tic();
00161 currentStep.setActualThresh(X.thresh(chosenThresh, normTruncation));
00162 currentStep.setTimeThresh(t.toc());
00163 stepComputeInfo(currentStep);
00164 if (niter > 5 &&
00165 currentStep.getHomo().low() > 0.9 &&
00166 currentStep.getLumo().upp() < 0.1 &&
00167 info(niter-5).getHomo().low() > 0.9 &&
00168 info(niter-5).getLumo().upp() < 0.1 &&
00169 ((1 - currentStep.getHomo().low()) >
00170 (1 - info(niter-5).getHomo().low()) -
00171 currentStep.getEigAccLoss() ||
00172 currentStep.getLumo().upp() >
00173 info(niter-5).getLumo().upp() -
00174 currentStep.getEigAccLoss())) {
00175 throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>"
00176 "::step(): HOMO-LUMO gap has not increased in the "
00177 "last five iterations. Desired accuracy too high for"
00178 "current floating point precision?");
00179 }
00180
00181 ++niter;
00182 currentStep.setTimeTotal(tStepTotal.toc());
00183 }
00184
00185 template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00186 void Purification<Treal, Tmatrix, TdebugPolicy>::purify() {
00187
00188 while (niter < info.getMaxSteps() - 1 && !info.converged() ) {
00189 step();
00190 }
00191 if ( info.converged() ) {
00192
00193
00194 info.forceCorrectOccupation();
00195 info.improveCorrectOccupation();
00196 info.improveInfo();
00197 info.improveHomoLumoF();
00198 }
00199 }
00200
00201 template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00202 void Purification<Treal, Tmatrix, TdebugPolicy>::
00203 stepComputeInfo(PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep) {
00204 Treal const XmX2ENIsSmallValue = 0.207106781186547;
00205 Treal const ONE = 1.0;
00206
00207 Time t;
00208 t.tic();
00209 X2 = ONE * X * X;
00210 currentStep.setTimeSquare(t.toc());
00211
00212 currentStep.setNnzX(X.nnz());
00213 currentStep.setNnzX2(X2.nnz());
00214 currentStep.computeEigAccLoss();
00215 currentStep.computeExactValues(X, X2);
00216 currentStep.setTraceX(X.trace());
00217 currentStep.setTraceX2(X2.trace());
00218
00219
00220
00221
00222
00223
00224 typename Tmatrix::VectorType * eigVecPtr = new typename Tmatrix::VectorType;
00225 Treal diffAcc;
00226 Interval<Treal> XmX2EN;
00227 t.tic();
00228 if (info.ShouldComputeXmX2EuclNormAccurately(diffAcc)) {
00229 if (normXmX2 == euclNorm)
00230 XmX2EN = Tmatrix::euclDiffIfSmall(X, X2,
00231 diffAcc,
00232 XmX2ENIsSmallValue,
00233 eigVecPtr);
00234 else
00235 XmX2EN = Tmatrix::diffIfSmall(X, X2,
00236 normXmX2, diffAcc,
00237 XmX2ENIsSmallValue);
00238 }
00239 else {
00240 XmX2EN = Tmatrix::diffIfSmall(X, X2,
00241 frobNorm, diffAcc,
00242 XmX2ENIsSmallValue);
00243 }
00244 currentStep.setTimeXmX2Norm(t.toc());
00245
00246 if (!eigVecPtr->is_empty())
00247 currentStep.setEigVecPtr(eigVecPtr);
00248 else
00249 delete eigVecPtr;
00250
00251 XmX2EN.increase(currentStep.getEigAccLoss());
00252 Interval<Treal> zeroInt(0.0, XmX2EN.upp());
00253 XmX2EN.intersect(zeroInt);
00254
00255
00256
00257 currentStep.setXmX2EuclNorm(XmX2EN);
00258 currentStep.checkIntervals("Purification::stepComputeInfo after setXmX2EuclNorm.");
00259
00260 Treal tmpVal = template_blas_sqrt(1 + 4 * XmX2EN.upp());
00261 currentStep.improveEigInterval
00262 (Interval<Treal>((1 - tmpVal) / 2, (1 + tmpVal) / 2));
00263
00264 currentStep.checkIntervals("Purification::stepComputeInfo B");
00265
00266 info.improveInfo();
00267 if (currentStep.getEigInterval().length() > 1.5)
00268 throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>"
00269 "::step() : "
00270 "Eigenvalues moved to far from [0, 1] interval.");
00271
00272 currentStep.setPoly();
00273 currentStep.checkIntervals("Purification::stepComputeInfo C");
00274 }
00275
00276
00277
00278
00279 }
00280 #undef __ID__
00281 #endif