dqrt

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

// Algorithms.
#include <myramath/dense/qrt.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/trmm.h>
#include <myramath/dense/geqrf.h>
#include <myramath/dense/ormqr.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/hermitian.h>

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

using namespace myra;

namespace {

// Tests qrt() on a tall/skinny A.
template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  // Useful typedefs/constants.
  typedef typename ReflectPrecision<Number>::type Precision;
  Number one(1);
  int K = I < J ? I : J;
  myra::out() << typestring<Number>() << std::endl;
  // Form random A/B, factor A=QR, form T in triu(A)
  auto A = Matrix<Number>::random(I,J);
  auto B1 = Matrix<Number>::random(I,J);
  auto B2 = B1;
  auto tau = geqrf_inplace(A);
  qrt_inplace(A,tau);
  // Apply scalar reflections, B1 = (I-tau*v*v')*B1, using ormqr_inplace() ..
  ormqr_inplace(A,tau,B1,'L','N');
  // Apply block reflection, B2 = (I-V*T*V')*B2, using qrt() ..
  // .. partition V and B.
  auto VtVb = A.split_top(K);
  auto Vt = VtVb.first;
  auto Vb = VtVb.second;
  auto BtBb = B2.split_top(K);
  auto Bt = BtBb.first;
  auto Bb = BtBb.second;
  // .. form W = V'*B
  Matrix<Number> W(Bt);
  trmm_inplace('L','L','H',Vt,W,'U');
  gemm_inplace(W,Vb,'H',Bb,'N',one,one);
  // .. update W = T*W
  trmm_inplace('L','U','N',Vt,W,'N',one);
  // .. update B -= V*W
  gemm_inplace(Bb,Vb,'N',W,'N',-one,one);
  trmm_inplace('L','L','N',Vt,W,'U');
  Bt -= W;
  // Check answers.
  Precision error = frobenius(B1-B2) / frobenius(B1);
  myra::out() << "|Q*B1(ormqr) - Q*B2(qrt)| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

} // namespace

ADD_TEST("sqrt","[dense][lapack]")
  { test<NumberS>(178,159,1.0e-4f); }

ADD_TEST("dqrt","[dense][lapack]")
  { test<NumberD>(178,159,1.0e-10); }

ADD_TEST("cqrt","[dense][lapack]")
  { test<NumberC>(178,159,1.0e-4f); }

ADD_TEST("zqrt","[dense][lapack]")
  { test<NumberZ>(178,159,1.0e-10); }

Results: [PASS]

double
|Q*B1(ormqr) - Q*B2(qrt)| = 1.36386e-15

Go back to Summary of /test programs.