All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
tql2.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_TQL2_HPP
#define VIENNACL_LINALG_TQL2_HPP
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include "viennacl/vector.hpp"
#include "viennacl/matrix.hpp"
#include <iomanip>
#include "viennacl/linalg/qr-method-common.hpp"
#include "viennacl/linalg/prod.hpp"
namespace viennacl
{
namespace linalg
{
// Symmetric tridiagonal QL algorithm.
// This is derived from the Algol procedures tql1, by Bowdler, Martin, Reinsch, and Wilkinson,
// Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding Fortran subroutine in EISPACK.
template <typename SCALARTYPE, typename VectorType>
void tql1(vcl_size_t n,
VectorType & d,
VectorType & e)
{
for (vcl_size_t i = 1; i < n; i++)
e[i - 1] = e[i];
e[n - 1] = 0;
SCALARTYPE f = 0.;
SCALARTYPE tst1 = 0.;
SCALARTYPE eps = static_cast<SCALARTYPE>(1e-6);
for (vcl_size_t l = 0; l < n; l++)
{
// Find small subdiagonal element.
tst1 = std::max<SCALARTYPE>(tst1, std::fabs(d[l]) + std::fabs(e[l]));
vcl_size_t m = l;
while (m < n)
{
if (std::fabs(e[m]) <= eps * tst1)
break;
m++;
}
// If m == l, d[l) is an eigenvalue, otherwise, iterate.
if (m > l)
{
vcl_size_t iter = 0;
do
{
iter = iter + 1; // (Could check iteration count here.)
// Compute implicit shift
SCALARTYPE g = d[l];
SCALARTYPE p = (d[l + 1] - g) / (2 * e[l]);
SCALARTYPE r = viennacl::linalg::detail::pythag<SCALARTYPE>(p, 1);
if (p < 0)
{
r = -r;
}
d[l] = e[l] / (p + r);
d[l + 1] = e[l] * (p + r);
SCALARTYPE h = g - d[l];
for (vcl_size_t i = l + 2; i < n; i++)
{
d[i] -= h;
}
f = f + h;
// Implicit QL transformation.
p = d[m];
SCALARTYPE c = 1;
SCALARTYPE s = 0;
for (int i = int(m - 1); i >= int(l); i--)
{
g = c * e[vcl_size_t(i)];
h = c * p;
r = viennacl::linalg::detail::pythag<SCALARTYPE>(p, e[vcl_size_t(i)]);
e[vcl_size_t(i) + 1] = s * r;
s = e[vcl_size_t(i)] / r;
c = p / r;
p = c * d[vcl_size_t(i)] - s * g;
d[vcl_size_t(i) + 1] = h + s * (c * g + s * d[vcl_size_t(i)]);
}
e[l] = s * p;
d[l] = c * p;
// Check for convergence.
}
while (std::fabs(e[l]) > eps * tst1);
}
d[l] = d[l] + f;
e[l] = 0;
}
}
// Symmetric tridiagonal QL algorithm.
// This is derived from the Algol procedures tql2, by Bowdler, Martin, Reinsch, and Wilkinson,
// Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding Fortran subroutine in EISPACK.
template <typename SCALARTYPE, typename VectorType, typename F>
void tql2(matrix_base<SCALARTYPE, F> & Q,
VectorType & d,
VectorType & e)
{
vcl_size_t n = static_cast<vcl_size_t>(viennacl::traits::size1(Q));
//boost::numeric::ublas::vector<SCALARTYPE> cs(n), ss(n);
std::vector<SCALARTYPE> cs(n), ss(n);
viennacl::vector<SCALARTYPE> tmp1(n), tmp2(n);
for (vcl_size_t i = 1; i < n; i++)
e[i - 1] = e[i];
e[n - 1] = 0;
SCALARTYPE f = 0;
SCALARTYPE tst1 = 0;
SCALARTYPE eps = static_cast<SCALARTYPE>(viennacl::linalg::detail::EPS);
for (vcl_size_t l = 0; l < n; l++)
{
// Find small subdiagonal element.
tst1 = std::max<SCALARTYPE>(tst1, std::fabs(d[l]) + std::fabs(e[l]));
vcl_size_t m = l;
while (m < n)
{
if (std::fabs(e[m]) <= eps * tst1)
break;
m++;
}
// If m == l, d[l) is an eigenvalue, otherwise, iterate.
if (m > l)
{
vcl_size_t iter = 0;
do
{
iter = iter + 1; // (Could check iteration count here.)
// Compute implicit shift
SCALARTYPE g = d[l];
SCALARTYPE p = (d[l + 1] - g) / (2 * e[l]);
SCALARTYPE r = viennacl::linalg::detail::pythag<SCALARTYPE>(p, 1);
if (p < 0)
{
r = -r;
}
d[l] = e[l] / (p + r);
d[l + 1] = e[l] * (p + r);
SCALARTYPE dl1 = d[l + 1];
SCALARTYPE h = g - d[l];
for (vcl_size_t i = l + 2; i < n; i++)
{
d[i] -= h;
}
f = f + h;
// Implicit QL transformation.
p = d[m];
SCALARTYPE c = 1;
SCALARTYPE c2 = c;
SCALARTYPE c3 = c;
SCALARTYPE el1 = e[l + 1];
SCALARTYPE s = 0;
SCALARTYPE s2 = 0;
for (int i = int(m - 1); i >= int(l); i--)
{
c3 = c2;
c2 = c;
s2 = s;
g = c * e[vcl_size_t(i)];
h = c * p;
r = viennacl::linalg::detail::pythag(p, e[vcl_size_t(i)]);
e[vcl_size_t(i) + 1] = s * r;
s = e[vcl_size_t(i)] / r;
c = p / r;
p = c * d[vcl_size_t(i)] - s * g;
d[vcl_size_t(i) + 1] = h + s * (c * g + s * d[vcl_size_t(i)]);
cs[vcl_size_t(i)] = c;
ss[vcl_size_t(i)] = s;
}
p = -s * s2 * c3 * el1 * e[l] / dl1;
e[l] = s * p;
d[l] = c * p;
viennacl::copy(cs, tmp1);
viennacl::copy(ss, tmp2);
viennacl::linalg::givens_next(Q, tmp1, tmp2, int(l), int(m));
// Check for convergence.
}
while (std::fabs(e[l]) > eps * tst1);
}
d[l] = d[l] + f;
e[l] = 0;
}
// Sort eigenvalues and corresponding vectors.
/*
for (int i = 0; i < n-1; i++) {
int k = i;
SCALARTYPE p = d[i];
for (int j = i+1; j < n; j++) {
if (d[j] > p) {
k = j;
p = d[j);
}
}
if (k != i) {
d[k] = d[i];
d[i] = p;
for (int j = 0; j < n; j++) {
p = Q(j, i);
Q(j, i) = Q(j, k);
Q(j, k) = p;
}
}
}
*/
}
} // namespace linalg
} // namespace viennacl
#endif
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
matrix.hpp
Implementation of the dense matrix class.
viennacl::traits::size1
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:163
viennacl::linalg::tql2
void tql2(matrix_base< SCALARTYPE, F > &Q, VectorType &d, VectorType &e)
Definition: tql2.hpp:131
viennacl::matrix_base
Definition: matrix_def.hpp:103
s2
viennacl::scalar< int > s2
Definition: global_variables.cpp:58
qr-method-common.hpp
Common routines used for the QR method and SVD. Experimental.
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::vector
Definition: forwards.h:266
viennacl::linalg::detail::pythag
SCALARTYPE pythag(SCALARTYPE a, SCALARTYPE b)
Definition: qr-method-common.hpp:65
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::copy
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) ...
Definition: circulant_matrix.hpp:150
viennacl::linalg::tql1
void tql1(vcl_size_t n, VectorType &d, VectorType &e)
Definition: tql2.hpp:41
viennacl::linalg::givens_next
void givens_next(matrix_base< NumericT > &Q, vector_base< NumericT > &tmp1, vector_base< NumericT > &tmp2, int l, int m)
This function updates the matrix Q. It is part of the tql2 algorithm.
Definition: matrix_operations.hpp:1092