ergo
MatrixTridiagSymmetric.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
36 #ifndef MAT_MATRIXTRIDIAGSYMMETRIC
37 #define MAT_MATRIXTRIDIAGSYMMETRIC
38 #include "gblas.h"
39 namespace mat { /* Matrix namespace */
40  namespace arn { /* Arnoldi type methods namespace */
44  template<typename Treal>
46  public:
47  explicit MatrixTridiagSymmetric(int k = 100)
48  : alphaVec(new Treal[k]), betaVec(new Treal[k]),
49  size(0), capacity(k) {}
50  void increase(Treal const & alpha, Treal const & beta) {
51  if (size + 1 > capacity)
53  ++size;
54  alphaVec[size - 1] = alpha;
55  betaVec[size - 1] = beta;
56  }
58  delete[] alphaVec;
59  delete[] betaVec;
60  }
61  void getEigsByInterval(Treal* eigVals,
62  Treal* eigVectors,
63  Treal* acc,
64  int & nEigsFound,
65  Treal const lowBound,
66  Treal const uppBound,
67  Treal const abstol = 0) const;
68  void getEigsByIndex(Treal* eigVals,
69  Treal* eigVectors,
70  Treal* acc,
71  int const lowInd,
72  int const uppInd,
73  Treal const abstol = 0) const;
74  inline void clear() {
75  size = 0;
76  }
77  protected:
78  Treal* alphaVec;
79  Treal* betaVec;
80  int size;
81  int capacity;
82  void increaseCapacity(int const newCapacity);
83  private:
84  };
85 
86  template<typename Treal>
88  getEigsByInterval(Treal* eigVals, /* length: >= nEigsFound */
89  Treal* eigVectors, /* length: >= size * nEigsFound */
90  Treal* acc, /* length: size */
91  int & nEigsFound, /* The number of found eigenpairs. */
92  Treal const lowBound,
93  Treal const uppBound,
94  Treal const absTol) const {
95  Treal* eigArray = new Treal[size];
96  Treal* alphaCopy = new Treal[size];
97  Treal* betaCopy = new Treal[size];
98  Treal* work = new Treal[5 * size];
99  int* iwork = new int[5 * size];
100  int* ifail = new int[size];
101  for (int ind = 0; ind < size; ind++){
102  alphaCopy[ind] = alphaVec[ind];
103  betaCopy[ind] = betaVec[ind];
104  }
105  int dummy = -1;
106  int info;
107  /* Find eigenvalues */
108  /* FIXME: change to stevr */
109  mat::stevx("V", "V", &size, alphaCopy, betaCopy,
110  &lowBound, &uppBound, &dummy, &dummy,
111  &absTol,
112  &nEigsFound, eigArray, eigVectors, &size,
113  work, iwork, ifail,
114  &info);
115  assert(info == 0);
116  for (int ind = 0; ind < nEigsFound; ind++) {
117  eigVals[ind] = eigArray[ind];
118  acc[ind] = betaCopy[size - 1] *
119  template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
120  }
121  delete[] eigArray;
122  delete[] alphaCopy;
123  delete[] betaCopy;
124  delete[] work;
125  delete[] iwork;
126  delete[] ifail;
127  }
128 
129  template<typename Treal>
131  getEigsByIndex(Treal* eigVals, /* length: uppInd-lowInd+1 */
132  Treal* eigVectors, /* length: size*(uppInd-lowInd+1) */
133  Treal* acc, /* length: size */
134  int const lowInd,
135  int const uppInd,
136  Treal const abstol) const {
137  Treal* eigArray = new Treal[size];
138  Treal* alphaCopy = new Treal[size];
139  Treal* betaCopy = new Treal[size];
140  for (int ind = 0; ind < size; ind++){
141  alphaCopy[ind] = alphaVec[ind];
142  betaCopy[ind] = betaVec[ind];
143  }
144 #if 1
145  // Emanuel note 2010-03-14:
146  // The following code uses stevr instead of stevx for two reasons:
147  // 1) Due to a bug in LAPACK we previously computed all
148  // eigenvalues (see Elias' note below) which turned out to be
149  // too time consuming in some cases.
150  // 2) Contrary to stevx, stevr should never fail to compute the
151  // desired eigenpairs unless there is a bug in the implementation
152  // or erroneous input.
153  int const lowIndNew(lowInd + 1);
154  int const uppIndNew(uppInd + 1);
155  int nEigsWanted = uppInd - lowInd + 1;
156  int nEigsFound = 0;
157  int* isuppz = new int[2*nEigsWanted];
158  Treal* work;
159  int lwork = -1;
160  int* iwork;
161  int liwork = -1;
162  Treal dummy = -1.0;
163  int info;
164  // First do a workspace query:
165  Treal work_query;
166  int iwork_query;
167  mat::stevr("V", "I", &size, alphaCopy, betaCopy,
168  &dummy, &dummy, &lowIndNew, &uppIndNew,
169  &abstol,
170  &nEigsFound, eigArray, eigVectors, &size,
171  isuppz,
172  &work_query, &lwork, &iwork_query, &liwork, &info);
173  lwork = int(work_query);
174  liwork = iwork_query;
175  work = new Treal[lwork];
176  iwork = new int[liwork];
177  mat::stevr("V", "I", &size, alphaCopy, betaCopy,
178  &dummy, &dummy, &lowIndNew, &uppIndNew,
179  &abstol,
180  &nEigsFound, eigArray, eigVectors, &size,
181  isuppz,
182  work, &lwork, iwork, &liwork, &info);
183  if (info)
184  std::cout << "info = " << info <<std::endl;
185  assert(info == 0);
186  assert(nEigsFound == nEigsWanted);
187  for (int ind = 0; ind < nEigsFound; ind++) {
188  eigVals[ind] = eigArray[ind];
189  acc[ind] = betaCopy[size - 1] *
190  template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
191  }
192  delete[] eigArray;
193  delete[] alphaCopy;
194  delete[] betaCopy;
195  delete[] isuppz;
196  delete[] work;
197  delete[] iwork;
198 
199 #else
200  Treal* work = new Treal[5 * size];
201  int* iwork = new int[5 * size];
202  int* ifail = new int[size];
203  Treal dummy = -1.0;
204  int info;
205  int nEigsFound = 0;
206  /*
207  Elias note 2007-07-02:
208  There have been some problems with stevx returning with info=0
209  at the same time as nEigsFound != uppInd - lowInd + 1.
210  This is due to a bug in LAPACK which has been reported and confirmed.
211  To avoid this problem Elias changed the code so that ALL eigenvectors
212  are computed and then the desired ones are selected from the
213  complete list.
214  */
215 #if 1
216  /* Original version of code calling stevx to get only the
217  desired eigenvalues/vectors. */
218  int const lowIndNew(lowInd + 1);
219  int const uppIndNew(uppInd + 1);
220  mat::stevx("V", "I", &size, alphaCopy, betaCopy,
221  &dummy, &dummy, &lowIndNew, &uppIndNew,
222  &abstol,
223  &nEigsFound, eigArray, eigVectors, &size,
224  work, iwork, ifail,
225  &info);
226  assert(info == 0);
227  assert(nEigsFound == uppInd - lowInd + 1);
228  for (int ind = 0; ind < nEigsFound; ind++) {
229  eigVals[ind] = eigArray[ind];
230  acc[ind] = betaCopy[size - 1] *
231  template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
232  }
233 #else
234  /* Modified version of code calling stevx to get ALL
235  eigenvalues/vectors, and then picking out the desired ones. */
236  Treal* eigVectorsTmp = new Treal[size*size];
237  int const lowIndNew(1);
238  int const uppIndNew(size);
239  mat::stevx("V", "A", &size, alphaCopy, betaCopy,
240  &dummy, &dummy, &lowIndNew, &uppIndNew,
241  &abstol,
242  &nEigsFound, eigArray, eigVectorsTmp, &size,
243  work, iwork, ifail,
244  &info);
245  assert(info == 0);
246  assert(nEigsFound == size);
247  int nEigsWanted = uppInd - lowInd + 1;
248  /* Copy desired eigenvectors from eigVectorsTmp to eigVectors. */
249  for(int i = 0; i < nEigsWanted; i++)
250  for(int j = 0; j < size; j++)
251  eigVectors[i*size+j] = eigVectorsTmp[(lowInd+i)*size+j];
252  delete [] eigVectorsTmp;
253  for (int ind = 0; ind < nEigsWanted; ind++) {
254  eigVals[ind] = eigArray[lowInd+ind];
255  acc[ind] = betaCopy[size - 1] *
256  template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
257  }
258 #endif
259  delete[] eigArray;
260  delete[] alphaCopy;
261  delete[] betaCopy;
262  delete[] work;
263  delete[] iwork;
264  delete[] ifail;
265 #endif
266  }
267 
268 
269 
270  template<typename Treal>
272  increaseCapacity(int const newCapacity) {
273  capacity = newCapacity;
274  Treal* alphaNew = new Treal[capacity];
275  Treal* betaNew = new Treal[capacity];
276  for (int ind = 0; ind < size; ind++){
277  alphaNew[ind] = alphaVec[ind];
278  betaNew[ind] = betaVec[ind];
279  }
280  delete[] alphaVec;
281  delete[] betaVec;
282  alphaVec = alphaNew;
283  betaVec = betaNew;
284  }
285 
286 
287  } /* end namespace arn */
288 } /* end namespace mat */
289 #endif
void getEigsByInterval(Treal *eigVals, Treal *eigVectors, Treal *acc, int &nEigsFound, Treal const lowBound, Treal const uppBound, Treal const abstol=0) const
Definition: MatrixTridiagSymmetric.h:88
int capacity
Definition: MatrixTridiagSymmetric.h:81
Treal * alphaVec
Definition: MatrixTridiagSymmetric.h:78
virtual ~MatrixTridiagSymmetric()
Definition: MatrixTridiagSymmetric.h:57
MatrixTridiagSymmetric(int k=100)
Definition: MatrixTridiagSymmetric.h:47
int size
Definition: MatrixTridiagSymmetric.h:80
Tridiagonal symmetric matrix class template.
Definition: MatrixTridiagSymmetric.h:45
void clear()
Definition: MatrixTridiagSymmetric.h:74
void getEigsByIndex(Treal *eigVals, Treal *eigVectors, Treal *acc, int const lowInd, int const uppInd, Treal const abstol=0) const
Definition: MatrixTridiagSymmetric.h:131
Treal template_blas_fabs(Treal x)
C++ interface to a subset of BLAS and LAPACK.
Treal * betaVec
Definition: MatrixTridiagSymmetric.h:79
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: gblas.h:367
void increaseCapacity(int const newCapacity)
Definition: MatrixTridiagSymmetric.h:272
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: gblas.h:356
void increase(Treal const &alpha, Treal const &beta)
Definition: MatrixTridiagSymmetric.h:50