sparse_gemm

Source: tests/sparse/gemm.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>
#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>

// Algorithms.
#include <myramath/sparse/gemm.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/frobenius.h>

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

using namespace myra;

namespace {

template<class Number> void test(int I, int J, int N, int K, typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  // Make random sparse A operand.
  auto A = SparseMatrix<Number>::random(I,J,N);
  // Walk every possible combination of op_A and op_B
  std::vector<char> ops;
  ops.push_back('N');
  ops.push_back('T');
  ops.push_back('H');
  ops.push_back('C');
  // Tests with the SparseMatrix appearing on the 'L'eft ..
  for (int a = 0; a < ops.size(); ++a)
    for (int b = 0; b < ops.size(); ++b)
      {
      char op_A = ops[a];
      char op_B = ops[b];
      // Deduce size of B operand.
            std::pair<int,int> A_size(I,J);
            if (op_A == 'T' || op_A == 'H')
                std::swap(A_size.first, A_size.second);
            std::pair<int,int> B_size(A_size.second,K);
            if (op_B == 'T' || op_B == 'H')
                std::swap(B_size.first, B_size.second);    
      // Construct random B operand of given size.
      auto B = Matrix<Number>::random(B_size.first,B_size.second);
      // Generate C = op(A)*op(B), sparsely and densely.
      Matrix<Number> C1 = gemm(A,op_A,B,op_B); 
      Matrix<Number> C2 = gemm(A.make_Matrix(),op_A,B,op_B);
      // Check error between the two.
      Precision error = frobenius(C1-C2);
      myra::out() << "  |A^" << op_A << "*B^" << op_B << " [sparse*dense]| = " << error << std::endl;
      REQUIRE(error < tolerance);
      }
  // Tests with the SparseMatrix appearing on the 'R'ight ..
  for (int a = 0; a < ops.size(); ++a)
    for (int b = 0; b < ops.size(); ++b)
      {
      char op_A = ops[a];
      char op_B = ops[b];
      // Deduce size of B operand.
            std::pair<int,int> A_size(I,J);
            if (op_A == 'T' || op_A == 'H')
                std::swap(A_size.first, A_size.second);
            std::pair<int,int> B_size(K,A_size.first);
            if (op_B == 'T' || op_B == 'H')
                std::swap(B_size.first, B_size.second);    
      // Construct random B operand of given size.
      auto B = Matrix<Number>::random(B_size.first,B_size.second);
      // Generate C = op(B)*op(A), sparsely and densely.
      Matrix<Number> C1 = gemm(B,op_B,A,op_A); 
      Matrix<Number> C2 = gemm(B,op_B,A.make_Matrix(),op_A);
      Precision error = frobenius(C1-C2);
      myra::out() << "  |B^" << op_B << "*A^" << op_A << " [dense*sparse]| = " << error << std::endl;
      REQUIRE(error < tolerance);
      }
  }

} // namespace 

ADD_TEST("sparse_gemm","[sparse]")
  {
  test<NumberS>(15,10,40,8,1.0e-4f);
  test<NumberD>(15,10,40,8,1.0e-9);
  test<NumberC>(15,10,40,8,1.0e-4f);
  test<NumberZ>(15,10,40,8,1.0e-9);
  }

Results: [PASS]

