zher2kLower

Source: tests/dense/her2k3.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/LowerMatrix.h>
#include <myramath/dense/LowerMatrixRange.h>

// Algorithms.
#include <myramath/utility/random.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/her2k.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/hermitian.h>
#include <myramath/dense/conjugate.h>
#include <myramath/dense/tril.h>
#include <myramath/dense/triu.h>

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

using namespace myra;

namespace {

template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  // Useful typedef.
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  // Make random nonsquare matrix A and B.
  auto A = Matrix<Number>::random(I,J);
  auto B = Matrix<Number>::random(I,J);
  Number alpha = random<Number>();
  Precision beta = random<Precision>();
  Number one(1);
  // Compare her2k('N') to gemm(A,'N',B,'H')+gemm(B,'N',A,'H')
    {
    auto C1 = Matrix<Number>::random(I,I);
    C1 = C1 + hermitian(C1);
    LowerMatrix<Number> C2(C1);
    gemm_inplace(C1,A,'N',B,'H',alpha,beta);
    gemm_inplace(C1,B,'N',A,'H',conjugate(alpha),one); 
    her2k_inplace(C2,A,B,'N',alpha,beta);
    Precision error = frobenius(C2.make_Matrix('H')-C1);
    myra::out() << "  |her2k('N')-gemm('N','H')| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Compare her2k('T') to gemm(A,'T',B,'C')+gemm(B,'T',A,'C')
    {
    auto C1 = Matrix<Number>::random(J,J);
    C1 = C1 + hermitian(C1);
    LowerMatrix<Number> C2(C1);
    gemm_inplace(C1,A,'T',B,'C',alpha,beta);
    gemm_inplace(C1,B,'T',A,'C',conjugate(alpha),one);    
    her2k_inplace(C2,A,B,'T',alpha,beta);
    Precision error = frobenius(C2.make_Matrix('H')-C1);
    myra::out() << "  |her2k('T')-gemm('T','C')| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Compare her2k('H') to gemm(A,'H',B,'N')+gemm(B,'H',A,'N')
    {
    auto C1 = Matrix<Number>::random(J,J);
    C1 = C1 + hermitian(C1);
    LowerMatrix<Number> C2(C1);
    gemm_inplace(C1,A,'H',B,'N',alpha,beta);
    gemm_inplace(C1,B,'H',A,'N',conjugate(alpha),one);    
    her2k_inplace(C2,A,B,'H',alpha,beta);
    Precision error = frobenius(C2.make_Matrix('H')-C1);
    myra::out() << "  |her2k('H')-gemm('H','N')| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Compare her2k('C') to gemm(A,'C',B,'T')+gemm(B,'C',A,'T')
    {
    auto C1 = Matrix<Number>::random(I,I);
    C1 = C1 + hermitian(C1);
    LowerMatrix<Number> C2(C1);
    gemm_inplace(C1,A,'C',B,'T',alpha,beta);
    gemm_inplace(C1,B,'C',A,'T',conjugate(alpha),one); 
    her2k_inplace(C2,A,B,'C',alpha,beta);
    Precision error = frobenius(C2.make_Matrix('H')-C1);
    myra::out() << "  |her2k('C')-gemm('C','T')| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

} // namespace

ADD_TEST("sher2kLower","[dense][blas]")
  { test<NumberS>(57,24,1.0e-4f); }

ADD_TEST("dher2kLower","[dense][blas]")
  { test<NumberD>(57,24,1.0e-8); }

ADD_TEST("cher2kLower","[dense][blas]")
  { test<NumberC>(57,24,1.0e-4f); }
  
ADD_TEST("zher2kLower","[dense][blas]")
  { test<NumberZ>(57,24,1.0e-8); }

Results: [PASS]

std::complex<double>
  |her2k('N')-gemm('N','H')| = 6.92643e-14
  |her2k('T')-gemm('T','C')| = 6.32414e-14
  |her2k('H')-gemm('H','N')| = 6.21473e-14
  |her2k('C')-gemm('C','T')| = 7.07046e-14

Go back to Summary of /test programs.