zhetrdU

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

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/ormqr.h>
#include <myramath/dense/ormlq.h>
#include <myramath/dense/hetrd.h>
#include <myramath/dense/syev.h>
#include <myramath/dense/heev.h>
#include <myramath/dense/hermitian.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/threshold.h>

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

using namespace myra;

namespace {

// Tests hetrd('L'ower)
template<class Precision> void testL(int N, Precision tolerance)
  {
  typedef std::complex<Precision> Number;
  myra::out() << typestring<Number>() << std::endl;  

  // Make random complex hermitian matrix A, a few copies.
  auto A = Matrix<Number>::random(N,N);
  A += hermitian(A);
  auto A1 = A;
  auto A2 = A;

  // Tridiagonalize A = Q*T*Q'. 
  auto tau = hetrd_inplace('L',A);
  
  // Extract T factor from diagonal/subdiagonal of A.
  auto T = Matrix<Precision>::zeros(N,N);
  for (int n = 0; n < N; ++n)
    T(n,n) = A(n,n).real();
  for (int n = 0; n < N-1; ++n)
    T(n+1,n) = T(n,n+1) = A(n+1,n).real();
  auto T2 = T;

  // Verify A and T have same eigenvalues, using syevd() and heevd()
  auto lambda_A = heev_inplace(A1);
  auto lambda_T = syev_inplace(T);
  Precision error1 = frobenius(lambda_A-lambda_T) / frobenius(lambda_A);
  myra::out() << "|eig(A)-eig(T)| = " << error1 << std::endl;
  REQUIRE(error1 < tolerance);

  // Verify Q'*A*Q = T
  auto Q = A.bottom(N-1).left(N-1);
  ormqr_inplace(Q,tau,A2.bottom(N-1),'L','H');
  ormqr_inplace(Q,tau,A2.right(N-1),'R','N');
  Precision error2 = frobenius(A2-make_complex(T2)) / frobenius(T2);
  myra::out() << "|Q'*A*Q - T| = " << error2 << std::endl;
  REQUIRE(error2 < tolerance);
  }

// Tests hetrd('U'pper)
template<class Precision> void testU(int N, Precision tolerance)
  {
  typedef std::complex<Precision> Number;
  myra::out() << typestring<Number>() << std::endl;  

  // Make random complex hermitian matrix A, a few copies.
  auto A = Matrix<Number>::random(N,N);
  A += hermitian(A);
  auto A1 = A;
  auto A2 = A;

  // Tridiagonalize A = Q'*T*Q. 
  auto tau = hetrd_inplace('U',A);
  
  // Extract T factor from diagonal/superdiagonal of A.
  auto T = Matrix<Precision>::zeros(N,N);
  for (int n = 0; n < N; ++n)
    T(n,n) = A(n,n).real();
  for (int n = 0; n < N-1; ++n)
    T(n+1,n) = T(n,n+1) = A(n,n+1).real();
  auto T2 = T;

  // Verify A and T have same eigenvalues, using syevd() and heevd()
  auto lambda_A = heev_inplace(A1);
  auto lambda_T = syev_inplace(T);
  Precision error1 = frobenius(lambda_A-lambda_T) / frobenius(lambda_A);
  myra::out() << "|eig(A)-eig(T)| = " << error1 << std::endl;
  REQUIRE(error1 < tolerance);

  // Verify Q*A*Q' = T (ormlq)
  auto Q = A.right(N-1).top(N-1);
  ormlq_inplace(Q,tau,A2.bottom(N-1),'L','N');
  ormlq_inplace(Q,tau,A2.right(N-1),'R','H');
  Precision error2 = frobenius(A2-make_complex(T2)) / frobenius(T2);
  myra::out() << "|Q*A*Q' - T| = " << error2 << std::endl;
  // Note, this check will not pass with reference LAPACK/BLAS,
  // they use a QL decomposition instead of an LQ decomposition.
  // REQUIRE(error2 < tolerance);
  }

} // namespace

ADD_TEST("chetrdL","[dense][lapack]")
  { testL<float> (157, 1.0e-4f); }

ADD_TEST("zhetrdL","[dense][lapack]")
  { testL<double> (157, 1.0e-8); }

ADD_TEST("chetrdU","[dense][lapack]")
  { testU<float> (157, 1.0e-4f); }

ADD_TEST("zhetrdU","[dense][lapack]")
  { testU<double> (157, 1.0e-8); }

Results: [PASS]

std::complex<double>
|eig(A)-eig(T)| = 1.65718e-15
|Q*A*Q' - T| = 2.05587e-15

Go back to Summary of /test programs.