dstev

Source: tests/dense/stev.cpp

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

// Containers.
#include <myramath/utility/Number.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/Vector.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/DiagonalMatrix.h>
#include <myramath/dense/DiagonalMatrixRange.h>

// Algorithms.
#include <myramath/dense/expr.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/dimm.h>
#include <myramath/dense/stev.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/frobenius.h>

// Reporting.
#include <tests/myratest.h>

using namespace myra;

namespace {

// Implementation detail of all tests.
template<class Precision> void test_detail(Vector<Precision>& T0, Vector<Precision>& T1, Precision tolerance)
  {
  myra::out() << typestring<Precision>() << std::endl;
  // Form T.
  int N = T0.size();
  Matrix<Precision> T(N,N);
  for (int n = 0; n < N; ++n)
    T(n,n) = T0(n);
  for (int n = 0; n < N-1; ++n)
    T(n+1,n) = T(n,n+1) = T1(n);
  // Decompose T = X*D*X'
  auto XD = stev(T0,T1);
  const auto& X = XD.first;
  const auto& D = XD.second;
  auto I = Matrix<Precision>::identity(N);
  // Check residual error |R|, where R = T-X*D*X'
  Precision R_error = frobenius(T-X*D*transpose(X))/frobenius(D);
  myra::out() << "  |T-X*D*X'| = " << R_error << std::endl;
  // Check eigenvector orthogonality, X'X = I.
  Precision X_error = frobenius(gemm(X,'T',X)-I);
  myra::out() << "  |X'X-I| = " << X_error << std::endl;
  // Check the layout assumptions of X (orthogonal columns)
  int K = N/2;
  auto Xk = X.left(K);
  auto Ik = Matrix<Precision>::identity(K);  
  Precision Xk_error = frobenius(gemm(Xk,'T',Xk)-Ik)/N;
  myra::out() << "  |Xk'Xk-Ik| = " << Xk_error << std::endl;  
  // Verify eigenvalues are sorted in ascending order.
  bool sorted = true;
  for (int n = 0; n < N-1; ++n)
    sorted &= D(n) <= D(n+1);
  myra::out() << std::boolalpha;
  myra::out() << "  sorted D(n) <= D(n+1) = " << sorted << std::endl;
  // Apply tests.
  REQUIRE(R_error < tolerance);
  REQUIRE(X_error < tolerance);
  REQUIRE(Xk_error < tolerance);
  REQUIRE(sorted);
  }

// Tests a random tridiagonal matrix.
template<class Precision> void test_random(int N, Precision tolerance)
  {
  auto T0 = Vector<Precision>::random(N);
  auto T1 = Vector<Precision>::random(N-1);
  test_detail(T0,T1,tolerance);
  }

// Tests a tridiagonal matrix with zero off-diagonal (that is, a zero within T1 vector)
template<class Precision> void test_split(int N, int S, Precision tolerance)
  {
  auto T0 = Vector<Precision>::random(N);
  auto T1 = Vector<Precision>::random(N-1);
  T1(S) = 0;
  test_detail(T0,T1,tolerance);
  }

} // namespace

// Random test, float precision.
ADD_TEST("sstev","[dense][lapack]")
  { test_random<float>(50,1.0e-4f); }

// Zero off-diagonal test, float precision.
ADD_TEST("sstev_split","[dense][lapack]")
  { test_split<float>(100,27,1.0e-4f); }

// Stress test for float precision, lots of sizes (N) and lots of trials (T).
ADD_TEST("sstev_stress","[dense][lapack][.]")
  {
  int N = 100;
  int T = 10;
  for (int n = 1; n < N; ++n)
    for (int t = 0; t < T; ++t)
      test_random<float>(n,1.0e-4f);
  }

// Random test, double precision.
ADD_TEST("dstev","[dense][lapack]")
  { test_random<double>(50,1.0e-12); }

// Zero off-diagonal test, double precision.
ADD_TEST("dstev_split","[dense][lapack]")
  { test_split<double>(100,27,1.0e-12); }

// Stress test for double precision, lots of sizes (N) and lots of trials (T).
ADD_TEST("dstev_stress","[dense][lapack][.]")
  {
  int N = 100;
  int T = 10;
  for (int n = 1; n < N; ++n)
    for (int t = 0; t < T; ++t)
      test_random<double>(n,1.0e-12);
  }

Results: [PASS]

double
  |T-X*D*X'| = 4.3845e-15
  |X'X-I| = 1.22716e-14
  |Xk'Xk-Ik| = 1.56473e-16
  sorted D(n) <= D(n+1) = true

Go back to Summary of /test programs.