ergo
MatrixTridiagSymmetric.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_MATRIXTRIDIAGSYMMETRIC
39#define MAT_MATRIXTRIDIAGSYMMETRIC
40#include "mat_gblas.h"
41namespace mat { /* Matrix namespace */
42 namespace arn { /* Arnoldi type methods namespace */
46 template<typename Treal>
48 public:
49 explicit MatrixTridiagSymmetric(int k = 100)
51 size(0), capacity(k) {}
52 void increase(Treal const & alpha, Treal const & beta) {
53 if (size + 1 > capacity)
55 ++size;
56 alphaVec[size - 1] = alpha;
57 betaVec[size - 1] = beta;
58 }
60 delete[] alphaVec;
61 delete[] betaVec;
62 }
63
64
65 void update_beta(Treal const & beta)
66 {
67 betaVec[size-1] = beta;
68 }
69
72 Treal* acc,
73 int & nEigsFound,
74 Treal const lowBound,
75 Treal const uppBound,
76 Treal const abstol = 0) const;
79 Treal* acc,
80 int const lowInd,
81 int const uppInd,
82 Treal const abstol = 0) const;
83 inline void clear() {
84 size = 0;
85 }
86 protected:
89 int size;
91 void increaseCapacity(int const newCapacity);
92 private:
93 };
94
95 template<typename Treal>
97 getEigsByInterval(Treal* eigVals, /* length: >= nEigsFound */
98 Treal* eigVectors, /* length: >= size * nEigsFound */
99 Treal* acc, /* length: size */
100 int & nEigsFound, /* The number of found eigenpairs. */
101 Treal const lowBound,
102 Treal const uppBound,
103 Treal const absTol) const {
104 Treal* eigArray = new Treal[size];
105 Treal* alphaCopy = new Treal[size];
106 Treal* betaCopy = new Treal[size];
107 Treal* work = new Treal[5 * size];
108 int* iwork = new int[5 * size];
109 int* ifail = new int[size];
110 for (int ind = 0; ind < size; ind++){
111 alphaCopy[ind] = alphaVec[ind];
112 betaCopy[ind] = betaVec[ind];
113 }
114 int dummy = -1;
115 int info;
116 /* Find eigenvalues */
117 /* FIXME: change to stevr */
118 mat::stevx("V", "V", &size, alphaCopy, betaCopy,
119 &lowBound, &uppBound, &dummy, &dummy,
120 &absTol,
121 &nEigsFound, eigArray, eigVectors, &size,
122 work, iwork, ifail,
123 &info);
124 assert(info == 0);
125 for (int ind = 0; ind < nEigsFound; ind++) {
127 acc[ind] = betaCopy[size - 1] *
128 template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9; // FIXME: WHY IS THERE A FACTOR 0.9 HERE ?!?
129 }
130 delete[] eigArray;
131 delete[] alphaCopy;
132 delete[] betaCopy;
133 delete[] work;
134 delete[] iwork;
135 delete[] ifail;
136 }
137
138 template<typename Treal>
140 getEigsByIndex(Treal* eigVals, /* length: uppInd-lowInd+1 */
141 Treal* eigVectors, /* length: size*(uppInd-lowInd+1) */
142 Treal* acc, /* length: size */
143 int const lowInd,
144 int const uppInd,
145 Treal const abstol) const {
146 Treal* eigArray = new Treal[size];
147 Treal* alphaCopy = new Treal[size];
148 Treal* betaCopy = new Treal[size];
149 for (int ind = 0; ind < size; ind++){
150 alphaCopy[ind] = alphaVec[ind];
151 betaCopy[ind] = betaVec[ind];
152 }
153#if 1
154 // Emanuel note 2010-03-14:
155 // The following code uses stevr instead of stevx for two reasons:
156 // 1) Due to a bug in LAPACK we previously computed all
157 // eigenvalues (see Elias' note below) which turned out to be
158 // too time consuming in some cases.
159 // 2) Contrary to stevx, stevr should never fail to compute the
160 // desired eigenpairs unless there is a bug in the implementation
161 // or erroneous input.
162 int const lowIndNew(lowInd + 1);
163 int const uppIndNew(uppInd + 1);
164 int nEigsWanted = uppInd - lowInd + 1;
165 int nEigsFound = 0;
166 int* isuppz = new int[2*nEigsWanted];
167 Treal* work;
168 int lwork = -1;
169 int* iwork;
170 int liwork = -1;
171 Treal dummy = -1.0;
172 int info;
173 // First do a workspace query:
175 int iwork_query;
176 mat::stevr("V", "I", &size, alphaCopy, betaCopy,
177 &dummy, &dummy, &lowIndNew, &uppIndNew,
178 &abstol,
179 &nEigsFound, eigArray, eigVectors, &size,
180 isuppz,
181 &work_query, &lwork, &iwork_query, &liwork, &info);
184 work = new Treal[lwork];
185 iwork = new int[liwork];
186 mat::stevr("V", "I", &size, alphaCopy, betaCopy,
187 &dummy, &dummy, &lowIndNew, &uppIndNew,
188 &abstol,
189 &nEigsFound, eigArray, eigVectors, &size,
190 isuppz,
191 work, &lwork, iwork, &liwork, &info);
192 if (info)
193 std::cout << "info = " << info <<std::endl;
194 assert(info == 0);
195 assert(nEigsFound == nEigsWanted);
196 for (int ind = 0; ind < nEigsFound; ind++) {
198 acc[ind] = betaCopy[size - 1] *
199 template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9; // FIXME: WHY IS THERE A FACTOR 0.9 HERE ?!?
200 }
201 delete[] eigArray;
202 delete[] alphaCopy;
203 delete[] betaCopy;
204 delete[] isuppz;
205 delete[] work;
206 delete[] iwork;
207
208#else
209 Treal* work = new Treal[5 * size];
210 int* iwork = new int[5 * size];
211 int* ifail = new int[size];
212 Treal dummy = -1.0;
213 int info;
214 int nEigsFound = 0;
215 /*
216 Elias note 2007-07-02:
217 There have been some problems with stevx returning with info=0
218 at the same time as nEigsFound != uppInd - lowInd + 1.
219 This is due to a bug in LAPACK which has been reported and confirmed.
220 To avoid this problem Elias changed the code so that ALL eigenvectors
221 are computed and then the desired ones are selected from the
222 complete list.
223 */
224#if 1
225 /* Original version of code calling stevx to get only the
226 desired eigenvalues/vectors. */
227 int const lowIndNew(lowInd + 1);
228 int const uppIndNew(uppInd + 1);
229 mat::stevx("V", "I", &size, alphaCopy, betaCopy,
230 &dummy, &dummy, &lowIndNew, &uppIndNew,
231 &abstol,
232 &nEigsFound, eigArray, eigVectors, &size,
233 work, iwork, ifail,
234 &info);
235 assert(info == 0);
236 assert(nEigsFound == uppInd - lowInd + 1);
237 for (int ind = 0; ind < nEigsFound; ind++) {
239 acc[ind] = betaCopy[size - 1] *
240 template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9; // FIXME: WHY IS THERE A FACTOR 0.9 HERE ?!?
241 }
242#else
243 /* Modified version of code calling stevx to get ALL
244 eigenvalues/vectors, and then picking out the desired ones. */
245 Treal* eigVectorsTmp = new Treal[size*size];
246 int const lowIndNew(1);
247 int const uppIndNew(size);
248 mat::stevx("V", "A", &size, alphaCopy, betaCopy,
249 &dummy, &dummy, &lowIndNew, &uppIndNew,
250 &abstol,
252 work, iwork, ifail,
253 &info);
254 assert(info == 0);
255 assert(nEigsFound == size);
256 int nEigsWanted = uppInd - lowInd + 1;
257 /* Copy desired eigenvectors from eigVectorsTmp to eigVectors. */
258 for(int i = 0; i < nEigsWanted; i++)
259 for(int j = 0; j < size; j++)
260 eigVectors[i*size+j] = eigVectorsTmp[(lowInd+i)*size+j];
261 delete [] eigVectorsTmp;
262 for (int ind = 0; ind < nEigsWanted; ind++) {
264 acc[ind] = betaCopy[size - 1] *
265 template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9; // FIXME: WHY IS THERE A FACTOR 0.9 HERE ?!?
266 }
267#endif
268 delete[] eigArray;
269 delete[] alphaCopy;
270 delete[] betaCopy;
271 delete[] work;
272 delete[] iwork;
273 delete[] ifail;
274#endif
275 }
276
277
278
279 template<typename Treal>
282 capacity = newCapacity;
283 Treal* alphaNew = new Treal[capacity];
284 Treal* betaNew = new Treal[capacity];
285 for (int ind = 0; ind < size; ind++){
286 alphaNew[ind] = alphaVec[ind];
287 betaNew[ind] = betaVec[ind];
288 }
289 delete[] alphaVec;
290 delete[] betaVec;
291 alphaVec = alphaNew;
292 betaVec = betaNew;
293 }
294
295
296 } /* end namespace arn */
297} /* end namespace mat */
298#endif
Tridiagonal symmetric matrix class template.
Definition MatrixTridiagSymmetric.h:47
void clear()
Definition MatrixTridiagSymmetric.h:83
int size
Definition MatrixTridiagSymmetric.h:89
MatrixTridiagSymmetric(int k=100)
Definition MatrixTridiagSymmetric.h:49
void increase(Treal const &alpha, Treal const &beta)
Definition MatrixTridiagSymmetric.h:52
void getEigsByIndex(Treal *eigVals, Treal *eigVectors, Treal *acc, int const lowInd, int const uppInd, Treal const abstol=0) const
Definition MatrixTridiagSymmetric.h:140
void getEigsByInterval(Treal *eigVals, Treal *eigVectors, Treal *acc, int &nEigsFound, Treal const lowBound, Treal const uppBound, Treal const abstol=0) const
Definition MatrixTridiagSymmetric.h:97
void increaseCapacity(int const newCapacity)
Definition MatrixTridiagSymmetric.h:281
Treal * betaVec
Definition MatrixTridiagSymmetric.h:88
virtual ~MatrixTridiagSymmetric()
Definition MatrixTridiagSymmetric.h:59
Treal * alphaVec
Definition MatrixTridiagSymmetric.h:87
int capacity
Definition MatrixTridiagSymmetric.h:90
void update_beta(Treal const &beta)
Definition MatrixTridiagSymmetric.h:65
C++ interface to a subset of BLAS and LAPACK.
Definition allocate.cc:39
static void stevr(const char *jobz, const char *range, const int *n, T *d, T *e, const T *vl, const T *vu, const int *il, const int *iu, const T *abstol, int *m, T *w, T *z, const int *ldz, int *isuppz, T *work, int *lwork, int *iwork, int *liwork, int *info)
Definition mat_gblas.h:369
static void stevx(const char *jobz, const char *range, const int *n, T *d, T *e, const T *vl, const T *vu, const int *il, const int *iu, const T *abstol, int *m, T *w, T *z, const int *ldz, T *work, int *iwork, int *ifail, int *info)
Definition mat_gblas.h:358
static Treal getMachineEpsilon()
Definition matInclude.h:147
Treal template_blas_fabs(Treal x)