float
  |A^N*B^N [sparse*dense]| = 1.94031e-07
  |A^N*B^T [sparse*dense]| = 1.76519e-07
  |A^N*B^H [sparse*dense]| = 2.05409e-07
  |A^N*B^C [sparse*dense]| = 2.078e-07
  |A^T*B^N [sparse*dense]| = 2.5154e-07
  |A^T*B^T [sparse*dense]| = 2.85827e-07
  |A^T*B^H [sparse*dense]| = 3.13821e-07
  |A^T*B^C [sparse*dense]| = 2.78602e-07
  |A^H*B^N [sparse*dense]| = 2.03185e-07
  |A^H*B^T [sparse*dense]| = 2.00622e-07
  |A^H*B^H [sparse*dense]| = 1.87717e-07
  |A^H*B^C [sparse*dense]| = 1.4769e-07
  |A^C*B^N [sparse*dense]| = 2.16122e-07
  |A^C*B^T [sparse*dense]| = 2.56853e-07
  |A^C*B^H [sparse*dense]| = 2.20496e-07
  |A^C*B^C [sparse*dense]| = 2.30649e-07
  |B^N*A^N [dense*sparse]| = 1.92069e-07
  |B^T*A^N [dense*sparse]| = 2.87789e-07
  |B^H*A^N [dense*sparse]| = 2.30528e-07
  |B^C*A^N [dense*sparse]| = 2.24809e-07
  |B^N*A^T [dense*sparse]| = 2.40501e-07
  |B^T*A^T [dense*sparse]| = 1.68636e-07
  |B^H*A^T [dense*sparse]| = 2.17397e-07
  |B^C*A^T [dense*sparse]| = 1.73469e-07
  |B^N*A^H [dense*sparse]| = 1.16892e-07
  |B^T*A^H [dense*sparse]| = 1.54104e-07
  |B^H*A^H [dense*sparse]| = 2.17707e-07
  |B^C*A^H [dense*sparse]| = 2.2582e-07
  |B^N*A^C [dense*sparse]| = 2.25941e-07
  |B^T*A^C [dense*sparse]| = 2.24849e-07
  |B^H*A^C [dense*sparse]| = 2.61129e-07
  |B^C*A^C [dense*sparse]| = 1.58164e-07
double
  |A^N*B^N [sparse*dense]| = 5.1608e-16
  |A^N*B^T [sparse*dense]| = 5.87731e-16
  |A^N*B^H [sparse*dense]| = 4.57628e-16
  |A^N*B^C [sparse*dense]| = 4.42586e-16
  |A^T*B^N [sparse*dense]| = 6.06509e-16
  |A^T*B^T [sparse*dense]| = 5.47404e-16
  |A^T*B^H [sparse*dense]| = 3.94016e-16
  |A^T*B^C [sparse*dense]| = 7.79754e-16
  |A^H*B^N [sparse*dense]| = 4.92296e-16
  |A^H*B^T [sparse*dense]| = 5.27354e-16
  |A^H*B^H [sparse*dense]| = 5.05635e-16
  |A^H*B^C [sparse*dense]| = 4.38635e-16
  |A^C*B^N [sparse*dense]| = 6.92146e-16
  |A^C*B^T [sparse*dense]| = 4.87183e-16
  |A^C*B^H [sparse*dense]| = 3.82395e-16
  |A^C*B^C [sparse*dense]| = 4.13042e-16
  |B^N*A^N [dense*sparse]| = 7.07147e-16
  |B^T*A^N [dense*sparse]| = 6.25456e-16
  |B^H*A^N [dense*sparse]| = 5.16818e-16
  |B^C*A^N [dense*sparse]| = 5.51151e-16
  |B^N*A^T [dense*sparse]| = 4.23514e-16
  |B^T*A^T [dense*sparse]| = 5.86953e-16
  |B^H*A^T [dense*sparse]| = 4.77263e-16
  |B^C*A^T [dense*sparse]| = 5.21753e-16
  |B^N*A^H [dense*sparse]| = 5.2208e-16
  |B^T*A^H [dense*sparse]| = 4.0231e-16
  |B^H*A^H [dense*sparse]| = 6.14986e-16
  |B^C*A^H [dense*sparse]| = 3.18103e-16
  |B^N*A^C [dense*sparse]| = 4.87597e-16
  |B^T*A^C [dense*sparse]| = 2.59653e-16
  |B^H*A^C [dense*sparse]| = 5.05234e-16
  |B^C*A^C [dense*sparse]| = 7.89419e-16
