lanczos_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file lanczos_base.h
24  \brief File defining \ref o2scl_linalg::lanczos
25 */
26 
27 #ifdef DOXYGEN
28 namespace o2scl_linalg {
29 #endif
30 
31  /** \brief Lanczos diagonalization
32 
33  This class approximates the largest eigenvalues of a symmetric
34  matrix.
35 
36  The vector and matrix types can be any type which provides access
37  via \c operator[] or \c operator(), given a suitable vector
38  allocation type.
39 
40  The tridiagonalization routine was rewritten from the EISPACK
41  routines \c imtql1.f (but uses \c gsl_hypot() instead of \c
42  pythag.f).
43 
44  \future The function eigen_tdiag() automatically sorts the
45  eigenvalues, which may not be necessary.
46 
47  \future Do something better than the simple matrix-vector product.
48  For example, use dgemm() and allow user to specify column or
49  row-major.
50 
51  \future Rework memory allocation to perform as needed.
52 
53  \comment
54  We need the o2scl:: prefix in the default template parameters
55  below because this class isn't in the o2scl namespace.
56  \endcomment
57  */
58  template<class vec_t=boost::numeric::ublas::vector<double>,
59  class mat_t=boost::numeric::ublas::matrix<double> > class lanczos {
60 
61  public:
62 
63  lanczos() {
64  td_iter=30;
65  td_lasteval=0;
66  }
67 
68  /** \brief Number of iterations for finding the eigenvalues of the
69  tridiagonal matrix (default 30)
70  */
71  size_t td_iter;
72 
73  /** \brief The index for the last eigenvalue not determined if
74  tridiagonalization fails
75  */
76  size_t td_lasteval;
77 
78  /** \brief Approximate the largest eigenvalues of a symmetric
79  matrix \c mat using the Lanczos method
80 
81  Given a square matrix \c mat with size \c size, this function
82  applies \c n_iter iterations of the Lanczos algorithm to
83  produce \c n_iter approximate eigenvalues stored in \c
84  eigen. As a by-product, this function also partially
85  tridiagonalizes the matrix placing the result in \c diag and
86  \c off_diag. Before calling this function, space must have
87  already been allocated for \c eigen, \c diag, and \c
88  off_diag. All three of these arrays must have at least enough
89  space for \c n_iter elements.
90 
91  Choosing /c n_iter = \c size will produce all of the exact
92  eigenvalues and the corresponding tridiagonal matrix, but this
93  may be slower than diagonalizing the matrix directly.
94  */
95  int eigenvalues(size_t size, mat_t &mat, size_t n_iter,
96  vec_t &eigen, vec_t &diag, vec_t &off_diag) {
97  double t;
98  bool cont=true;
99  size_t i, j, k;
100 
101  vec_t v;
102  vec_t w;
103  vec_t b3;
104  vec_t prod;
105 
106  v.resize(size);
107  w.resize(size);
108  b3.resize(size);
109  prod.resize(size);
110 
111  // Pick a unit vector
112  O2SCL_IX(w,0)=1.0;
113  for(i=1;i<size;i++) O2SCL_IX(w,i)=0.0;
114 
115  for(i=0;i<size;i++) O2SCL_IX(v,i)=0;
116  j=0;
117  while (cont) {
118  if (j!=0) {
119  for(i=0;i<size;i++) {
120  t=O2SCL_IX(w,i);
121  O2SCL_IX(w,i)=O2SCL_IX(v,i)/O2SCL_IX(off_diag,j-1);
122  O2SCL_IX(v,i)=-O2SCL_IX(off_diag,j-1)*t;
123  }
124  }
125  product(size,mat,w,prod);
126  for(k=0;k<size;k++) O2SCL_IX(v,k)+=O2SCL_IX(prod,k);
127  O2SCL_IX(diag,j)=0.0;
128  O2SCL_IX(off_diag,j)=0.0;
129 
130  for(k=0;k<size;k++) O2SCL_IX(diag,j)+=O2SCL_IX(w,k)*O2SCL_IX(v,k);
131  for(k=0;k<size;k++) O2SCL_IX(v,k)-=O2SCL_IX(diag,j)*O2SCL_IX(w,k);
132  for(k=0;k<size;k++) O2SCL_IX(off_diag,j)+=O2SCL_IX(v,k)*O2SCL_IX(v,k);
133  O2SCL_IX(off_diag,j)=sqrt(O2SCL_IX(off_diag,j));
134  j++;
135 
136  if (j>=n_iter || O2SCL_IX(off_diag,j-1)==0.0) cont=false;
137 
138  if (j>0) {
139  for(k=0;k<size-1;k++) {
140  O2SCL_IX(b3,k+1)=O2SCL_IX(off_diag,k);
141  O2SCL_IX(eigen,k)=O2SCL_IX(diag,k);
142  }
143  O2SCL_IX(eigen,size-1)=O2SCL_IX(diag,size-1);
144  if (eigen_tdiag(j,eigen,b3)!=0) {
145  O2SCL_ERR2("Call to eigen_tdiag() in lanczos::",
146  "eigenvalues() failed.",o2scl::exc_efailed);
147  }
148  }
149  }
150 
151  return 0;
152  }
153 
154  /** \brief In-place diagonalization of a tridiagonal matrix
155 
156  On input, the vectors \c diag and \c off_diag should both be
157  vectors of size \c n. The diagonal entries stored in \c diag,
158  and the \f$ n-1 \f$ off-diagonal entries should be stored in
159  \c off_diag, starting with \c off_diag[1]. The value in \c
160  off_diag[0] is unused. The vector \c off_diag is destroyed by
161  the computation.
162 
163  This uses an implict QL method from the EISPACK routine \c
164  imtql1. The value of \c ierr from the original Fortran routine
165  is stored in \ref td_lasteval.
166 
167  */
168  int eigen_tdiag(size_t n, vec_t &diag, vec_t &off_diag) {
169 
170  // The variable 'i' is set to zero here because Cygwin complained
171  // about uninitialized variables. This is probably ok, but it
172  // would be nice to double check that there isn't a problem with
173  // setting i=0 here.
174 
175  int i=0,j,l,m,mml;
176  double b,c,f,g,p,r,s,tst1,tst2;
177 
178  if (n==1) return 0;
179 
180  for(size_t ij=1;ij<n;ij++) {
181  O2SCL_IX(off_diag,ij-1)=O2SCL_IX(off_diag,ij);
182  }
183  O2SCL_IX(off_diag,n-1)=0.0;
184 
185  bool done=false;
186 
187  l=1;
188  j=0;
189 
190  while (done==false && l<=((int)n)) {
191 
192  // Look for small sub-diagonal element
193  bool idone=false;
194  for(m=l;m<((int)n) && idone==false;m++) {
195  tst1=fabs(O2SCL_IX(diag,m-1))+fabs(O2SCL_IX(diag,m));
196  tst2=tst1+fabs(O2SCL_IX(off_diag,m-1));
197  if (tst2==tst1) {
198  m--;
199  idone=true;
200  }
201  }
202 
203  p=O2SCL_IX(diag,l-1);
204 
205  if (m!=l && j==((int)td_iter)) {
206 
207  // Set error. No convergence after td_iter iterations
208  td_lasteval=l;
209  O2SCL_ERR("No convergence in lanczos::eigen_tdiag()",
211  }
212 
213  if (m!=l) {
214 
215  j++;
216 
217  // Form shift
218  g=(O2SCL_IX(diag,l)-p)/(2.0*O2SCL_IX(off_diag,l-1));
219  r=gsl_hypot(g,1.0);
220 
221  g=O2SCL_IX(diag,m-1)-p+O2SCL_IX(off_diag,l-1)/
222  (g+(g>=0.0 ? fabs(r) : -fabs(r)));
223  s=1.0;
224  c=1.0;
225  p=0.0;
226  mml=m-l;
227 
228  for(int ii=1;ii<=mml;ii++) {
229 
230  i=m-ii;
231  f=s*O2SCL_IX(off_diag,i-1);
232  b=c*O2SCL_IX(off_diag,i-1);
233  r=gsl_hypot(f,g);
234  O2SCL_IX(off_diag,i)=r;
235 
236  if (r==0.0) {
237 
238  // Recover from underflow
239  O2SCL_IX(diag,i)-=p;
240  O2SCL_IX(off_diag,m-1)=0.0;
241  ii=mml+1;
242 
243  } else {
244 
245  s=f/r;
246  c=g/r;
247  g=O2SCL_IX(diag,i)-p;
248  r=(O2SCL_IX(diag,i-1)-g)*s+2.0*c*b;
249  p=s*r;
250  O2SCL_IX(diag,i)=g+p;
251  g=c*r-b;
252 
253  }
254 
255  }
256 
257  O2SCL_IX(diag,l-1)-=p;
258  O2SCL_IX(off_diag,l-1)=g;
259  O2SCL_IX(off_diag,m-1)=0.0;
260 
261 
262  } else {
263 
264  // Order eigenvalues
265 
266  if (l==1) {
267 
268  i=1;
269  O2SCL_IX(diag,i-1)=p;
270 
271  } else {
272 
273  bool skip=false;
274  for(int ii=2;ii<=l;ii++) {
275  i=l+2-ii;
276  if (p>=O2SCL_IX(diag,i-2)) {
277  ii=l+1;
278  skip=true;
279  } else {
280  O2SCL_IX(diag,i-1)=O2SCL_IX(diag,i-2);
281  }
282  }
283 
284  if (skip==false) i=1;
285  O2SCL_IX(diag,i-1)=p;
286  }
287 
288  j=0;
289  l++;
290  }
291 
292  }
293 
294  return 0;
295  }
296 
297 #ifndef DOXYGEN_INTERNAL
298 
299  protected:
300 
301  /** \brief Simple matrix-vector product
302 
303  It is assumed that memory is already allocated for \c prod.
304  */
305  void product(size_t n, mat_t &a, vec_t &w, vec_t &prod) {
306  size_t i, j;
307  for(i=0;i<n;i++) {
308  O2SCL_IX(prod,i)=0.0;
309  for(j=0;j<n;j++) {
310  O2SCL_IX(prod,i)+=O2SCL_IX2(a,i,j)*O2SCL_IX(w,j);
311  }
312  }
313  return;
314  }
315 
316 #endif
317 
318  };
319 
320 #ifdef DOXYGEN
321 }
322 #endif
int eigen_tdiag(size_t n, vec_t &diag, vec_t &off_diag)
In-place diagonalization of a tridiagonal matrix.
Definition: lanczos_base.h:168
generic failure
Definition: err_hnd.h:61
int eigenvalues(size_t size, mat_t &mat, size_t n_iter, vec_t &eigen, vec_t &diag, vec_t &off_diag)
Approximate the largest eigenvalues of a symmetric matrix mat using the Lanczos method.
Definition: lanczos_base.h:95
size_t td_lasteval
The index for the last eigenvalue not determined if tridiagonalization fails.
Definition: lanczos_base.h:76
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
Lanczos diagonalization.
Definition: lanczos_base.h:59
size_t td_iter
Number of iterations for finding the eigenvalues of the tridiagonal matrix (default 30) ...
Definition: lanczos_base.h:71
void product(size_t n, mat_t &a, vec_t &w, vec_t &prod)
Simple matrix-vector product.
Definition: lanczos_base.h:305

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).