cherkU

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

// Algorithms.
#include <myramath/utility/random.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/herk.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>
#include <myramath/dense/threshold.h>

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

using namespace myra;

namespace {

template<class Precision> void test(int I, int J, Precision tolerance)
  {
  // Useful typedef.
  typedef std::complex<Precision> Number;
  myra::out() << typestring<Number>() << std::endl;
  
  // Make random nonsquare matrix A.
  auto A = Matrix<Number>::random(I,J);
  Precision alpha = random<Precision>();
  Precision beta = random<Precision>();

  // Compare herk_inplace('N') to gemm_inplace(A,'N',A,'H')
    {
    auto C1 = Matrix<Number>::random(I,I);
    for (int i = 0; i < I; ++i)
      C1(i,i) = std::real(C1(i,i));
    Matrix<Number> C2 = C1;
    gemm_inplace(C1,A,'N',A,'H',Number(alpha),Number(beta));
    herk_inplace(C2,'U',A,'N',alpha,beta);
    Precision error = frobenius(triu(C1-C2));
    myra::out() << " |herk('N')-gemm('N','H')| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }

  // Compare herk_inplace('T') to gemm_inplace(A,'T',A,'C')
    {
    auto C1 = Matrix<Number>::random(J,J);
    for (int j = 0; j < J; ++j)
      C1(j,j) = std::real(C1(j,j));
    Matrix<Number> C2 = C1;
    gemm_inplace(C1,A,'T',A,'C',Number(alpha),Number(beta));
    herk_inplace(C2,'U',A,'T',alpha,beta);
    Precision error = frobenius(triu(C1-C2));
    myra::out() << " |herk('T')-gemm('T','C')| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  
  // Compare herk_inplace('H') to gemm_inplace(A,'H',A,'N')
    {
    auto C1 = Matrix<Number>::random(J,J);
    for (int j = 0; j < J; ++j)
      C1(j,j) = std::real(C1(j,j));
    Matrix<Number> C2 = C1;
    gemm_inplace(C1,A,'H',A,'N',Number(alpha),Number(beta));
    herk_inplace(C2,'U',A,'H',alpha,beta);
    Precision error = frobenius(triu(C1-C2));
    myra::out() << " |herk('H')-gemm('H','N')| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  
  // Compare herk_inplace('C') to gemm(A,'C',A,'T')
    {
    auto C1 = Matrix<Number>::random(I,I);
    for (int i = 0; i < I; ++i)
      C1(i,i) = std::real(C1(i,i));
    Matrix<Number> C2 = C1;
    gemm_inplace(C1,A,'C',A,'T',Number(alpha),Number(beta)); 
    herk_inplace(C2,'U',A,'C',alpha,beta);
    Precision error = frobenius(triu(C1-C2));
    myra::out() << " |herk('C')-gemm('C','T')| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
    
  }

} // namespace

ADD_TEST("cherkU","[dense][blas]")
  { test<float>(54,37,1.0e-4f); }
  
ADD_TEST("zherkU","[dense][blas]")
  { test<double>(54,37,1.0e-8); }

Results: [PASS]

std::complex<float>
 |herk('N')-gemm('N','H')| = 5.83482e-07
 |herk('T')-gemm('T','C')| = 5.3e-07
 |herk('H')-gemm('H','N')| = 5.3e-07
 |herk('C')-gemm('C','T')| = 5.83482e-07

Go back to Summary of /test programs.