36 #ifndef MAT_PERTURBATION
37 #define MAT_PERTURBATION
39 template<
typename Treal,
typename Tmatrix,
typename Tvector>
44 std::vector<Tmatrix *> & D,
64 template<
typename TmatNoSymm>
66 TmatNoSymm
const & dummyMat);
71 std::vector<Tmatrix *>
const &
F;
72 std::vector<Tmatrix *> &
X;
101 template<
typename Treal,
typename Tmatrix,
typename Tvector>
104 std::vector<Tmatrix *> & X_in,
107 Treal
const deltaMax_in,
108 Treal
const errorTol_in,
111 : F(F_in), X(X_in), gap(gap_in), allEigs(allEigs_in),
112 deltaMax(deltaMax_in), errorTol(errorTol_in), norm(norm_in),
115 throw "Perturbation constructor: D vector is expected to be empty (size==0)";
116 for (
unsigned int iMat = 0; iMat <
F.size(); ++iMat)
117 X.push_back(
new Tmatrix(*
F[iMat]));
123 typename std::vector<Tmatrix *>::iterator matIt =
X.begin();
125 (*matIt)->add_identity(-lmax);
126 *(*matIt) *= ((Treal)1.0 / (lmin - lmax));
129 for ( ; matIt !=
X.end(); matIt++ )
130 *(*matIt) *= ((Treal)-1.0 / (lmin - lmax));
132 gap = (
gap - lmax) / (lmin - lmax);
135 template<
typename Treal,
typename Tmatrix,
typename Tvector>
137 Treal errorTolPerIter;
144 while (nIterGuess < nIter) {
146 errorTolPerIter = 0.5 * errorTol /nIterGuess;
151 while (gapTmp.
low() > 0.5 * errorTol || gapTmp.
upp() < 0.5 * errorTol) {
162 lumo = 2*lumo - lumo*lumo;
163 homo = 2*homo - homo*homo;
170 Treal subspaceThresh = errorTolPerIter * (homo-lumo) / (1+errorTolPerIter);
172 threshVal.push_back(forceConvThresh < subspaceThresh ?
173 forceConvThresh : subspaceThresh);
174 homo -= threshVal.back();
175 lumo += threshVal.back();
178 throw "Perturbation<Treal, Tmatrix, Tvector>::dryRun() : Perturbation iterations will fail to converge; Gap is too small or desired accuracy too high.";
185 template<
typename Treal,
typename Tmatrix,
typename Tvector>
187 Treal
const ONE = 1.0;
189 X.front()->getRows(rowsCopy);
191 X.front()->getCols(colsCopy);
195 Treal threshValPerOrder;
197 for (
int iter = 0; iter < nIter; iter++) {
198 std::cout<<
"\n\nInside outer loop iter = "<<iter
200 <<
" sigma = "<<sigma[iter]<<std::endl;
202 X.push_back(
new Tmatrix);
203 nMatrices = X.size();
204 X[nMatrices-1]->resetSizesAndBlocks(rowsCopy, colsCopy);
206 threshValPerOrder = threshVal[iter] / nMatrices;
208 std::cout<<
"Entering inner loop nMatrices = "<<nMatrices<<std::endl;
209 for (
int j = nMatrices-1 ; j >= 0 ; --j) {
210 std::cout<<
"Inside inner loop j = "<<j<<std::endl;
211 std::cout<<
"X[j]->eucl() (before compute) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
212 std::cout<<
"X[j]->frob() (before compute) = "<<X[j]->frob()<<std::endl;
213 tmpMat = Treal(Treal(1.0)+sigma[iter]) * (*X[j]);
214 std::cout<<
"tmpMat.eucl() (before for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl;
215 std::cout<<
"tmpMat.frob() (before for) = "<<tmpMat.frob()<<std::endl;
216 for (
int k = 0; k <= j; k++) {
218 Tmatrix::ssmmUpperTriangleOnly(-sigma[iter], *X[k], *X[j-k],
221 std::cout<<
"tmpMat.eucl() (after for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl;
225 chosenThresh = threshValPerOrder / pow(deltaMax, Treal(j));
226 std::cout<<
"X[j]->eucl() (before thresh) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
227 std::cout<<
"Chosen thresh: "<<chosenThresh<<std::endl;
228 Treal actualThresh = X[j]->thresh(chosenThresh, vect, norm);
229 std::cout<<
"X[j]->eucl() (after thresh) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
234 if (*X[j] == 0 &&
int(X.size()) == j+1) {
235 std::cout<<
"DELETION: j = "<<j<<
" X.size() = "<<X.size()
236 <<
" X[j] = "<<X[j]<<
" X[j]->frob() = "<<X[j]->frob()
242 std::cout<<
"NO DELETION: j = "<<j<<
" X.size() = "<<X.size()
243 <<
" X[j] = "<<X[j]<<
" X[j]->frob() = "<<X[j]->frob()
252 template<
typename Treal,
typename Tmatrix,
typename Tvector>
256 Treal
const ONE = 1.0;
258 for (
unsigned int m = 0; m < X.size(); ++m) {
259 tmpMat = (-ONE) * (*X[m]);
260 for (
unsigned int i = 0; i <= m; ++i) {
263 Tmatrix::ssmmUpperTriangleOnly(ONE, *X[i], *X[j], ONE, tmpMat);
266 idemErrors.push_back(tmpMat.eucl(vect,1e-10));
270 template<
typename Treal,
typename Tmatrix,
typename Tvector>
271 template<
typename TmatNoSymm>
274 TmatNoSymm
const & dummyMat) {
276 X.front()->getRows(rowsCopy);
278 X.front()->getCols(colsCopy);
280 tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy);
281 Treal
const ONE = 1.0;
283 for (
unsigned int m = 0; m < X.size(); ++m) {
285 std::cout<<
"New loop\n";
286 for (
unsigned int i = 0; i <= m && i < F.size(); ++i) {
288 std::cout<<i<<
", "<<j<<std::endl;
290 tmpMat += ONE * (*F[i]) * (*X[j]);
291 tmpMat += -ONE * (*X[j]) * (*F[i]);
294 commErrors.push_back(tmpMat.frob());
299 template<
typename Treal,
typename Tmatrix,
typename Tvector>
302 Treal
const ONE = 1.0;
303 Tmatrix XdeltaMax(*F.front());
304 for (
unsigned int ind = 1; ind < F.size(); ++ind)
305 XdeltaMax += pow(deltaMax, Treal(ind)) * (*F[ind]);
307 Treal lmin = allEigs.low();
308 Treal lmax = allEigs.upp();
310 XdeltaMax.add_identity(-lmax);
311 XdeltaMax *= ((Treal)1.0 / (lmin - lmax));
314 for (
int iter = 0; iter < nIter; iter++) {
315 X2 = ONE * XdeltaMax * XdeltaMax;
316 if (sigma[iter] == Treal(1.0)) {
325 Tmatrix DdeltaMax(*X.front());
326 for (
unsigned int ind = 1; ind < X.size(); ++ind)
327 DdeltaMax += pow(deltaMax, Treal(ind)) * (*X[ind]);
328 subsError = Tmatrix::eucl_diff(XdeltaMax,DdeltaMax,
329 vect, errorTol *1e-2);
Tvector & vect
Definition: Perturbation.h:78
void dryRun()
Dry run to obtain some needed numbers.
Definition: Perturbation.h:136
Treal upp() const
Definition: Interval.h:143
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
std::vector< Tmatrix * > const & F
Definition: Perturbation.h:71
Perturbation(std::vector< Tmatrix * > const &F, std::vector< Tmatrix * > &D, mat::Interval< Treal > const &gap, mat::Interval< Treal > const &allEigs, Treal const deltaMax, Treal const errorTol, mat::normType const norm, Tvector &vect)
Definition: Perturbation.h:103
mat::Interval< Treal > const & allEigs
Definition: Perturbation.h:74
void run()
Definition: Perturbation.h:186
bool empty() const
Definition: Interval.h:49
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
Treal midPoint() const
Definition: Interval.h:113
Treal template_blas_fabs(Treal x)
Definition: Interval.h:44
void checkCommutators(std::vector< Treal > &commErrors, TmatNoSymm const &dummyMat)
Definition: Perturbation.h:273
Treal errorTol
Definition: Perturbation.h:76
void perturb()
Definition: Perturbation.h:58
std::vector< Tmatrix * > & X
Definition: Perturbation.h:72
void checkMaxSubspaceError(Treal &subsError)
Definition: Perturbation.h:301
mat::normType const norm
Definition: Perturbation.h:77
Treal low() const
Definition: Interval.h:142
Definition: Perturbation.h:40
std::vector< Treal > sigma
Definition: Perturbation.h:83
void checkIdempotencies(std::vector< Treal > &idemErrors)
Definition: Perturbation.h:254
Treal deltaMax
Definition: Perturbation.h:75
int nIter
Definition: Perturbation.h:81
std::vector< Treal > threshVal
Definition: Perturbation.h:82
mat::Interval< Treal > gap
Definition: Perturbation.h:73
normType
Definition: matInclude.h:135