gebrd_square

Source: tests/dense/gebrd.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>

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/ormqr.h>
#include <myramath/dense/ormlq.h>
#include <myramath/dense/gebrd.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/threshold.h>

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

using namespace myra;

namespace {

// Tests tall/skinny gebrd(). Also handles the square case (see fine print on gebrd).
template<class Number> void test1(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  // Useful constants.
  typedef typename ReflectPrecision<Number>::type Precision;
  int K = I < J ? I : J;
  // Make random A, and a few copy.  
  auto A = Matrix<Number>::random(I,J);
  auto B = A;
  // Factor A = Q*B*P
  auto tau = gebrd_inplace(A);
  auto tauq = tau.first;
  auto taup = tau.second;
  // Multiply out (Q'*A)*P', order matters.
  auto Q = A.range();
  ormqr_inplace(Q,tauq,B,'L','H');
  auto P = A.cut_left(1).top(J);
  ormlq_inplace(P,taup.cut_last(1),B.cut_left(1).top(J),'R','H'); 
  // Subtract out D and E from B.
  for (int k = 0; k < K; ++k)
    B(k,k) -= A(k,k);
  for (int k = 0; k < K-1; ++k)    
    B(k,k+1) -= A(k,k+1);
  // Verify Q'*A*P' = B.
  Precision error = frobenius(B);  
  myra::out() << "|Q'*A*P' - B| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

// Tests short/fat gebrd()
template<class Number> void test2(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  // Useful constants.
  typedef typename ReflectPrecision<Number>::type Precision;
  int K = I < J ? I : J;
  // Make random A, and a few copy.  
  auto A = Matrix<Number>::random(I,J);
  auto B = A;
  // Factor A = Q*B*P
  auto tau = gebrd_inplace(A);
  auto tauq = tau.first;
  auto taup = tau.second;
  // Multiply out Q'*(A*P'), order matters.
  auto P = A.range();
  ormlq_inplace(P,taup,B,'R','H');
  auto Q = A.cut_top(1).left(I);
  ormqr_inplace(Q,tauq.cut_last(1),B.cut_top(1).left(I),'L','H'); 
  // Subtract out D and E from B.
  for (int k = 0; k < K; ++k)
    B(k,k) -= A(k,k);
  for (int k = 0; k < K-1; ++k)    
    B(k+1,k) -= A(k+1,k);  
  // Verify Q'*A*P' = B.
  Precision error = frobenius(B);  
  myra::out() << "|Q'*A*P' - B| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

} // namespace

ADD_TEST("gebrd_tallskinny","[dense][lapack]")
  {
  int N1 = 17;
  int N2 = 13;
  test1<NumberS>(N1,N2,1.0e-5f);
  test1<NumberD>(N1,N2,1.0e-10);
  test1<NumberC>(N1,N2,1.0e-5f);
  test1<NumberZ>(N1,N2,1.0e-10);
  }

ADD_TEST("gebrd_shortfat","[dense][lapack]")
  {
  int N1 = 17;
  int N2 = 13;
  test2<NumberS>(N2,N1,1.0e-5f);
  test2<NumberD>(N2,N1,1.0e-10);
  test2<NumberC>(N2,N1,1.0e-5f);
  test2<NumberZ>(N2,N1,1.0e-10);
  }

ADD_TEST("gebrd_square","[dense][lapack]")
  {
  int N1 = 17;
  int N2 = 13;
  test1<NumberS>(N1,N1,1.0e-5f);
  test1<NumberD>(N1,N1,1.0e-10);
  test1<NumberC>(N1,N1,1.0e-5f);
  test1<NumberZ>(N1,N1,1.0e-10);
  }

Results: [PASS]

|Q'*A*P' - B| = 5.33681e-06
|Q'*A*P' - B| = 8.38749e-15
|Q'*A*P' - B| = 9.12197e-06
|Q'*A*P' - B| = 1.32992e-14

Go back to Summary of /test programs.