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,
113 : F(F_in), X(X_in), gap(gap_in), allEigs(allEigs_in),
114 deltaMax(deltaMax_in), errorTol(errorTol_in), norm(norm_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;
153 while (gapTmp.
low() > 0.5 * errorTol || gapTmp.
upp() < 0.5 * errorTol) {
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);
176 homo -= threshVal.back();
177 lumo += threshVal.back();
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]);
309 Treal lmin = allEigs.low();
310 Treal lmax = allEigs.upp();
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,
331 vect, errorTol *1e-2);
Treal upp() const
Definition: Interval.h:145
Tvector & vect
Definition: Perturbation.h:80
void dryRun()
Dry run to obtain some needed numbers.
Definition: Perturbation.h:138
Treal low() const
Definition: Interval.h:144
std::vector< Tmatrix * > const & F
Definition: Perturbation.h:73
Vector class.
Definition: Matrix.h:78
mat::Interval< Treal > const & allEigs
Definition: Perturbation.h:76
void run()
Definition: Perturbation.h:188
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
Treal template_blas_fabs(Treal x)
Definition: Interval.h:46
void checkCommutators(std::vector< Treal > &commErrors, TmatNoSymm const &dummyMat)
Definition: Perturbation.h:275
Treal errorTol
Definition: Perturbation.h:78
void perturb()
Definition: Perturbation.h:60
std::vector< Tmatrix * > & X
Definition: Perturbation.h:74
void checkMaxSubspaceError(Treal &subsError)
Definition: Perturbation.h:303
mat::normType const norm
Definition: Perturbation.h:79
Treal midPoint() const
Definition: Interval.h:115
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
Definition: Perturbation.h:42
std::vector< Treal > sigma
Definition: Perturbation.h:85
void checkIdempotencies(std::vector< Treal > &idemErrors)
Definition: Perturbation.h:256
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
Treal deltaMax
Definition: Perturbation.h:77
int nIter
Definition: Perturbation.h:83
bool empty() const
Definition: Interval.h:51
std::vector< Treal > threshVal
Definition: Perturbation.h:84
mat::Interval< Treal > gap
Definition: Perturbation.h:75
Definition: Perturbation.h:40
normType
Definition: matInclude.h:139