38#ifndef MAT_PERTURBATION
39#define MAT_PERTURBATION
41 template<
typename Treal,
typename Tmatrix,
typename Tvector>
46 std::vector<Tmatrix *> & D,
66 template<
typename TmatNoSymm>
68 TmatNoSymm
const & dummyMat);
73 std::vector<Tmatrix *>
const &
F;
74 std::vector<Tmatrix *> &
X;
103 template<
typename Treal,
typename Tmatrix,
typename Tvector>
106 std::vector<Tmatrix *> & X_in,
109 Treal
const deltaMax_in,
110 Treal
const errorTol_in,
117 throw "Perturbation constructor: D vector is expected to be empty (size==0)";
118 for (
unsigned int iMat = 0; iMat <
F.size(); ++iMat)
119 X.push_back(
new Tmatrix(*
F[iMat]));
125 typename std::vector<Tmatrix *>::iterator matIt =
X.begin();
127 (*matIt)->add_identity(-lmax);
128 *(*matIt) *= ((Treal)1.0 / (lmin - lmax));
131 for ( ; matIt !=
X.end(); matIt++ )
132 *(*matIt) *= ((Treal)-1.0 / (lmin - lmax));
134 gap = (
gap - lmax) / (lmin - lmax);
137 template<
typename Treal,
typename Tmatrix,
typename Tvector>
139 Treal errorTolPerIter;
146 while (nIterGuess <
nIter) {
148 errorTolPerIter = 0.5 *
errorTol /nIterGuess;
164 lumo = 2*lumo - lumo*lumo;
165 homo = 2*homo - homo*homo;
172 Treal subspaceThresh = errorTolPerIter * (homo-lumo) / (1+errorTolPerIter);
174 threshVal.push_back(forceConvThresh < subspaceThresh ?
175 forceConvThresh : subspaceThresh);
180 throw "Perturbation<Treal, Tmatrix, Tvector>::dryRun() : Perturbation iterations will fail to converge; Gap is too small or desired accuracy too high.";
187 template<
typename Treal,
typename Tmatrix,
typename Tvector>
189 Treal
const ONE = 1.0;
191 X.front()->getRows(rowsCopy);
193 X.front()->getCols(colsCopy);
197 Treal threshValPerOrder;
199 for (
int iter = 0; iter <
nIter; iter++) {
200 std::cout<<
"\n\nInside outer loop iter = "<<iter
202 <<
" sigma = "<<
sigma[iter]<<std::endl;
204 X.push_back(
new Tmatrix);
205 nMatrices =
X.size();
206 X[nMatrices-1]->resetSizesAndBlocks(rowsCopy, colsCopy);
208 threshValPerOrder =
threshVal[iter] / nMatrices;
210 std::cout<<
"Entering inner loop nMatrices = "<<nMatrices<<std::endl;
211 for (
int j = nMatrices-1 ; j >= 0 ; --j) {
212 std::cout<<
"Inside inner loop j = "<<j<<std::endl;
213 std::cout<<
"X[j]->eucl() (before compute) = "<<
X[j]->eucl(
vect,1e-7)<<std::endl;
214 std::cout<<
"X[j]->frob() (before compute) = "<<
X[j]->frob()<<std::endl;
215 tmpMat = Treal(Treal(1.0)+
sigma[iter]) * (*
X[j]);
216 std::cout<<
"tmpMat.eucl() (before for) = "<<tmpMat.eucl(
vect,1e-7)<<std::endl;
217 std::cout<<
"tmpMat.frob() (before for) = "<<tmpMat.frob()<<std::endl;
218 for (
int k = 0; k <= j; k++) {
220 Tmatrix::ssmmUpperTriangleOnly(-
sigma[iter], *
X[k], *
X[j-k],
223 std::cout<<
"tmpMat.eucl() (after for) = "<<tmpMat.eucl(
vect,1e-7)<<std::endl;
227 chosenThresh = threshValPerOrder / pow(
deltaMax, Treal(j));
228 std::cout<<
"X[j]->eucl() (before thresh) = "<<
X[j]->eucl(
vect,1e-7)<<std::endl;
229 std::cout<<
"Chosen thresh: "<<chosenThresh<<std::endl;
230 Treal actualThresh =
X[j]->thresh(chosenThresh,
vect,
norm);
231 std::cout<<
"X[j]->eucl() (after thresh) = "<<
X[j]->eucl(
vect,1e-7)<<std::endl;
236 if (*
X[j] == 0 &&
int(
X.size()) == j+1) {
237 std::cout<<
"DELETION: j = "<<j<<
" X.size() = "<<
X.size()
238 <<
" X[j] = "<<
X[j]<<
" X[j]->frob() = "<<
X[j]->frob()
244 std::cout<<
"NO DELETION: j = "<<j<<
" X.size() = "<<
X.size()
245 <<
" X[j] = "<<
X[j]<<
" X[j]->frob() = "<<
X[j]->frob()
254 template<
typename Treal,
typename Tmatrix,
typename Tvector>
258 Treal
const ONE = 1.0;
260 for (
unsigned int m = 0; m <
X.size(); ++m) {
261 tmpMat = (-ONE) * (*
X[m]);
262 for (
unsigned int i = 0; i <= m; ++i) {
265 Tmatrix::ssmmUpperTriangleOnly(ONE, *
X[i], *
X[j], ONE, tmpMat);
268 idemErrors.push_back(tmpMat.eucl(
vect,1e-10));
272 template<
typename Treal,
typename Tmatrix,
typename Tvector>
273 template<
typename TmatNoSymm>
276 TmatNoSymm
const & dummyMat) {
278 X.front()->getRows(rowsCopy);
280 X.front()->getCols(colsCopy);
282 tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy);
283 Treal
const ONE = 1.0;
285 for (
unsigned int m = 0; m <
X.size(); ++m) {
287 std::cout<<
"New loop\n";
288 for (
unsigned int i = 0; i <= m && i <
F.size(); ++i) {
290 std::cout<<i<<
", "<<j<<std::endl;
292 tmpMat += ONE * (*
F[i]) * (*
X[j]);
293 tmpMat += -ONE * (*
X[j]) * (*
F[i]);
296 commErrors.push_back(tmpMat.frob());
301 template<
typename Treal,
typename Tmatrix,
typename Tvector>
304 Treal
const ONE = 1.0;
305 Tmatrix XdeltaMax(*
F.front());
306 for (
unsigned int ind = 1; ind <
F.size(); ++ind)
307 XdeltaMax += pow(
deltaMax, Treal(ind)) * (*
F[ind]);
312 XdeltaMax.add_identity(-lmax);
313 XdeltaMax *= ((Treal)1.0 / (lmin - lmax));
316 for (
int iter = 0; iter <
nIter; iter++) {
317 X2 = ONE * XdeltaMax * XdeltaMax;
318 if (
sigma[iter] == Treal(1.0)) {
327 Tmatrix DdeltaMax(*
X.front());
328 for (
unsigned int ind = 1; ind <
X.size(); ++ind)
329 DdeltaMax += pow(
deltaMax, Treal(ind)) * (*
X[ind]);
330 subsError = Tmatrix::eucl_diff(XdeltaMax,DdeltaMax,
Treal low() const
Definition Interval.h:144
Treal midPoint() const
Definition Interval.h:115
Treal length() const
Returns the length of the interval.
Definition Interval.h:109
Treal upp() const
Definition Interval.h:145
bool empty() const
Definition Interval.h:51
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
std::vector< Treal > threshVal
Definition Perturbation.h:84
Treal errorTol
Definition Perturbation.h:78
int nIter
Definition Perturbation.h:83
void checkIdempotencies(std::vector< Treal > &idemErrors)
Definition Perturbation.h:256
mat::Interval< Treal > const & allEigs
Definition Perturbation.h:76
Treal deltaMax
Definition Perturbation.h:77
mat::Interval< Treal > gap
Definition Perturbation.h:75
std::vector< Treal > sigma
Definition Perturbation.h:85
void run()
Definition Perturbation.h:188
std::vector< Tmatrix * > & X
Definition Perturbation.h:74
mat::normType const norm
Definition Perturbation.h:79
void checkMaxSubspaceError(Treal &subsError)
Definition Perturbation.h:303
void dryRun()
Dry run to obtain some needed numbers.
Definition Perturbation.h:138
std::vector< Tmatrix * > const & F
Definition Perturbation.h:73
void perturb()
Definition Perturbation.h:60
void checkCommutators(std::vector< Treal > &commErrors, TmatNoSymm const &dummyMat)
Definition Perturbation.h:275
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:105
Tvector & vect
Definition Perturbation.h:80
normType
Definition matInclude.h:139
Definition Perturbation.h:40
Treal template_blas_fabs(Treal x)