1 #ifndef VIENNACL_LINALG_LANCZOS_HPP_
2 #define VIENNACL_LINALG_LANCZOS_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>
76 vcl_size_t krylov = 100) : factor_(
factor), num_eigenvalues_(numeig), method_(met), krylov_size_(krylov) {}
85 void factor(
double fct) { factor_ = fct; }
88 double factor()
const {
return factor_; }
97 void method(
int met){ method_ = met; }
124 template<
typename MatrixT,
typename VectorT >
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);
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);
146 long i, j, k, index, retry, reorths;
147 std::vector<long> l_bound(size/2), u_bound(size/2);
149 CPU_ScalarType squ_eps, eta, temp, eps, retry_th;
151 std::vector< std::vector<CPU_ScalarType> > w(2, std::vector<CPU_ScalarType>(size));
152 CPU_ScalarType cpu_beta;
154 boost::numeric::ublas::vector<CPU_ScalarType> s(n);
157 CPU_ScalarType inner_rt;
159 ScalarType vcl_alpha;
160 std::vector<CPU_ScalarType> alphas, betas;
161 boost::numeric::ublas::matrix<CPU_ScalarType> Q(n, size);
164 eps = std::numeric_limits<CPU_ScalarType>::epsilon();
165 squ_eps = std::sqrt(eps);
167 eta = std::exp(std::log(eps) * tag.
factor());
180 alphas.push_back(vcl_alpha);
182 betas.push_back(vcl_beta);
185 for(i = 1;i < static_cast<long>(
size); i++)
187 r = u - vcl_alpha * r;
190 betas.push_back(vcl_beta);
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);
198 for(j = 1;j < i - 1;j++)
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);
202 w[index][i - 1] = 0.6 * eps * n * get_N() * betas[1] / vcl_beta;
206 for(j = 0;j < batches;j++)
211 for(k = l_bound[j];k < u_bound[j];k++)
215 r = r - inner_rt * t;
216 w[index][k] = 1.5 * eps * get_N();
222 vcl_beta = vcl_beta * temp;
229 if(std::fabs(w[index][j]) >= squ_eps)
233 r = r - inner_rt * t;
234 w[index][j] = 1.5 * eps * get_N();
237 while(k >= 0 && std::fabs(w[index][k]) > eta)
241 r = r - inner_rt * t;
242 w[index][k] = 1.5 * eps * get_N();
246 l_bound[batches] = k + 1;
249 while(k < i && std::fabs(w[index][k]) > eta)
253 r = r - inner_rt * t;
254 w[index][k] = 1.5 * eps * get_N();
258 u_bound[batches] = k - 1;
268 vcl_beta = vcl_beta * temp;
271 while(temp < retry_th)
277 r = r - inner_rt * t;
283 vcl_beta = vcl_beta * temp;
295 alphas.push_back(vcl_alpha);
298 return bisect(alphas, betas);
311 template<
typename MatrixT,
typename VectorT >
322 ScalarType vcl_alpha;
323 std::vector<CPU_ScalarType> alphas, betas;
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);
330 u_zero = boost::numeric::ublas::zero_vector<CPU_ScalarType>(n);
344 r = u - vcl_alpha * r;
351 alphas.push_back(vcl_alpha);
352 betas.push_back(vcl_beta);
356 return bisect(alphas, betas);
367 template<
typename MatrixT,
typename VectorT >
380 ScalarType vcl_alpha;
381 std::vector<CPU_ScalarType> alphas, betas;
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);
401 r = r - inner_rt * t;
406 vcl_beta = temp * norm;
412 r = u - vcl_alpha * r;
417 alphas.push_back(vcl_alpha);
418 betas.push_back(vcl_beta);
421 return bisect(alphas, betas);
433 template<
typename MatrixT >
434 std::vector< typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type >
439 typedef typename viennacl::result_of::vector_for_matrix<MatrixT>::type VectorT;
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);
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);
450 std::vector<CPU_ScalarType> eigenvalues;
452 VectorT r(matrix_size);
453 std::vector<CPU_ScalarType> s(matrix_size);
456 s[i] = 3.0 * get_B() + get_T() - 1.5;
476 std::vector<CPU_ScalarType> largest_eigenvalues;
479 largest_eigenvalues.push_back(eigenvalues[size_krylov-i]);
482 return largest_eigenvalues;
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
Definition: lanczos.hpp:61
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
Definition: lanczos.hpp:62
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 ...
Definition: lanczos.hpp:60
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