zheev

Source: tests/dense/heev.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/MatrixRange.h>
#include <myramath/dense/DiagonalMatrix.h>
#include <myramath/dense/DiagonalMatrixRange.h>

// Algorithms.
#include <myramath/dense/dimm.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/heev.h>
#include <myramath/dense/hermitian.h>
#include <myramath/dense/frobenius.h>

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

using namespace myra;

namespace {

template<class Precision> void test(int N, Precision tolerance)
  {
  typedef typename ReflectComplex<Precision>::type Number;  
  myra::out() << typestring<Number>() << std::endl;
  // Make random hermitian matrix A.
  auto A = Matrix<Number>::random(N,N);
  A += hermitian(A);
  // Eigendecompose A = X*D*X'
  auto XD = heev(A);
  const auto& X = XD.first;
  const auto& D = XD.second;
  auto I = Matrix<Number>::identity(N);
  // Check residual error |R|, where R = A-X*D*X'
  Precision R_error = frobenius(A-X*D*hermitian(X))/frobenius(D);
  myra::out() << "  |A-X*D*X'| = " << R_error << std::endl; 
  // Check eigenvector orthogonality, X'*X = I
  Precision X_error = frobenius(gemm(X,'H',X)-I)/N;
  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<Number>::identity(K);  
  Precision Xk_error = frobenius(gemm(Xk,'H',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);
  }

} // namespace

// Single test, float precision.
ADD_TEST("cheev","[dense][lapack]")
  { test<float>(50,2.0e-5f); }

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

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

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

Results: [PASS]

std::complex<double>
  |A-X*D*X'| = 5.67864e-15
  |X'*X -I| = 2.8205e-16
  |Xk'Xk-Ik| = 1.6696e-16
  sorted D(n) <= D(n+1) = true

Go back to Summary of /test programs.