slqt

Source: tests/dense/lqt.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/lqt.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/trmm.h>
#include <myramath/dense/gelqf.h>
#include <myramath/dense/ormlq.h>
#include <myramath/dense/hermitian.h>
#include <myramath/dense/frobenius.h>

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

using namespace myra;

namespace {

// Tests lqt() 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=LQ, form T in tril(A)
  auto A = Matrix<Number>::random(I,J);
  auto B1 = Matrix<Number>::random(I,J);
  auto B2 = B1;
  auto tau = gelqf_inplace(A);
  lqt_inplace(A,tau);
  // Apply scalar, B1 = B1*(I-tau*v*v'), using ormlq_inplace() ..
  ormlq_inplace(A,tau,B1,'R','N');
  // Apply block reflections, B2 = B2*(I-V'*T*V), using lqt() ..
  // .. partition V and B.
  auto VlVr = A.split_left(K);
  auto Vl = VlVr.first;
  auto Vr = VlVr.second;
  auto BlBr = B2.split_left(K);
  auto Bl = BlBr.first;
  auto Br = BlBr.second;
  // .. form W = B2*V'
  Matrix<Number> W(Bl);
  trmm_inplace('R','U','H',Vl,W,'U');
  gemm_inplace(W,Br,'N',Vr,'H',one,one);
  // .. update W = W*T
  trmm_inplace('R','L','N',Vl,W,'N',one);
  // .. update B -= W*V
  gemm_inplace(Br,W,'N',Vr,'N',-one,one);
  trmm_inplace('R','U','N',Vl,W,'U');
  Bl -= W;
  // Check answers.
  Precision error = frobenius(B1-B2) / frobenius(B1);
  myra::out() << "|Q*B1(ormlq) - Q*B2(lqt)| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

} // namespace

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

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

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

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

Results: [PASS]

float
|Q*B1(ormlq) - Q*B2(lqt)| = 8.18313e-07

Go back to Summary of /test programs.