ViennaCL - The Vienna Computing Library  1.5.1
lanczos.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_LANCZOS_HPP_
2 #define VIENNACL_LINALG_LANCZOS_HPP_
3 
4 /* =========================================================================
5  Copyright (c) 2010-2014, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the PDF manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
27 #include <cmath>
28 #include <vector>
29 #include "viennacl/vector.hpp"
31 #include "viennacl/linalg/prod.hpp"
36 #include <boost/random.hpp>
37 #include <boost/random/mersenne_twister.hpp>
38 #include <boost/numeric/ublas/matrix.hpp>
39 #include <boost/numeric/ublas/matrix_proxy.hpp>
40 #include <boost/numeric/ublas/matrix_expression.hpp>
41 #include <boost/numeric/ublas/matrix_sparse.hpp>
42 #include <boost/numeric/ublas/vector.hpp>
43 #include <boost/numeric/ublas/operation.hpp>
44 #include <boost/numeric/ublas/vector_expression.hpp>
45 #include <boost/numeric/ublas/io.hpp>
46 
47 namespace viennacl
48 {
49  namespace linalg
50  {
51 
55  {
56  public:
57 
58  enum
59  {
63  };
64 
73  lanczos_tag(double factor = 0.75,
74  vcl_size_t numeig = 10,
75  int met = 0,
76  vcl_size_t krylov = 100) : factor_(factor), num_eigenvalues_(numeig), method_(met), krylov_size_(krylov) {}
77 
79  void num_eigenvalues(int numeig){ num_eigenvalues_ = numeig; }
80 
82  vcl_size_t num_eigenvalues() const { return num_eigenvalues_; }
83 
85  void factor(double fct) { factor_ = fct; }
86 
88  double factor() const { return factor_; }
89 
91  void krylov_size(int max) { krylov_size_ = max; }
92 
94  vcl_size_t krylov_size() const { return krylov_size_; }
95 
97  void method(int met){ method_ = met; }
98 
100  int method() const { return method_; }
101 
102 
103  private:
104  double factor_;
105  vcl_size_t num_eigenvalues_;
106  int method_; // see enum defined above for possible values
107  vcl_size_t krylov_size_;
108 
109  };
110 
111 
112  namespace detail
113  {
124  template< typename MatrixT, typename VectorT >
125  std::vector<
127  >
128  lanczosPRO (MatrixT const& A, VectorT & r, vcl_size_t size, lanczos_tag const & tag)
129  {
130 
131  typedef typename viennacl::result_of::value_type<MatrixT>::type ScalarType;
132  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
133 
134 
135  // generation of some random numbers, used for lanczos PRO algorithm
136  boost::mt11213b mt;
137  boost::normal_distribution<CPU_ScalarType> N(0, 1);
138  boost::bernoulli_distribution<CPU_ScalarType> B(0.5);
139  boost::triangle_distribution<CPU_ScalarType> T(-1, 0, 1);
140 
141  boost::variate_generator<boost::mt11213b&, boost::normal_distribution<CPU_ScalarType> > get_N(mt, N);
142  boost::variate_generator<boost::mt11213b&, boost::bernoulli_distribution<CPU_ScalarType> > get_B(mt, B);
143  boost::variate_generator<boost::mt11213b&, boost::triangle_distribution<CPU_ScalarType> > get_T(mt, T);
144 
145 
146  long i, j, k, index, retry, reorths;
147  std::vector<long> l_bound(size/2), u_bound(size/2);
148  bool second_step;
149  CPU_ScalarType squ_eps, eta, temp, eps, retry_th;
150  vcl_size_t n = r.size();
151  std::vector< std::vector<CPU_ScalarType> > w(2, std::vector<CPU_ScalarType>(size));
152  CPU_ScalarType cpu_beta;
153 
154  boost::numeric::ublas::vector<CPU_ScalarType> s(n);
155 
156  VectorT t(n);
157  CPU_ScalarType inner_rt;
158  ScalarType vcl_beta;
159  ScalarType vcl_alpha;
160  std::vector<CPU_ScalarType> alphas, betas;
161  boost::numeric::ublas::matrix<CPU_ScalarType> Q(n, size);
162 
163  second_step = false;
164  eps = std::numeric_limits<CPU_ScalarType>::epsilon();
165  squ_eps = std::sqrt(eps);
166  retry_th = 1e-2;
167  eta = std::exp(std::log(eps) * tag.factor());
168  reorths = 0;
169  retry = 0;
170 
171  vcl_beta = viennacl::linalg::norm_2(r);
172 
173  r /= vcl_beta;
174 
177 
178  VectorT u = viennacl::linalg::prod(A, r);
179  vcl_alpha = viennacl::linalg::inner_prod(u, r);
180  alphas.push_back(vcl_alpha);
181  w[0][0] = 1;
182  betas.push_back(vcl_beta);
183 
184  long batches = 0;
185  for(i = 1;i < static_cast<long>(size); i++)
186  {
187  r = u - vcl_alpha * r;
188  vcl_beta = viennacl::linalg::norm_2(r);
189 
190  betas.push_back(vcl_beta);
191  r = r / vcl_beta;
192 
193  index = i % 2;
194  w[index][i] = 1;
195  k = (i + 1) % 2;
196  w[index][0] = (betas[1] * w[k][1] + (alphas[0] - vcl_alpha) * w[k][0] - betas[i - 1] * w[index][0]) / vcl_beta + eps * 0.3 * get_N() * (betas[1] + vcl_beta);
197 
198  for(j = 1;j < i - 1;j++)
199  {
200  w[index][j] = (betas[j + 1] * w[k][j + 1] + (alphas[j] - vcl_alpha) * w[k][j] + betas[j] * w[k][j - 1] - betas[i - 1] * w[index][j]) / vcl_beta + eps * 0.3 * get_N() * (betas[j + 1] + vcl_beta);
201  }
202  w[index][i - 1] = 0.6 * eps * n * get_N() * betas[1] / vcl_beta;
203 
204  if(second_step)
205  {
206  for(j = 0;j < batches;j++)
207  {
208  l_bound[j]++;
209  u_bound[j]--;
210 
211  for(k = l_bound[j];k < u_bound[j];k++)
212  {
214  inner_rt = viennacl::linalg::inner_prod(r,t);
215  r = r - inner_rt * t;
216  w[index][k] = 1.5 * eps * get_N();
217  reorths++;
218  }
219  }
220  temp = viennacl::linalg::norm_2(r);
221  r = r / temp;
222  vcl_beta = vcl_beta * temp;
223  second_step = false;
224  }
225  batches = 0;
226 
227  for(j = 0;j < i;j++)
228  {
229  if(std::fabs(w[index][j]) >= squ_eps)
230  {
232  inner_rt = viennacl::linalg::inner_prod(r,t);
233  r = r - inner_rt * t;
234  w[index][j] = 1.5 * eps * get_N();
235  k = j - 1;
236  reorths++;
237  while(k >= 0 && std::fabs(w[index][k]) > eta)
238  {
240  inner_rt = viennacl::linalg::inner_prod(r,t);
241  r = r - inner_rt * t;
242  w[index][k] = 1.5 * eps * get_N();
243  k--;
244  reorths++;
245  }
246  l_bound[batches] = k + 1;
247  k = j + 1;
248 
249  while(k < i && std::fabs(w[index][k]) > eta)
250  {
252  inner_rt = viennacl::linalg::inner_prod(r,t);
253  r = r - inner_rt * t;
254  w[index][k] = 1.5 * eps * get_N();
255  k++;
256  reorths++;
257  }
258  u_bound[batches] = k - 1;
259  batches++;
260  j = k;
261  }
262  }
263 
264  if(batches > 0)
265  {
266  temp = viennacl::linalg::norm_2(r);
267  r = r / temp;
268  vcl_beta = vcl_beta * temp;
269  second_step = true;
270 
271  while(temp < retry_th)
272  {
273  for(j = 0;j < i;j++)
274  {
276  inner_rt = viennacl::linalg::inner_prod(r,t);
277  r = r - inner_rt * t;
278  reorths++;
279  }
280  retry++;
281  temp = viennacl::linalg::norm_2(r);
282  r = r / temp;
283  vcl_beta = vcl_beta * temp;
284  }
285  }
286 
289 
290  cpu_beta = vcl_beta;
291  s = - cpu_beta * boost::numeric::ublas::column(Q, i - 1);
293  u += viennacl::linalg::prod(A, r);
294  vcl_alpha = viennacl::linalg::inner_prod(u, r);
295  alphas.push_back(vcl_alpha);
296  }
297 
298  return bisect(alphas, betas);
299 
300  }
301 
302 
311  template< typename MatrixT, typename VectorT >
312  std::vector<
314  >
315  lanczos (MatrixT const& A, VectorT & r, vcl_size_t size, lanczos_tag)
316  {
317 
318  typedef typename viennacl::result_of::value_type<MatrixT>::type ScalarType;
319  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
320 
321  ScalarType vcl_beta;
322  ScalarType vcl_alpha;
323  std::vector<CPU_ScalarType> alphas, betas;
324  CPU_ScalarType norm;
325  vcl_size_t n = r.size();
326  VectorT u(n), t(n);
327  boost::numeric::ublas::vector<CPU_ScalarType> s(r.size()), u_zero(n), q(n);
328  boost::numeric::ublas::matrix<CPU_ScalarType> Q(n, size);
329 
330  u_zero = boost::numeric::ublas::zero_vector<CPU_ScalarType>(n);
331  detail::copy_vec_to_vec(u_zero, u);
332  norm = norm_2(r);
333 
334  for(vcl_size_t i = 0;i < size; i++)
335  {
336  r /= norm;
337  vcl_beta = norm;
338 
341 
342  u += prod(A, r);
343  vcl_alpha = inner_prod(u, r);
344  r = u - vcl_alpha * r;
345  norm = norm_2(r);
346 
349 
350  u = - norm * t;
351  alphas.push_back(vcl_alpha);
352  betas.push_back(vcl_beta);
353  s.clear();
354  }
355 
356  return bisect(alphas, betas);
357  }
358 
367  template< typename MatrixT, typename VectorT >
368  std::vector<
370  >
371  lanczosFRO (MatrixT const& A, VectorT & r, vcl_size_t size, lanczos_tag)
372  {
373 
374  typedef typename viennacl::result_of::value_type<MatrixT>::type ScalarType;
375  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
376 
377  CPU_ScalarType temp;
378  CPU_ScalarType norm;
379  ScalarType vcl_beta;
380  ScalarType vcl_alpha;
381  std::vector<CPU_ScalarType> alphas, betas;
382  vcl_size_t n = r.size();
383  VectorT u(n), t(n);
384  ScalarType inner_rt;
385  boost::numeric::ublas::vector<CPU_ScalarType> u_zero(n), s(r.size()), q(n);
386  boost::numeric::ublas::matrix<CPU_ScalarType> Q(n, size);
387 
388  long reorths = 0;
389  norm = norm_2(r);
390 
391 
392  for(vcl_size_t i = 0; i < size; i++)
393  {
394  r /= norm;
395 
396  for(vcl_size_t j = 0; j < i; j++)
397  {
400  inner_rt = viennacl::linalg::inner_prod(r,t);
401  r = r - inner_rt * t;
402  reorths++;
403  }
404  temp = viennacl::linalg::norm_2(r);
405  r = r / temp;
406  vcl_beta = temp * norm;
409 
410  u += viennacl::linalg::prod(A, r);
411  vcl_alpha = viennacl::linalg::inner_prod(u, r);
412  r = u - vcl_alpha * r;
413  norm = viennacl::linalg::norm_2(r);
416  u = - norm * t;
417  alphas.push_back(vcl_alpha);
418  betas.push_back(vcl_beta);
419  }
420 
421  return bisect(alphas, betas);
422  }
423 
424  } // end namespace detail
425 
433  template< typename MatrixT >
434  std::vector< typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type >
435  eig(MatrixT const & matrix, lanczos_tag const & tag)
436  {
437  typedef typename viennacl::result_of::value_type<MatrixT>::type ScalarType;
438  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
439  typedef typename viennacl::result_of::vector_for_matrix<MatrixT>::type VectorT;
440 
441  boost::mt11213b mt;
442  boost::normal_distribution<CPU_ScalarType> N(0, 1);
443  boost::bernoulli_distribution<CPU_ScalarType> B(0.5);
444  boost::triangle_distribution<CPU_ScalarType> T(-1, 0, 1);
445 
446  boost::variate_generator<boost::mt11213b&, boost::normal_distribution<CPU_ScalarType> > get_N(mt, N);
447  boost::variate_generator<boost::mt11213b&, boost::bernoulli_distribution<CPU_ScalarType> > get_B(mt, B);
448  boost::variate_generator<boost::mt11213b&, boost::triangle_distribution<CPU_ScalarType> > get_T(mt, T);
449 
450  std::vector<CPU_ScalarType> eigenvalues;
451  vcl_size_t matrix_size = matrix.size1();
452  VectorT r(matrix_size);
453  std::vector<CPU_ScalarType> s(matrix_size);
454 
455  for(vcl_size_t i=0; i<s.size(); ++i)
456  s[i] = 3.0 * get_B() + get_T() - 1.5;
457 
459 
460  vcl_size_t size_krylov = (matrix_size < tag.krylov_size()) ? matrix_size
461  : tag.krylov_size();
462 
463  switch(tag.method())
464  {
466  eigenvalues = detail::lanczosPRO(matrix, r, size_krylov, tag);
467  break;
469  eigenvalues = detail::lanczosFRO(matrix, r, size_krylov, tag);
470  break;
472  eigenvalues = detail::lanczos(matrix, r, size_krylov, tag);
473  break;
474  }
475 
476  std::vector<CPU_ScalarType> largest_eigenvalues;
477 
478  for(vcl_size_t i = 1; i<=tag.num_eigenvalues(); i++)
479  largest_eigenvalues.push_back(eigenvalues[size_krylov-i]);
480 
481 
482  return largest_eigenvalues;
483  }
484 
485 
486 
487 
488  } // end namespace linalg
489 } // end namespace viennacl
490 #endif
std::size_t vcl_size_t
Definition: forwards.h:58
T norm_2(std::vector< T, A > const &v1)
Definition: norm_2.hpp:86
A reader and writer for the matrix market format is implemented here.
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
void factor(double fct)
Sets the exponent of epsilon.
Definition: lanczos.hpp:85
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixT::value_type >::type > lanczosPRO(MatrixT const &A, VectorT &r, vcl_size_t size, lanczos_tag const &tag)
Implementation of the Lanczos PRO algorithm.
Definition: lanczos.hpp:128
A dense matrix class.
Definition: forwards.h:293
viennacl::enable_if< viennacl::is_stl< typename viennacl::traits::tag_of< VectorT1 >::type >::value, typename VectorT1::value_type >::type inner_prod(VectorT1 const &v1, VectorT2 const &v2)
Definition: inner_prod.hpp:89
void num_eigenvalues(int numeig)
Sets the number of eigenvalues.
Definition: lanczos.hpp:79
Generic interface for the computation of inner products. See viennacl/linalg/vector_operations.hpp for implementations.
std::vector< typename viennacl::result_of::cpu_value_type< typename VectorT::value_type >::type > bisect(VectorT const &alphas, VectorT const &betas)
Implementation of the bisect-algorithm for the calculation of the eigenvalues of a tridiagonal matrix...
Definition: bisect.hpp:74
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
vcl_size_t krylov_size() const
Returns the size of the kylov space.
Definition: lanczos.hpp:94
std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixT::value_type >::type > lanczos(MatrixT const &A, VectorT &r, vcl_size_t size, lanczos_tag)
Implementation of the lanczos algorithm without reorthogonalization.
Definition: lanczos.hpp:315
Implementation of the compressed_matrix class.
void krylov_size(int max)
Sets the size of the kylov space.
Definition: lanczos.hpp:91
void method(int met)
Sets the reorthogonalization method.
Definition: lanczos.hpp:97
T::value_type type
Definition: result_of.hpp:230
lanczos_tag(double factor=0.75, vcl_size_t numeig=10, int met=0, vcl_size_t krylov=100)
The constructor.
Definition: lanczos.hpp:73
std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixT::value_type >::type > lanczosFRO(MatrixT const &A, VectorT &r, vcl_size_t size, lanczos_tag)
Implementation of the Lanczos FRO algorithm.
Definition: lanczos.hpp:371
std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixT::value_type >::type > eig(MatrixT const &matrix, lanczos_tag const &tag)
Implementation of the calculation of eigenvalues using lanczos.
Definition: lanczos.hpp:435
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:276
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
double factor() const
Returns the exponent.
Definition: lanczos.hpp:88
Implementation of the algorithm for finding eigenvalues of a tridiagonal matrix.
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_column > column(const matrix_base< NumericT, F > &A, unsigned int j)
Definition: matrix.hpp:918
void copy_vec_to_vec(viennacl::vector< T > const &src, OtherVectorType &dest)
overloaded function for copying vectors
Definition: bisect.hpp:44
A tag for the lanczos algorithm.
Definition: lanczos.hpp:54
int method() const
Returns the reorthogonalization method.
Definition: lanczos.hpp:100
vcl_size_t num_eigenvalues() const
Returns the number of eigenvalues.
Definition: lanczos.hpp:82