ergo
MatrixTridiagSymmetric.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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"
41 namespace mat { /* Matrix namespace */
42  namespace arn { /* Arnoldi type methods namespace */
46  template<typename Treal>
48  public:
49  explicit MatrixTridiagSymmetric(int k = 100)
50  : alphaVec(new Treal[k]), betaVec(new Treal[k]),
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 
70  void getEigsByInterval(Treal* eigVals,
71  Treal* eigVectors,
72  Treal* acc,
73  int & nEigsFound,
74  Treal const lowBound,
75  Treal const uppBound,
76  Treal const abstol = 0) const;
77  void getEigsByIndex(Treal* eigVals,
78  Treal* eigVectors,
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:
87  Treal* alphaVec;
88  Treal* betaVec;
89  int size;
90  int capacity;
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++) {
126  eigVals[ind] = eigArray[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:
174  Treal work_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);
182  lwork = int(work_query);
183  liwork = iwork_query;
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++) {
197  eigVals[ind] = eigArray[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++) {
238  eigVals[ind] = eigArray[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,
251  &nEigsFound, eigArray, eigVectorsTmp, &size,
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++) {
263  eigVals[ind] = eigArray[lowInd+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>
281  increaseCapacity(int const newCapacity) {
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
int capacity
Definition: MatrixTridiagSymmetric.h:90
Treal * alphaVec
Definition: MatrixTridiagSymmetric.h:87
virtual ~MatrixTridiagSymmetric()
Definition: MatrixTridiagSymmetric.h:59
MatrixTridiagSymmetric(int k=100)
Definition: MatrixTridiagSymmetric.h:49
int size
Definition: MatrixTridiagSymmetric.h:89
Tridiagonal symmetric matrix class template.
Definition: MatrixTridiagSymmetric.h:47
void clear()
Definition: MatrixTridiagSymmetric.h:83
Definition: allocate.cc:39
Treal template_blas_fabs(Treal x)
C++ interface to a subset of BLAS and LAPACK.
void update_beta(Treal const &beta)
Definition: MatrixTridiagSymmetric.h:65
void getEigsByIndex(Treal *eigVals, Treal *eigVectors, Treal *acc, int const lowInd, int const uppInd, Treal const abstol=0) const
Definition: MatrixTridiagSymmetric.h:140
Treal * betaVec
Definition: MatrixTridiagSymmetric.h:88
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
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
void increaseCapacity(int const newCapacity)
Definition: MatrixTridiagSymmetric.h:281
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
void increase(Treal const &alpha, Treal const &beta)
Definition: MatrixTridiagSymmetric.h:52