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