std::complex<float>
  |A^N*B^N [sparse*dense]| = 5.51267e-07
  |A^N*B^T [sparse*dense]| = 6.13317e-07
  |A^N*B^H [sparse*dense]| = 5.57662e-07
  |A^N*B^C [sparse*dense]| = 0
  |A^T*B^N [sparse*dense]| = 6.85747e-07
  |A^T*B^T [sparse*dense]| = 6.28747e-07
  |A^T*B^H [sparse*dense]| = 7.4517e-07
  |A^T*B^C [sparse*dense]| = 6.51574e-07
  |A^H*B^N [sparse*dense]| = 6.52851e-07
  |A^H*B^T [sparse*dense]| = 5.81575e-07
  |A^H*B^H [sparse*dense]| = 5.63049e-07
  |A^H*B^C [sparse*dense]| = 7.19415e-07
  |A^C*B^N [sparse*dense]| = 0
  |A^C*B^T [sparse*dense]| = 5.22738e-07
  |A^C*B^H [sparse*dense]| = 5.46933e-07
  |A^C*B^C [sparse*dense]| = 5.96861e-07
  |B^N*A^N [dense*sparse]| = 6.3997e-07
  |B^T*A^N [dense*sparse]| = 6.33915e-07
  |B^H*A^N [dense*sparse]| = 6.50381e-07
  |B^C*A^N [dense*sparse]| = 0
  |B^N*A^T [dense*sparse]| = 5.54095e-07
  |B^T*A^T [dense*sparse]| = 4.53339e-07
  |B^H*A^T [dense*sparse]| = 5.13361e-07
  |B^C*A^T [dense*sparse]| = 6.15344e-07
  |B^N*A^H [dense*sparse]| = 4.86578e-07
  |B^T*A^H [dense*sparse]| = 5.87265e-07
  |B^H*A^H [dense*sparse]| = 5.96212e-07
  |B^C*A^H [dense*sparse]| = 5.36546e-07
  |B^N*A^C [dense*sparse]| = 0
  |B^T*A^C [dense*sparse]| = 5.52512e-07
  |B^H*A^C [dense*sparse]| = 6.91511e-07
  |B^C*A^C [dense*sparse]| = 6.3388e-07
std::complex<double>
  |A^N*B^N [sparse*dense]| = 1.36183e-15
  |A^N*B^T [sparse*dense]| = 1.38855e-15
  |A^N*B^H [sparse*dense]| = 1.24241e-15
  |A^N*B^C [sparse*dense]| = 0
  |A^T*B^N [sparse*dense]| = 1.40783e-15
  |A^T*B^T [sparse*dense]| = 1.7345e-15
  |A^T*B^H [sparse*dense]| = 1.90833e-15
  |A^T*B^C [sparse*dense]| = 1.4234e-15
  |A^H*B^N [sparse*dense]| = 1.61862e-15
  |A^H*B^T [sparse*dense]| = 1.69424e-15
  |A^H*B^H [sparse*dense]| = 1.3333e-15
  |A^H*B^C [sparse*dense]| = 1.40197e-15
  |A^C*B^N [sparse*dense]| = 0
  |A^C*B^T [sparse*dense]| = 1.28503e-15
  |A^C*B^H [sparse*dense]| = 1.50084e-15
  |A^C*B^C [sparse*dense]| = 1.47281e-15
  |B^N*A^N [dense*sparse]| = 1.82692e-15
  |B^T*A^N [dense*sparse]| = 1.41146e-15
  |B^H*A^N [dense*sparse]| = 1.66941e-15
  |B^C*A^N [dense*sparse]| = 0
  |B^N*A^T [dense*sparse]| = 1.30283e-15
  |B^T*A^T [dense*sparse]| = 1.42823e-15
  |B^H*A^T [dense*sparse]| = 1.47206e-15
  |B^C*A^T [dense*sparse]| = 1.22227e-15
  |B^N*A^H [dense*sparse]| = 1.22042e-15
  |B^T*A^H [dense*sparse]| = 1.31776e-15
  |B^H*A^H [dense*sparse]| = 1.31069e-15
  |B^C*A^H [dense*sparse]| = 1.33438e-15
  |B^N*A^C [dense*sparse]| = 0
  |B^T*A^C [dense*sparse]| = 1.47503e-15
  |B^H*A^C [dense*sparse]| = 1.49687e-15
  |B^C*A^C [dense*sparse]| = 1.48648e-15

Go back to Summary of /test programs.