dsytrdL

Source: tests/dense/sytrd.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/diag.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/ormqr.h>
#include <myramath/dense/ormlq.h>
#include <myramath/dense/sytrd.h>
#include <myramath/dense/syev.h>
#include <myramath/dense/stev.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/threshold.h>

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

using namespace myra;

namespace {

// Tests sytrd('L'ower)
template<class Precision> void testL(int N, Precision tolerance)
  {
  myra::out() << typestring<Precision>() << std::endl;
  // Make random real symmetric matrix A, a few copies.
  auto A = Matrix<Precision>::random(N,N);
  A += transpose(A);
  auto A1 = A;
  auto A2 = A;
  // Tridiagonalize A = Q*T*Q'. 
  auto tau = sytrd_inplace('L',A);
  // Extract T factor from diagonal/subdiagonal of A.
  auto T = Matrix<Precision>::zeros(N,N);
  auto T0 = Vector<Precision>::zeros(N);
  auto T1 = Vector<Precision>::zeros(N-1);  
  for (int n = 0; n < N; ++n)
    T0(n) = T(n,n) = A(n,n);
  for (int n = 0; n < N-1; ++n)
    T1(n) = T(n+1,n) = T(n,n+1) = A(n+1,n);
  // Verify A and T have same eigenvalues, using syev() and stev()
  auto lambda_A = syev_inplace(A1);
  auto lambda_T = stev(T0,T1).second;
  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 (ormqr)
  auto Q = A.bottom(N-1).left(N-1);
  ormqr_inplace(Q,tau,A2.bottom(N-1),'L','T');
  ormqr_inplace(Q,tau,A2.right(N-1),'R','N');
  Precision error2 = frobenius(A2-T) / frobenius(T);
  myra::out() << "|Q'*A*Q - T| = " << error2 << std::endl;
  REQUIRE(error2 < tolerance);
  }

// Tests sytrd('U'pper)
template<class Precision> void testU(int N, Precision tolerance)
  {
  myra::out() << typestring<Precision>() << std::endl;
  // Make random real symmetric matrix A, a few copies.
  auto A = Matrix<Precision>::random(N,N);
  A += transpose(A);
  auto A1 = A;
  auto A2 = A;
  // Tridiagonalize A = Q'*T*Q. 
  auto tau = sytrd_inplace('U',A);
  // Extract T factor from diagonal/superdiagonal of A.
  auto T = Matrix<Precision>::zeros(N,N);
  auto T0 = Vector<Precision>::zeros(N);
  auto T1 = Vector<Precision>::zeros(N-1);  
  for (int n = 0; n < N; ++n)
    T0(n) = T(n,n) = A(n,n);
  for (int n = 0; n < N-1; ++n)
    T1(n) = T(n+1,n) = T(n,n+1) = A(n,n+1);
  // Verify A and T have same eigenvalues, using syev() and stev()
  auto lambda_A = syev_inplace(A1);
  auto lambda_T = stev(T0,T1).second;
  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','T');
  Precision error2 = frobenius(A2-T) / frobenius(T);
  myra::out() << "|Q*A*Q' - T| = " << error2 << std::endl;
  // Note this check will fail with reference LAPACK/BLAS, since
  // they use a QL decomposition instead of an LQ decomposition.
  // REQUIRE(error2 < tolerance);
  }

} // namespace

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

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

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

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

Results: [PASS]

double
|eig(A)-eig(T)| = 0
|Q'*A*Q - T| = 1.48196e-15

Go back to Summary of /test programs.