ViennaCL - The Vienna Computing Library  1.7.1
Free open-source GPU-accelerated linear algebra and solver library.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
tql.cpp
Go to the documentation of this file.
1 /* =========================================================================
2  Copyright (c) 2010-2016, Institute for Microelectronics,
3  Institute for Analysis and Scientific Computing,
4  TU Wien.
5  Portions of this software are copyright by UChicago Argonne, LLC.
6 
7  -----------------
8  ViennaCL - The Vienna Computing Library
9  -----------------
10 
11  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12 
13  (A list of authors and contributors can be found in the PDF manual)
14 
15  License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17 
18 
19 
24 /*
25 *
26 * Test file for tql-algorithm
27 *
28 */
29 
30 // include necessary system headers
31 #include <iostream>
32 
33 
34 //include basic scalar and vector types of ViennaCL
35 #include "viennacl/scalar.hpp"
36 #include "viennacl/vector.hpp"
37 #include "viennacl/linalg/tql2.hpp"
38 
39 #define EPS 10.0e-5
40 
41 
42 typedef float ScalarType;
43 
44 
45 // Test the eigenvectors
46 // Perform the multiplication (T - lambda * I) * Q, with the original tridiagonal matrx T, the
47 // eigenvalues lambda and the eigenvectors in Q. Result has to be 0.
48 
49 template <typename MatrixLayout>
51  std::vector<ScalarType> & eigenvalues,
52  std::vector<ScalarType> & d,
53  std::vector<ScalarType> & e)
54 {
55  std::size_t Q_size = Q.size2();
56  ScalarType value = 0;
57 
58  for(std::size_t j = 0; j < Q_size; j++)
59  {
60  // calculate first row
61  value = (d[0]- eigenvalues[j]) * Q(0, j) + e[1] * Q(1, j);
62  if (value > EPS)
63  return false;
64 
65  // calcuate inner rows
66  for(std::size_t i = 1; i < Q_size - 1; i++)
67  {
68  value = e[i] * Q(i - 1, j) + (d[i]- eigenvalues[j]) * Q(i, j) + e[i + 1] * Q(i + 1, j);
69  if (value > EPS)
70  return false;
71  }
72 
73  // calculate last row
74  value = e[Q_size - 1] * Q(Q_size - 2, j) + (d[Q_size - 1] - eigenvalues[j]) * Q(Q_size - 1, j);
75  if (value > EPS)
76  return false;
77  }
78  return true;
79 }
80 
81 
86 template <typename MatrixLayout>
88 {
89  std::size_t sz = 220;
90 
92  std::vector<ScalarType> d(sz), e(sz), d_ref(sz), e_ref(sz);
93 
94  std::cout << "Testing matrix of size " << sz << "-by-" << sz << std::endl << std::endl;
95 
96  // Initialize diagonal and superdiagonal elements
97  for(unsigned int i = 0; i < sz; ++i)
98  {
99  d[i] = ((float)(i % 9)) - 4.5f;
100  e[i] = ((float)(i % 5)) - 4.5f;
101  }
102  e[0] = 0.0f;
103  d_ref = d;
104  e_ref = e;
105 
106 //---Run the tql2 algorithm-----------------------------------
107  viennacl::linalg::tql2(Q, d, e);
108 
109 
110 // ---Test the computed eigenvalues and eigenvectors
111  if(!test_eigen_val_vec<MatrixLayout>(Q, d, d_ref, e_ref))
112  exit(EXIT_FAILURE);
113 /*
114  for( unsigned int i = 0; i < sz; ++i)
115  std::cout << "Eigenvalue " << i << "= " << d[i] << std::endl;
116  */
117 }
118 
119 int main()
120 {
121 
122  std::cout << std::endl << "COMPUTATION OF EIGENVALUES AND EIGENVECTORS" << std::endl;
123  std::cout << std::endl << "Testing QL algorithm for symmetric tridiagonal row-major matrices..." << std::endl;
124  test_qr_method_sym<viennacl::row_major>();
125 
126  std::cout << std::endl << "Testing QL algorithm for symmetric tridiagonal column-major matrices..." << std::endl;
127  test_qr_method_sym<viennacl::column_major>();
128 
129  std::cout << std::endl <<"--------TEST SUCCESSFULLY COMPLETED----------" << std::endl;
130 }
void tql2(matrix_base< SCALARTYPE, F > &Q, VectorType &d, VectorType &e)
Definition: tql2.hpp:131
Implementation of the tql2-algorithm for eigenvalue computations.
A dense matrix class.
Definition: forwards.h:375
Represents a vector consisting of 1 at a given index and zeros otherwise. To be used as an initialize...
Definition: matrix_def.hpp:69
void test_qr_method_sym()
Definition: tql.cpp:87
bool test_eigen_val_vec(viennacl::matrix< ScalarType, MatrixLayout > &Q, std::vector< ScalarType > &eigenvalues, std::vector< ScalarType > &d, std::vector< ScalarType > &e)
Definition: tql.cpp:50
float ScalarType
Definition: tql.cpp:42
size_type size2() const
Returns the number of columns.
Definition: matrix_def.hpp:226
int main()
Definition: tql.cpp:119
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
float ScalarType
Definition: fft_1d.cpp:42
#define EPS
Definition: tql.cpp:39
Implementation of the ViennaCL scalar class.