pgemm

Source: tests/pdense/pgemm.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/dense/Matrix.h>
#include <myramath/dense/MatrixRange.h>

// Serial algorithms.
#include <myramath/utility/random.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/frobenius.h>

// Parallel algorithms.
#include <myramath/pdense/pgemm.h>
#include <myramath/pdense/Options.h>

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

using namespace myra;
typedef pdense::Options Options;

namespace {

template<class Number> void test_detail(int I, int J, int K, char op_A, char op_B, typename ReflectPrecision<Number>::type tolerance)
  {
  // Construct random A.
  std::pair<int,int> A_size(I,J);
  if (op_A == 'T' || op_A == 'H')
    std::swap(A_size.first,A_size.second);
  auto A = Matrix<Number>::random(A_size.first,A_size.second);
  // Construct random B.
  std::pair<int,int> B_size(J,K);
  if (op_B == 'T' || op_B == 'H')
    std::swap(B_size.first,B_size.second);  
  auto B = Matrix<Number>::random(B_size.first,B_size.second);
  // Construct random C1 such that C = beta*C + alpha*op(A)*op(B) is well formed.
  auto C1 = Matrix<Number>::random(I,K);
  // Copy into C2.
  Matrix<Number> C2 = C1;
  // Randomize alpha/beta.
  Number alpha = random<Number>();
  Number beta  = random<Number>();
  Options options = Options::create().set_blocksize(64);
  // Update C1/C2 using gemm/pgemm.
  gemm_inplace(C1, A, op_A, B, op_B, alpha, beta);
  pgemm_inplace(C2, A, op_A, B, op_B, alpha, beta, options);
  // Check error between C1 and C2.
  typedef typename ReflectPrecision<Number>::type Precision;
  Precision error = frobenius(C1-C2)/frobenius(C1);
  myra::out() << "  |A^" << op_A << " * B^" << op_B << "|, gemm-pgemm = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

template<class Number> void test(int I, int J, int K,  typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;   
  test_detail<Number>(I,J,K,'N','N',tolerance); // C = A   * B
  test_detail<Number>(I,J,K,'N','T',tolerance); // C = A   * B^T
  test_detail<Number>(I,J,K,'N','H',tolerance); // C = A   * B^H
  test_detail<Number>(I,J,K,'T','N',tolerance); // C = A^T * B
  test_detail<Number>(I,J,K,'T','T',tolerance); // C = A^T * B^T
  test_detail<Number>(I,J,K,'T','H',tolerance); // C = A^T * B^H
  test_detail<Number>(I,J,K,'H','N',tolerance); // C = A^H * B
  test_detail<Number>(I,J,K,'H','T',tolerance); // C = A^H * B^T
  test_detail<Number>(I,J,K,'H','H',tolerance); // C = A^H * B^H
  }

} // namespace

ADD_TEST("pgemm","[pdense][parallel]")
  {
  int I = 207;
  int J = 139;
  int K = 162;
  test<NumberS> (I,J,K,1.0e-4f);
  test<NumberD> (I,J,K,1.0e-9);
  test<NumberC> (I,J,K,1.0e-4f);
  test<NumberZ> (I,J,K,1.0e-9);   
  }

Results: [PASS]

float
  |A^N * B^N|, gemm-pgemm = 1.91876e-07
  |A^N * B^T|, gemm-pgemm = 1.90779e-07
  |A^N * B^H|, gemm-pgemm = 1.91628e-07
  |A^T * B^N|, gemm-pgemm = 6.6212e-08
  |A^T * B^T|, gemm-pgemm = 1.91143e-07
  |A^T * B^H|, gemm-pgemm = 1.89745e-07
  |A^H * B^N|, gemm-pgemm = 1.83065e-07
  |A^H * B^T|, gemm-pgemm = 1.75062e-07
  |A^H * B^H|, gemm-pgemm = 1.92316e-07
double
  |A^N * B^N|, gemm-pgemm = 0
  |A^N * B^T|, gemm-pgemm = 0
  |A^N * B^H|, gemm-pgemm = 0
  |A^T * B^N|, gemm-pgemm = 0
  |A^T * B^T|, gemm-pgemm = 0
  |A^T * B^H|, gemm-pgemm = 0
  |A^H * B^N|, gemm-pgemm = 0
  |A^H * B^T|, gemm-pgemm = 0
  |A^H * B^H|, gemm-pgemm = 0
std::complex<float>
  |A^N * B^N|, gemm-pgemm = 1.94059e-07
  |A^N * B^T|, gemm-pgemm = 1.94295e-07
  |A^N * B^H|, gemm-pgemm = 1.93059e-07
  |A^T * B^N|, gemm-pgemm = 1.92857e-07
  |A^T * B^T|, gemm-pgemm = 1.95841e-07
  |A^T * B^H|, gemm-pgemm = 1.96164e-07
  |A^H * B^N|, gemm-pgemm = 1.93078e-07
  |A^H * B^T|, gemm-pgemm = 1.92529e-07
  |A^H * B^H|, gemm-pgemm = 1.95202e-07
std::complex<double>
  |A^N * B^N|, gemm-pgemm = 0
  |A^N * B^T|, gemm-pgemm = 0
  |A^N * B^H|, gemm-pgemm = 0
  |A^T * B^N|, gemm-pgemm = 0
  |A^T * B^T|, gemm-pgemm = 0
  |A^T * B^H|, gemm-pgemm = 0
  |A^H * B^N|, gemm-pgemm = 0
  |A^H * B^T|, gemm-pgemm = 0
  |A^H * B^H|, gemm-pgemm = 0

Go back to Summary of /test programs.