ptrmm2

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

// Serial algorithms.
#include <myramath/utility/random.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/tril.h>
#include <myramath/dense/triu.h>
#include <myramath/dense/trmm.h>

// Parallel algorithms.
#include <myramath/pdense/ptrmm.h>
#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(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  // Make random L
  auto L = LowerMatrix<Number>::random(I);
  Number one(1);
  auto alpha = random<Number>();
  // Initialize options.
  auto options = Options::create().set_nthreads(4).set_blocksize(64);
  // Try all variations of side/op.
  std::vector<char> sides;
  sides.push_back('L');
  sides.push_back('R');
  std::vector<char> ops;
  ops.push_back('N');
  ops.push_back('T');
  ops.push_back('H');
  ops.push_back('C');
  for (int s = 0; s < sides.size(); ++s)
    for (int o = 0; o < ops.size(); ++o)
      {
      char side = sides[s];
      char op = ops[o];
      // Generate random X/B
      auto X = (side == 'L') ? Matrix<Number>::random(I,J) : Matrix<Number>::random(J,I);
      auto B = (side == 'L') ? alpha*pgemm(L.make_Matrix(),op,X) : alpha*pgemm(X,L.make_Matrix(),op);
      // Multiply using both trmm() and ptrmm() 
      auto X_serial = X;
      auto X_parallel = X;
      trmm_inplace(side,op,L,X_serial,'N',alpha);
      ptrmm_inplace(side,op,L,X_parallel,'N',alpha,options);
      // Compare errors.
      Precision e_serial = frobenius(X_serial-B) / frobenius(B);
      Precision e_parallel = frobenius(X_parallel-B) / frobenius(B);
      if (side == 'L')
        {
        myra::out() << "  |L^" << op << "*B-X| serial = " << e_serial << std::endl;
        myra::out() << "  |L^" << op << "*B-X| parallel = " << e_parallel << std::endl;
        }
      else if (side == 'R')
        {
        myra::out() << "  |B*L^" << op << "-X| serial = " << e_serial << std::endl;
        myra::out() << "  |B*L^" << op << "-X| parallel = " << e_parallel << std::endl;
        }
      REQUIRE(e_serial < tolerance);
      REQUIRE(e_parallel < tolerance);
      }
  }

} // namespace

ADD_TEST("ptrmm2","[pdense][parallel]")
  {
  int I = 256;
  int J = 128;
  test<NumberS>(I,J,1.0e-4f);
  test<NumberD>(I,J,1.0e-8);
  test<NumberC>(I,J,1.0e-4f);
  test<NumberZ>(I,J,1.0e-8);
  }
  

Results: [PASS]

float
  |L^N*B-X| serial = 2.53095e-07
  |L^N*B-X| parallel = 2.18184e-07
  |L^T*B-X| serial = 2.56949e-07
  |L^T*B-X| parallel = 2.02036e-07
  |L^H*B-X| serial = 2.5948e-07
  |L^H*B-X| parallel = 2.02457e-07
  |L^C*B-X| serial = 2.57776e-07
  |L^C*B-X| parallel = 2.19654e-07
  |B*L^N-X| serial = 2.6558e-07
  |B*L^N-X| parallel = 1.98404e-07
  |B*L^T-X| serial = 2.65203e-07
  |B*L^T-X| parallel = 2.19372e-07
  |B*L^H-X| serial = 2.70392e-07
  |B*L^H-X| parallel = 2.21607e-07
  |B*L^C-X| serial = 2.66407e-07
  |B*L^C-X| parallel = 1.97841e-07
double
  |L^N*B-X| serial = 3.87107e-16
  |L^N*B-X| parallel = 2.92475e-16
  |L^T*B-X| serial = 3.81256e-16
  |L^T*B-X| parallel = 2.81406e-16
  |L^H*B-X| serial = 3.85273e-16
  |L^H*B-X| parallel = 2.82781e-16
  |L^C*B-X| serial = 3.84984e-16
  |L^C*B-X| parallel = 2.9615e-16
  |B*L^N-X| serial = 4.04396e-16
  |B*L^N-X| parallel = 2.66533e-16
  |B*L^T-X| serial = 4.07734e-16
  |B*L^T-X| parallel = 2.91352e-16
  |B*L^H-X| serial = 4.02137e-16
  |B*L^H-X| parallel = 2.87477e-16
  |B*L^C-X| serial = 4.00793e-16
  |B*L^C-X| parallel = 2.63748e-16
std::complex<float>
  |L^N*B-X| serial = 2.6978e-07
  |L^N*B-X| parallel = 2.2849e-07
  |L^T*B-X| serial = 2.66903e-07
  |L^T*B-X| parallel = 2.15101e-07
  |L^H*B-X| serial = 2.70349e-07
  |L^H*B-X| parallel = 2.15684e-07
  |L^C*B-X| serial = 2.70909e-07
  |L^C*B-X| parallel = 2.29597e-07
  |B*L^N-X| serial = 2.74094e-07
  |B*L^N-X| parallel = 2.18348e-07
  |B*L^T-X| serial = 2.76707e-07
  |B*L^T-X| parallel = 2.30205e-07
  |B*L^H-X| serial = 2.74062e-07
  |B*L^H-X| parallel = 2.28124e-07
  |B*L^C-X| serial = 2.72584e-07
  |B*L^C-X| parallel = 2.17341e-07
std::complex<double>
  |L^N*B-X| serial = 4.1294e-16
  |L^N*B-X| parallel = 3.16149e-16
  |L^T*B-X| serial = 4.09071e-16
  |L^T*B-X| parallel = 3.10439e-16
  |L^H*B-X| serial = 4.11319e-16
  |L^H*B-X| parallel = 3.10745e-16
  |L^C*B-X| serial = 4.12567e-16
  |L^C*B-X| parallel = 3.15023e-16
  |B*L^N-X| serial = 4.23415e-16
  |B*L^N-X| parallel = 3.17058e-16
  |B*L^T-X| serial = 4.20723e-16
  |B*L^T-X| parallel = 3.19332e-16
  |B*L^H-X| serial = 4.2498e-16
  |B*L^H-X| parallel = 3.21324e-16
  |B*L^C-X| serial = 4.17804e-16
  |B*L^C-X| parallel = 3.16839e-16

Go back to Summary of /test programs.