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_PERTURBATION
00039 #define MAT_PERTURBATION
00040 namespace per {
00041 template<typename Treal, typename Tmatrix, typename Tvector>
00042 class Perturbation {
00043 public:
00044 Perturbation(std::vector<Tmatrix *> const & F,
00046 std::vector<Tmatrix *> & D,
00048 mat::Interval<Treal> const & gap,
00049 mat::Interval<Treal> const & allEigs,
00055 Treal const deltaMax,
00056 Treal const errorTol,
00057 mat::normType const norm,
00058 Tvector & vect
00059 );
00060 void perturb() {
00061 dryRun();
00062 run();
00063 }
00064
00065 void checkIdempotencies(std::vector<Treal> & idemErrors);
00066 template<typename TmatNoSymm>
00067 void checkCommutators(std::vector<Treal> & commErrors,
00068 TmatNoSymm const & dummyMat);
00069 void checkMaxSubspaceError(Treal & subsError);
00070
00071 protected:
00072
00073 std::vector<Tmatrix *> const & F;
00074 std::vector<Tmatrix *> & X;
00075 mat::Interval<Treal> gap;
00076 mat::Interval<Treal> const & allEigs;
00077 Treal deltaMax;
00078 Treal errorTol;
00079 mat::normType const norm;
00080 Tvector & vect;
00081
00082
00083 int nIter;
00084 std::vector<Treal> threshVal;
00085 std::vector<Treal> sigma;
00086
00097 void dryRun();
00098 void run();
00099 private:
00100
00101 };
00102
00103 template<typename Treal, typename Tmatrix, typename Tvector>
00104 Perturbation<Treal, Tmatrix, Tvector>::
00105 Perturbation(std::vector<Tmatrix *> const & F_in,
00106 std::vector<Tmatrix *> & X_in,
00107 mat::Interval<Treal> const & gap_in,
00108 mat::Interval<Treal> const & allEigs_in,
00109 Treal const deltaMax_in,
00110 Treal const errorTol_in,
00111 mat::normType const norm_in,
00112 Tvector & vect_in)
00113 : F(F_in), X(X_in), gap(gap_in), allEigs(allEigs_in),
00114 deltaMax(deltaMax_in), errorTol(errorTol_in), norm(norm_in),
00115 vect(vect_in) {
00116 if (!X.empty())
00117 throw "Perturbation constructor: D vector is expected to be empty (size==0)";
00118 for (unsigned int iMat = 0; iMat < F.size(); ++iMat)
00119 X.push_back(new Tmatrix(*F[iMat]));
00120
00121 Treal lmin = allEigs.low();
00122 Treal lmax = allEigs.upp();
00123
00124
00125 typename std::vector<Tmatrix *>::iterator matIt = X.begin();
00126
00127 (*matIt)->add_identity(-lmax);
00128 *(*matIt) *= ((Treal)1.0 / (lmin - lmax));
00129 matIt++;
00130
00131 for ( ; matIt != X.end(); matIt++ )
00132 *(*matIt) *= ((Treal)-1.0 / (lmin - lmax));
00133
00134 gap = (gap - lmax) / (lmin - lmax);
00135 }
00136
00137 template<typename Treal, typename Tmatrix, typename Tvector>
00138 void Perturbation<Treal, Tmatrix, Tvector>::dryRun() {
00139 Treal errorTolPerIter;
00140 int nIterGuess = 0;
00141 nIter = 1;
00142 Treal lumo;
00143 Treal homo;
00144 Treal m;
00145 Treal g;
00146 while (nIterGuess < nIter) {
00147 nIterGuess++;
00148 errorTolPerIter = 0.5 * errorTol /nIterGuess;
00149 nIter = 0;
00150 mat::Interval<Treal> gapTmp(gap);
00151 sigma.resize(0);
00152 threshVal.resize(0);
00153 while (gapTmp.low() > 0.5 * errorTol || gapTmp.upp() < 0.5 * errorTol) {
00154 lumo = gapTmp.low();
00155 homo = gapTmp.upp();
00156 m = gapTmp.midPoint();
00157 g = gapTmp.length();
00158 if (m > 0.5) {
00159 lumo = lumo*lumo;
00160 homo = homo*homo;
00161 sigma.push_back(-1);
00162 }
00163 else {
00164 lumo = 2*lumo - lumo*lumo;
00165 homo = 2*homo - homo*homo;
00166 sigma.push_back(1);
00167 }
00168
00169 Treal forceConvThresh = template_blas_fabs(1-2*m) * g / 10;
00170
00171
00172 Treal subspaceThresh = errorTolPerIter * (homo-lumo) / (1+errorTolPerIter);
00173
00174 threshVal.push_back(forceConvThresh < subspaceThresh ?
00175 forceConvThresh : subspaceThresh);
00176 homo -= threshVal.back();
00177 lumo += threshVal.back();
00178 gapTmp = mat::Interval<Treal>(lumo, homo);
00179 if (gapTmp.empty())
00180 throw "Perturbation<Treal, Tmatrix, Tvector>::dryRun() : Perturbation iterations will fail to converge; Gap is too small or desired accuracy too high.";
00181 nIter++;
00182 }
00183 }
00184
00185 }
00186
00187 template<typename Treal, typename Tmatrix, typename Tvector>
00188 void Perturbation<Treal, Tmatrix, Tvector>::run() {
00189 Treal const ONE = 1.0;
00190 mat::SizesAndBlocks rowsCopy;
00191 X.front()->getRows(rowsCopy);
00192 mat::SizesAndBlocks colsCopy;
00193 X.front()->getCols(colsCopy);
00194 Tmatrix tmpMat;
00195
00196 int nMatrices;
00197 Treal threshValPerOrder;
00198 Treal chosenThresh;
00199 for (int iter = 0; iter < nIter; iter++) {
00200 std::cout<<"\n\nInside outer loop iter = "<<iter
00201 <<" nIter = "<<nIter
00202 <<" sigma = "<<sigma[iter]<<std::endl;
00203
00204 X.push_back(new Tmatrix);
00205 nMatrices = X.size();
00206 X[nMatrices-1]->resetSizesAndBlocks(rowsCopy, colsCopy);
00207
00208 threshValPerOrder = threshVal[iter] / nMatrices;
00209
00210 std::cout<<"Entering inner loop nMatrices = "<<nMatrices<<std::endl;
00211 for (int j = nMatrices-1 ; j >= 0 ; --j) {
00212 std::cout<<"Inside inner loop j = "<<j<<std::endl;
00213 std::cout<<"X[j]->eucl() (before compute) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
00214 std::cout<<"X[j]->frob() (before compute) = "<<X[j]->frob()<<std::endl;
00215 tmpMat = Treal(Treal(1.0)+sigma[iter]) * (*X[j]);
00216 std::cout<<"tmpMat.eucl() (before for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl;
00217 std::cout<<"tmpMat.frob() (before for) = "<<tmpMat.frob()<<std::endl;
00218 for (int k = 0; k <= j; k++) {
00219
00220 Tmatrix::ssmmUpperTriangleOnly(-sigma[iter], *X[k], *X[j-k],
00221 ONE, tmpMat);
00222 }
00223 std::cout<<"tmpMat.eucl() (after for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl;
00224 *X[j] = tmpMat;
00225
00226
00227 chosenThresh = threshValPerOrder / pow(deltaMax, Treal(j));
00228 std::cout<<"X[j]->eucl() (before thresh) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
00229 std::cout<<"Chosen thresh: "<<chosenThresh<<std::endl;
00230 Treal actualThresh = X[j]->thresh(chosenThresh, vect, norm);
00231 std::cout<<"X[j]->eucl() (after thresh) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
00232 #if 1
00233
00234
00235
00236 if (*X[j] == 0 && int(X.size()) == j+1) {
00237 std::cout<<"DELETION: j = "<<j<<" X.size() = "<<X.size()
00238 <<" X[j] = "<<X[j]<< " X[j]->frob() = "<<X[j]->frob()
00239 <<std::endl;
00240 delete X[j];
00241 X.pop_back();
00242 }
00243 else
00244 std::cout<<"NO DELETION: j = "<<j<<" X.size() = "<<X.size()
00245 <<" X[j] = "<<X[j]<< " X[j]->frob() = "<<X[j]->frob()
00246 <<std::endl;
00247 #endif
00248
00249 }
00250 }
00251 }
00252
00253
00254 template<typename Treal, typename Tmatrix, typename Tvector>
00255 void Perturbation<Treal, Tmatrix, Tvector>::
00256 checkIdempotencies(std::vector<Treal> & idemErrors) {
00257 Tmatrix tmpMat;
00258 Treal const ONE = 1.0;
00259 unsigned int j;
00260 for (unsigned int m = 0; m < X.size(); ++m) {
00261 tmpMat = (-ONE) * (*X[m]);
00262 for (unsigned int i = 0; i <= m; ++i) {
00263 j = m - i;
00264
00265 Tmatrix::ssmmUpperTriangleOnly(ONE, *X[i], *X[j], ONE, tmpMat);
00266 }
00267
00268 idemErrors.push_back(tmpMat.eucl(vect,1e-10));
00269 }
00270 }
00271
00272 template<typename Treal, typename Tmatrix, typename Tvector>
00273 template<typename TmatNoSymm>
00274 void Perturbation<Treal, Tmatrix, Tvector>::
00275 checkCommutators(std::vector<Treal> & commErrors,
00276 TmatNoSymm const & dummyMat) {
00277 mat::SizesAndBlocks rowsCopy;
00278 X.front()->getRows(rowsCopy);
00279 mat::SizesAndBlocks colsCopy;
00280 X.front()->getCols(colsCopy);
00281 TmatNoSymm tmpMat;
00282 tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy);
00283 Treal const ONE = 1.0;
00284 unsigned int j;
00285 for (unsigned int m = 0; m < X.size(); ++m) {
00286 tmpMat = 0;
00287 std::cout<<"New loop\n";
00288 for (unsigned int i = 0; i <= m && i < F.size(); ++i) {
00289 j = m - i;
00290 std::cout<<i<<", "<<j<<std::endl;
00291
00292 tmpMat += ONE * (*F[i]) * (*X[j]);
00293 tmpMat += -ONE * (*X[j]) * (*F[i]);
00294 }
00295
00296 commErrors.push_back(tmpMat.frob());
00297 }
00298 }
00299
00300
00301 template<typename Treal, typename Tmatrix, typename Tvector>
00302 void Perturbation<Treal, Tmatrix, Tvector>::
00303 checkMaxSubspaceError(Treal & subsError) {
00304 Treal const ONE = 1.0;
00305 Tmatrix XdeltaMax(*F.front());
00306 for (unsigned int ind = 1; ind < F.size(); ++ind)
00307 XdeltaMax += pow(deltaMax, Treal(ind)) * (*F[ind]);
00308
00309 Treal lmin = allEigs.low();
00310 Treal lmax = allEigs.upp();
00311
00312 XdeltaMax.add_identity(-lmax);
00313 XdeltaMax *= ((Treal)1.0 / (lmin - lmax));
00314
00315 Tmatrix X2;
00316 for (int iter = 0; iter < nIter; iter++) {
00317 X2 = ONE * XdeltaMax * XdeltaMax;
00318 if (sigma[iter] == Treal(1.0)) {
00319 XdeltaMax *= 2.0;
00320 XdeltaMax -= X2;
00321 }
00322 else {
00323 XdeltaMax = X2;
00324 }
00325 }
00326
00327 Tmatrix DdeltaMax(*X.front());
00328 for (unsigned int ind = 1; ind < X.size(); ++ind)
00329 DdeltaMax += pow(deltaMax, Treal(ind)) * (*X[ind]);
00330 subsError = Tmatrix::eucl_diff(XdeltaMax,DdeltaMax,
00331 vect, errorTol *1e-2);
00332 }
00333
00334
00335
00336 }
00337 #endif
00338