ergo
Perturbation.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
38#ifndef MAT_PERTURBATION
39#define MAT_PERTURBATION
40namespace per {
41 template<typename Treal, typename Tmatrix, typename Tvector>
43 public:
44 Perturbation(std::vector<Tmatrix *> const & F,
46 std::vector<Tmatrix *> & D,
55 Treal const deltaMax,
56 Treal const errorTol,
57 mat::normType const norm,
58 Tvector & vect
59 );
60 void perturb() {
61 dryRun();
62 run();
63 }
64
65 void checkIdempotencies(std::vector<Treal> & idemErrors);
66 template<typename TmatNoSymm>
67 void checkCommutators(std::vector<Treal> & commErrors,
68 TmatNoSymm const & dummyMat);
69 void checkMaxSubspaceError(Treal & subsError);
70
71 protected:
72 /* This is input from the beginning */
73 std::vector<Tmatrix *> const & F;
74 std::vector<Tmatrix *> & X;
77 Treal deltaMax;
78 Treal errorTol;
80 Tvector & vect;
81
82 /* These variables are set in the dry run. */
83 int nIter;
84 std::vector<Treal> threshVal;
85 std::vector<Treal> sigma;
86
97 void dryRun();
98 void run();
99 private:
100
101 };
102
103 template<typename Treal, typename Tmatrix, typename Tvector>
105 Perturbation(std::vector<Tmatrix *> const & F_in,
106 std::vector<Tmatrix *> & X_in,
107 mat::Interval<Treal> const & gap_in,
108 mat::Interval<Treal> const & allEigs_in,
109 Treal const deltaMax_in,
110 Treal const errorTol_in,
111 mat::normType const norm_in,
112 Tvector & vect_in)
113 : F(F_in), X(X_in), gap(gap_in), allEigs(allEigs_in),
114 deltaMax(deltaMax_in), errorTol(errorTol_in), norm(norm_in),
115 vect(vect_in) {
116 if (!X.empty())
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]));
120
121 Treal lmin = allEigs.low();
122 Treal lmax = allEigs.upp();
123
124 /***** Initial linear transformation of matrix sequence. */
125 typename std::vector<Tmatrix *>::iterator matIt = X.begin();
126 /* Scale to [0, 1] interval and negate */
127 (*matIt)->add_identity(-lmax);
128 *(*matIt) *= ((Treal)1.0 / (lmin - lmax));
129 matIt++;
130 /* ...and derivatives: */
131 for ( ; matIt != X.end(); matIt++ )
132 *(*matIt) *= ((Treal)-1.0 / (lmin - lmax));
133 /* Compute transformed gap */
134 gap = (gap - lmax) / (lmin - lmax);
135 }
136
137 template<typename Treal, typename Tmatrix, typename Tvector>
139 Treal errorTolPerIter;
140 int nIterGuess = 0;
141 nIter = 1;
142 Treal lumo;
143 Treal homo;
144 Treal m;
145 Treal g;
146 while (nIterGuess < nIter) {
147 nIterGuess++;
148 errorTolPerIter = 0.5 * errorTol /nIterGuess;
149 nIter = 0;
150 mat::Interval<Treal> gapTmp(gap);
151 sigma.resize(0);
152 threshVal.resize(0);
153 while (gapTmp.low() > 0.5 * errorTol || gapTmp.upp() < 0.5 * errorTol) {
154 lumo = gapTmp.low();
155 homo = gapTmp.upp();
156 m = gapTmp.midPoint();
157 g = gapTmp.length();
158 if (m > 0.5) {
159 lumo = lumo*lumo;
160 homo = homo*homo;
161 sigma.push_back(-1);
162 }
163 else {
164 lumo = 2*lumo - lumo*lumo;
165 homo = 2*homo - homo*homo;
166 sigma.push_back(1);
167 }
168 /* Compute threshold value necessary to converge. */
169 Treal forceConvThresh = template_blas_fabs(1-2*m) * g / 10;
170 /* We divide by 10 > 2 so that this loop converges at some point. */
171 /* Compute threshold value necessary to maintain accuracy in subspace.*/
172 Treal subspaceThresh = errorTolPerIter * (homo-lumo) / (1+errorTolPerIter);
173 /* Choose the most restrictive threshold of the two. */
174 threshVal.push_back(forceConvThresh < subspaceThresh ?
175 forceConvThresh : subspaceThresh);
176 homo -= threshVal.back();
177 lumo += threshVal.back();
178 gapTmp = mat::Interval<Treal>(lumo, homo);
179 if (gapTmp.empty())
180 throw "Perturbation<Treal, Tmatrix, Tvector>::dryRun() : Perturbation iterations will fail to converge; Gap is too small or desired accuracy too high.";
181 nIter++;
182 } /* end 2nd while */
183 } /* end 1st while */
184 /* Now, we have nIter, threshVal, and sigma. */
185 }
186
187 template<typename Treal, typename Tmatrix, typename Tvector>
189 Treal const ONE = 1.0;
190 mat::SizesAndBlocks rowsCopy;
191 X.front()->getRows(rowsCopy);
192 mat::SizesAndBlocks colsCopy;
193 X.front()->getCols(colsCopy);
194 Tmatrix tmpMat;
195 // tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy);
196 int nMatrices;
197 Treal threshValPerOrder;
198 Treal chosenThresh;
199 for (int iter = 0; iter < nIter; iter++) {
200 std::cout<<"\n\nInside outer loop iter = "<<iter
201 <<" nIter = "<<nIter
202 <<" sigma = "<<sigma[iter]<<std::endl;
203 /* Number of matrices increases by 1 in each iteration: */
204 X.push_back(new Tmatrix);
205 nMatrices = X.size();
206 X[nMatrices-1]->resetSizesAndBlocks(rowsCopy, colsCopy);
207 /* Compute threshold value for each order. */
208 threshValPerOrder = threshVal[iter] / nMatrices;
209 /* Loop through all nonzero orders. */
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++) {
219 /* X[j] = X[j] - sigma * X[k] * X[j-k] */
220 Tmatrix::ssmmUpperTriangleOnly(-sigma[iter], *X[k], *X[j-k],
221 ONE, tmpMat);
222 } /* End 3rd for */
223 std::cout<<"tmpMat.eucl() (after for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl;
224 *X[j] = tmpMat;
225
226 /* Truncate tmpMat, remove if it becomes zero. */
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;
232#if 1
233 /* If the current matrix is zero AND
234 * it is the last matrix
235 */
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()
239 <<std::endl;
240 delete X[j];
241 X.pop_back();
242 }
243 else
244 std::cout<<"NO DELETION: j = "<<j<<" X.size() = "<<X.size()
245 <<" X[j] = "<<X[j]<< " X[j]->frob() = "<<X[j]->frob()
246 <<std::endl;
247#endif
248
249 } /* End 2nd for (Loop through orders) */
250 } /* End 1st for (Loop through iterations) */
251 } /* End run() */
252
253
254 template<typename Treal, typename Tmatrix, typename Tvector>
256 checkIdempotencies(std::vector<Treal> & idemErrors) {
257 Tmatrix tmpMat;
258 Treal const ONE = 1.0;
259 unsigned int j;
260 for (unsigned int m = 0; m < X.size(); ++m) {
261 tmpMat = (-ONE) * (*X[m]);
262 for (unsigned int i = 0; i <= m; ++i) {
263 j = m - i;
264 /* TMP = TMP + X[i] * X[j] */
265 Tmatrix::ssmmUpperTriangleOnly(ONE, *X[i], *X[j], ONE, tmpMat);
266 }
267 /* TMP is symmetric! */
268 idemErrors.push_back(tmpMat.eucl(vect,1e-10));
269 }
270 }
271
272 template<typename Treal, typename Tmatrix, typename Tvector>
273 template<typename TmatNoSymm>
275 checkCommutators(std::vector<Treal> & commErrors,
276 TmatNoSymm const & dummyMat) {
277 mat::SizesAndBlocks rowsCopy;
278 X.front()->getRows(rowsCopy);
279 mat::SizesAndBlocks colsCopy;
280 X.front()->getCols(colsCopy);
281 TmatNoSymm tmpMat;
282 tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy);
283 Treal const ONE = 1.0;
284 unsigned int j;
285 for (unsigned int m = 0; m < X.size(); ++m) {
286 tmpMat = 0;
287 std::cout<<"New loop\n";
288 for (unsigned int i = 0; i <= m && i < F.size(); ++i) {
289 j = m - i;
290 std::cout<<i<<", "<<j<<std::endl;
291 /* TMP = TMP + F[i] * X[j] - X[j] * F[i] */
292 tmpMat += ONE * (*F[i]) * (*X[j]);
293 tmpMat += -ONE * (*X[j]) * (*F[i]);
294 }
295 /* TMP is not symmetric! */
296 commErrors.push_back(tmpMat.frob());
297 }
298 }
299
300
301 template<typename Treal, typename Tmatrix, typename Tvector>
303 checkMaxSubspaceError(Treal & subsError) {
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]);
308 /***** Initial linear transformation of matrix. */
309 Treal lmin = allEigs.low();
310 Treal lmax = allEigs.upp();
311 /* Scale to [0, 1] interval and negate */
312 XdeltaMax.add_identity(-lmax);
313 XdeltaMax *= ((Treal)1.0 / (lmin - lmax));
314
315 Tmatrix X2;
316 for (int iter = 0; iter < nIter; iter++) {
317 X2 = ONE * XdeltaMax * XdeltaMax;
318 if (sigma[iter] == Treal(1.0)) {
319 XdeltaMax *= 2.0;
320 XdeltaMax -= X2;
321 }
322 else {
323 XdeltaMax = X2;
324 }
325 } /* End of for (Loop through iterations) */
326
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);
332 }
333
334
335
336} /* end namespace mat */
337#endif
338
Definition Interval.h:46
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
Vector class.
Definition Vector.h:54
Definition Perturbation.h:42
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)