ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
bisect.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_BISECT_HPP_
2 #define VIENNACL_LINALG_BISECT_HPP_
3 
4 /* =========================================================================
5  Copyright (c) 2010-2015, 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 manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
27 #include <vector>
28 #include <cmath>
29 #include <limits>
30 #include <cstddef>
32 
33 namespace viennacl
34 {
35 namespace linalg
36 {
37 
38 namespace detail
39 {
43  template<typename NumericT, typename OtherVectorT>
44  void copy_vec_to_vec(viennacl::vector<NumericT> const & src, OtherVectorT & dest)
45  {
46  viennacl::copy(src, dest);
47  }
48 
49  template<typename OtherVectorT, typename NumericT>
50  void copy_vec_to_vec(OtherVectorT const & src, viennacl::vector<NumericT> & dest)
51  {
52  viennacl::copy(src, dest);
53  }
54 
55  template<typename VectorT1, typename VectorT2>
56  void copy_vec_to_vec(VectorT1 const & src, VectorT2 & dest)
57  {
58  for (vcl_size_t i=0; i<src.size(); ++i)
59  dest[i] = src[i];
60  }
61 
62 } //namespace detail
63 
74 template<typename VectorT>
75 std::vector<
77  >
78 bisect(VectorT const & alphas, VectorT const & betas)
79 {
80  typedef typename viennacl::result_of::value_type<VectorT>::type NumericType;
81  typedef typename viennacl::result_of::cpu_value_type<NumericType>::type CPU_NumericType;
82 
83  vcl_size_t size = betas.size();
84  std::vector<CPU_NumericType> x_temp(size);
85 
86 
87  std::vector<CPU_NumericType> beta_bisect;
88  std::vector<CPU_NumericType> wu;
89 
90  double rel_error = std::numeric_limits<CPU_NumericType>::epsilon();
91  beta_bisect.push_back(0);
92 
93  for (vcl_size_t i = 1; i < size; i++)
94  beta_bisect.push_back(betas[i] * betas[i]);
95 
96  double xmin = alphas[size - 1] - std::fabs(betas[size - 1]);
97  double xmax = alphas[size - 1] + std::fabs(betas[size - 1]);
98 
99  for (vcl_size_t i = 0; i < size - 1; i++)
100  {
101  double h = std::fabs(betas[i]) + std::fabs(betas[i + 1]);
102  if (alphas[i] + h > xmax)
103  xmax = alphas[i] + h;
104  if (alphas[i] - h < xmin)
105  xmin = alphas[i] - h;
106  }
107 
108 
109  double eps1 = 1e-6;
110  /*double eps2 = (xmin + xmax > 0) ? (rel_error * xmax) : (-rel_error * xmin);
111  if (eps1 <= 0)
112  eps1 = eps2;
113  else
114  eps2 = 0.5 * eps1 + 7.0 * eps2; */
115 
116  double x0 = xmax;
117 
118  for (vcl_size_t i = 0; i < size; i++)
119  {
120  x_temp[i] = xmax;
121  wu.push_back(xmin);
122  }
123 
124  for (long k = static_cast<long>(size) - 1; k >= 0; --k)
125  {
126  double xu = xmin;
127  for (long i = k; i >= 0; --i)
128  {
129  if (xu < wu[vcl_size_t(k-i)])
130  {
131  xu = wu[vcl_size_t(i)];
132  break;
133  }
134  }
135 
136  if (x0 > x_temp[vcl_size_t(k)])
137  x0 = x_temp[vcl_size_t(k)];
138 
139  double x1 = (xu + x0) / 2.0;
140  while (x0 - xu > 2.0 * rel_error * (std::fabs(xu) + std::fabs(x0)) + eps1)
141  {
142  vcl_size_t a = 0;
143  double q = 1;
144  for (vcl_size_t i = 0; i < size; i++)
145  {
146  if (q > 0 || q < 0)
147  q = alphas[i] - x1 - beta_bisect[i] / q;
148  else
149  q = alphas[i] - x1 - std::fabs(betas[i] / rel_error);
150 
151  if (q < 0)
152  a++;
153  }
154 
155  if (a <= static_cast<vcl_size_t>(k))
156  {
157  xu = x1;
158  if (a < 1)
159  wu[0] = x1;
160  else
161  {
162  wu[a] = x1;
163  if (x_temp[a - 1] > x1)
164  x_temp[a - 1] = x1;
165  }
166  }
167  else
168  x0 = x1;
169 
170  x1 = (xu + x0) / 2.0;
171  }
172  x_temp[vcl_size_t(k)] = x1;
173  }
174  return x_temp;
175 }
176 
177 } // end namespace linalg
178 } // end namespace viennacl
179 #endif
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:78
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:34
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:235
Definition: blas3.hpp:36
void copy_vec_to_vec(viennacl::vector< NumericT > const &src, OtherVectorT &dest)
overloaded function for copying vectors
Definition: bisect.hpp:44
std::size_t vcl_size_t
Definition: forwards.h:75
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:271
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
A collection of compile time type deductions.