MyraMath
phemm2


Source: tests/pdense/phemm2.cpp

1 // ========================================================================= //
2 // This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton //
3 // and distributed by MyraCore, LLC. See LICENSE.txt for license terms. //
4 // ========================================================================= //
5 
11 // Containers.
13 #include <myramath/dense/Matrix.h>
17 
18 // Serial algorithms.
21 #include <myramath/dense/tril.h>
22 #include <myramath/dense/triu.h>
23 #include <myramath/dense/hemm.h>
25 
26 // Parallel algorithms.
27 #include <myramath/pdense/phemm.h>
29 
30 // Reporting.
31 #include <tests/myratest.h>
32 
33 using namespace myra;
34 typedef pdense::Options Options;
35 
36 namespace {
37 
38 template<class Precision> void test(int I, int J, Precision tolerance)
39  {
40  typedef std::complex<Precision> Number;
41  myra::out() << typestring<Number>() << std::endl;
42  // Make random hermitian matrix A.
43  auto A = Matrix<Number>::random(I,I);
44  A = A + hermitian(A);
45  Number alpha = random<Number>();
46  Number beta = random<Number>();
47  // Grab A's lower triangle into L.
49  // Initialize options.
50  auto options = Options::create().set_nthreads(4).set_blocksize(128);
51  // Check phemm('L'eft)
52  {
53  auto B = Matrix<Number>::random(I,J);
54  auto C1 = Matrix<Number>::random(I,J);
55  auto C2 = C1;
56  hemm_inplace('L', C1, L, B, alpha, beta);
57  phemm_inplace('L', C2, L, B, alpha, beta, options);
58  Precision error = frobenius(C1-C2)/frobenius(C1);
59  myra::out() << " |hemm('L')-phemm('L')| = " << error << std::endl;
60  REQUIRE(error < tolerance);
61  }
62  // Check phemm('R'ight)
63  {
64  auto B = Matrix<Number>::random(J,I);
65  auto C1 = Matrix<Number>::random(J,I);
66  auto C2 = C1;
67  hemm_inplace('R', C1, L, B, alpha, beta);
68  phemm_inplace('R', C2, L, B, alpha, beta, options);
69  Precision error = frobenius(C1-C2)/frobenius(C1);
70  myra::out() << " |hemm('R')-phemm('R')| = " << error << std::endl;
71  REQUIRE(error < tolerance);
72  }
73  }
74 
75 } // namespace
76 
77 ADD_TEST("phemm2","[pdense][parallel]")
78  {
79  int I = 512;
80  int J = 256;
81  test<float > (I,J,1.0e-4f);
82  test<double> (I,J,1.0e-8);
83  }
Interface class for representing subranges of dense Matrix&#39;s.
Thread-parallel version of dense/phemm.h, hermitian Matrix * dense Matrix multiplication.
Routines for computing Frobenius norms of various algebraic containers.
static Matrix< Number > random(int I, int J)
Generates a random Matrix of specified size.
Definition: Matrix.cpp:353
Options pack for routines in /pdense.
Definition: Options.h:24
Routines for hermitian Matrix * dense Matrix multiplication.
Range construct for a lower triangular matrix stored in rectangular packed format.
Definition: syntax.dox:1
Specialized container for a lower triangular matrix, O(N^2/2) storage. Used by symmetry exploiting ma...
Returns the lower triangle of a dense Matrix.
Returns the upper triangle of a dense Matrix.
Various utility functions/classes related to scalar Number types.
General purpose dense matrix container, O(i*j) storage.
Options pack for routines in /pdense.
Returns a hermitian copy of a Matrix. The inplace version only works on a square operand.
Simplistic random number functions.
Stores a lower triangular matrix in rectangular packed format.
Definition: conjugate.h:22


Results: [PASS]

std::complex<float>
|hemm('L')-phemm('L')| = 2.26573e-07
|hemm('R')-phemm('R')| = 2.72961e-07
std::complex<double>
|hemm('L')-phemm('L')| = 4.0634e-16
|hemm('R')-phemm('R')| = 4.94636e-16


Go back to Summary of /test programs.