sparse_gemv

Source: tests/sparse/gemv.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/Vector.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>

// Algorithms.
#include <myramath/sparse/gemv.h>
#include <myramath/dense/gemv.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, typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  Number one(1);
  Number zero(0);
  // Make random sparse A operand.
  auto A = SparseMatrix<Number>::random(I,J,N);
  // Walk every possible combination of op.
  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 o = 0; o < ops.size(); ++o)
    {
    char op = ops[0];
    // Deduce size of b operand.
        std::pair<int,int> A_size(I,J);
        if (op == 'T' || op == 'H')
            std::swap(A_size.first, A_size.second);
    int b_size = A_size.second;
    // Construct random b operand of given size.
    auto b = Vector<Number>::random(b_size);
    // Generate c = op(A)*b, sparsely and densely.
    Vector<Number> c1 = gemv(A,op,b);    
    Vector<Number> c2 = gemv(A.make_Matrix(),op,b);
    // Check error between the two.
    Precision error = frobenius(c1-c2);
    myra::out() << "  |A^" << op << "*b [sparse*vector]| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Tests with the SparseMatrix appearing on the 'R'ight ..
  for (int o = 0; o < ops.size(); ++o)
    {
    char op = ops[0];
    // Deduce size of b operand.
      std::pair<int,int> A_size(I,J);
      if (op == 'T' || op == 'H')
          std::swap(A_size.first, A_size.second);
        int b_size = A_size.first;
    // Construct random b operand of given size.
    auto b = Vector<Number>::random(b_size);
    // Generate C = b*op(A), sparsely and densely.
    Vector<Number> c1 = gemv(b,A,op); 
    Vector<Number> c2 = gemv(b,A.make_Matrix(),op);
    Precision error = frobenius(c1-c2);
    myra::out() << "  |b*A^" << op << " [vector*sparse]| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

} // namespace

ADD_TEST("sparse_gemv","[sparse]")
  {
  test< float >(15,10,40,1.0e-4f);
  test< double >(15,10,40,1.0e-9);
  test< std::complex<float > >(15,10,40,1.0e-4f);
  test< std::complex<double> >(15,10,40,1.0e-9);
  }

Results: [PASS]

float
  |A^N*b [sparse*vector]| = 6.83872e-08
  |A^N*b [sparse*vector]| = 3.72529e-09
  |A^N*b [sparse*vector]| = 1.46191e-07
  |A^N*b [sparse*vector]| = 1.62058e-07
  |b*A^N [vector*sparse]| = 1.8193e-07
  |b*A^N [vector*sparse]| = 8.59293e-08
  |b*A^N [vector*sparse]| = 1.06354e-07
  |b*A^N [vector*sparse]| = 1.35149e-07
double
  |A^N*b [sparse*vector]| = 2.75118e-16
  |A^N*b [sparse*vector]| = 1.33111e-16
  |A^N*b [sparse*vector]| = 7.85046e-17
  |A^N*b [sparse*vector]| = 2.77556e-17
  |b*A^N [vector*sparse]| = 0
  |b*A^N [vector*sparse]| = 2.31909e-16
  |b*A^N [vector*sparse]| = 2.22478e-16
  |b*A^N [vector*sparse]| = 1.30923e-16
std::complex<float>
  |A^N*b [sparse*vector]| = 2.76777e-07
  |A^N*b [sparse*vector]| = 2.43943e-07
  |A^N*b [sparse*vector]| = 2.60877e-07
  |A^N*b [sparse*vector]| = 1.51963e-07
  |b*A^N [vector*sparse]| = 2.93897e-07
  |b*A^N [vector*sparse]| = 3.14693e-07
  |b*A^N [vector*sparse]| = 2.28552e-07
  |b*A^N [vector*sparse]| = 2.22022e-07
std::complex<double>
  |A^N*b [sparse*vector]| = 5.69254e-16
  |A^N*b [sparse*vector]| = 3.67434e-16
  |A^N*b [sparse*vector]| = 3.98367e-16
  |A^N*b [sparse*vector]| = 4.78331e-16
  |b*A^N [vector*sparse]| = 4.91046e-16
  |b*A^N [vector*sparse]| = 2.8441e-16
  |b*A^N [vector*sparse]| = 3.49503e-16
  |b*A^N [vector*sparse]| = 5.2461e-16

Go back to Summary of /test programs.