ViennaCL - The Vienna Computing Library  1.5.2
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-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 <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 T, typename OtherVectorType>
44  void copy_vec_to_vec(viennacl::vector<T> const & src, OtherVectorType & dest)
45  {
46  viennacl::copy(src, dest);
47  }
48 
49  template <typename OtherVectorType, typename T>
50  void copy_vec_to_vec(OtherVectorType const & src, viennacl::vector<T> & dest)
51  {
52  viennacl::copy(src, dest);
53  }
54 
55  template <typename VectorType1, typename VectorType2>
56  void copy_vec_to_vec(VectorType1 const & src, VectorType2 & dest)
57  {
58  for (vcl_size_t i=0; i<src.size(); ++i)
59  dest[i] = src[i];
60  }
61  }
62 
70  template< typename VectorT >
71  std::vector<
73  >
74  bisect(VectorT const & alphas, VectorT const & betas)
75  {
76  typedef typename viennacl::result_of::value_type<VectorT>::type ScalarType;
77  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
78 
79  vcl_size_t size = betas.size();
80  std::vector<CPU_ScalarType> x_temp(size);
81 
82 
83  std::vector<CPU_ScalarType> beta_bisect;
84  std::vector<CPU_ScalarType> wu;
85 
86  double rel_error = std::numeric_limits<CPU_ScalarType>::epsilon();
87  beta_bisect.push_back(0);
88 
89  for(vcl_size_t i = 1; i < size; i++){
90  beta_bisect.push_back(betas[i] * betas[i]);
91  }
92 
93  double xmin = alphas[size - 1] - std::fabs(betas[size - 1]);
94  double xmax = alphas[size - 1] + std::fabs(betas[size - 1]);
95 
96  for(vcl_size_t i = 0; i < size - 1; i++)
97  {
98  double h = std::fabs(betas[i]) + std::fabs(betas[i + 1]);
99  if (alphas[i] + h > xmax)
100  xmax = alphas[i] + h;
101  if (alphas[i] - h < xmin)
102  xmin = alphas[i] - h;
103  }
104 
105 
106  double eps1 = 1e-6;
107  /*double eps2 = (xmin + xmax > 0) ? (rel_error * xmax) : (-rel_error * xmin);
108  if(eps1 <= 0)
109  eps1 = eps2;
110  else
111  eps2 = 0.5 * eps1 + 7.0 * eps2; */
112 
113  double x0 = xmax;
114 
115  for(vcl_size_t i = 0; i < size; i++)
116  {
117  x_temp[i] = xmax;
118  wu.push_back(xmin);
119  }
120 
121  for(long k = static_cast<long>(size) - 1; k >= 0; --k)
122  {
123  double xu = xmin;
124  for(long i = k; i >= 0; --i)
125  {
126  if(xu < wu[k-i])
127  {
128  xu = wu[i];
129  break;
130  }
131  }
132 
133  if(x0 > x_temp[k])
134  x0 = x_temp[k];
135 
136  double x1 = (xu + x0) / 2.0;
137  while (x0 - xu > 2.0 * rel_error * (std::fabs(xu) + std::fabs(x0)) + eps1)
138  {
139  vcl_size_t a = 0;
140  double q = 1;
141  for(vcl_size_t i = 0; i < size; i++)
142  {
143  if(q != 0)
144  q = alphas[i] - x1 - beta_bisect[i] / q;
145  else
146  q = alphas[i] - x1 - std::fabs(betas[i] / rel_error);
147 
148  if(q < 0)
149  a++;
150  }
151 
152  if (a <= static_cast<vcl_size_t>(k))
153  {
154  xu = x1;
155  if(a < 1)
156  wu[0] = x1;
157  else
158  {
159  wu[a] = x1;
160  if(x_temp[a - 1] > x1)
161  x_temp[a - 1] = x1;
162  }
163  }
164  else
165  x0 = x1;
166 
167  x1 = (xu + x0) / 2.0;
168  }
169  x_temp[k] = x1;
170  }
171  return x_temp;
172  }
173 
174  } // end namespace linalg
175 } // end namespace viennacl
176 #endif
std::size_t vcl_size_t
Definition: forwards.h:58
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
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
void copy(std::vector< SCALARTYPE > &cpu_vec, circulant_matrix< SCALARTYPE, ALIGNMENT > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
T::value_type type
Definition: result_of.hpp:230
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:276
void copy_vec_to_vec(viennacl::vector< T > const &src, OtherVectorType &dest)
overloaded function for copying vectors
Definition: bisect.hpp:44
A collection of compile time type deductions.