MyraMath
pherk1


Source: tests/pdense/pherk1.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>
15 
16 // Serial algorithms.
18 #include <myramath/dense/herk.h>
21 
22 // Parallel algorithms.
23 #include <myramath/pdense/pherk.h>
25 
26 // Reporting.
27 #include <tests/myratest.h>
28 
29 using namespace myra;
30 typedef pdense::Options Options;
31 
32 namespace {
33 
34 // Tests pherk_inplace() against herk_inplace()
35 template<class Precision> void test(int I, int J, Precision tolerance)
36  {
37  typedef std::complex<Precision> Number;
38  myra::out() << typestring<Number>() << std::endl;
39  // Generate random matrix A.
40  auto A = Matrix<Number>::random(I,J);
41  Precision alpha = random<Precision>();
42  Precision beta = random<Precision>();
43  // Initialize options.
44  auto options = Options::create().set_nthreads(4).set_blocksize(128);
45  // Check herk('U','N')
46  {
47  auto C_serial = Matrix<Number>::random(I,I);
48  C_serial = C_serial + hermitian(C_serial);
49  auto C_parallel = C_serial;
50  herk_inplace(C_serial,'U',A,'N',alpha,beta);
51  pherk_inplace(C_parallel,'U',A,'N',alpha,beta,options);
52  Precision error = frobenius(C_serial-C_parallel) / frobenius(C_serial);
53  myra::out() << " |herk('U','N')-pherk('U','N')| = " << error << std::endl;
54  REQUIRE(error < tolerance);
55  }
56  // Check herk('U','T')
57  {
58  auto C_serial = Matrix<Number>::random(J,J);
59  C_serial = C_serial + hermitian(C_serial);
60  auto C_parallel = C_serial;
61  herk_inplace(C_serial,'U',A,'T',alpha,beta);
62  pherk_inplace(C_parallel,'U',A,'T',alpha,beta,options);
63  Precision error = frobenius(C_serial-C_parallel) / frobenius(C_serial);
64  myra::out() << " |herk('U','T')-pherk('U','T')| = " << error << std::endl;
65  REQUIRE(error < tolerance);
66  }
67  // Check herk('U','H')
68  {
69  auto C_serial = Matrix<Number>::random(J,J);
70  C_serial = C_serial + hermitian(C_serial);
71  auto C_parallel = C_serial;
72  herk_inplace(C_serial,'U',A,'H',alpha,beta);
73  pherk_inplace(C_parallel,'U',A,'H',alpha,beta,options);
74  Precision error = frobenius(C_serial-C_parallel) / frobenius(C_serial);
75  myra::out() << " |herk('U','H')-pherk('U','H')| = " << error << std::endl;
76  REQUIRE(error < tolerance);
77  }
78  // Check herk('U','C')
79  {
80  auto C_serial = Matrix<Number>::random(I,I);
81  C_serial = C_serial + hermitian(C_serial);
82  auto C_parallel = C_serial;
83  herk_inplace(C_serial,'U',A,'C',alpha,beta);
84  pherk_inplace(C_parallel,'U',A,'C',alpha,beta,options);
85  Precision error = frobenius(C_serial-C_parallel) / frobenius(C_serial);
86  myra::out() << " |herk('U','C')-pherk('U','C')| = " << error << std::endl;
87  REQUIRE(error < tolerance);
88  }
89  // Check herk('L','N')
90  {
91  auto C_serial = Matrix<Number>::random(I,I);
92  C_serial = C_serial + hermitian(C_serial);
93  auto C_parallel = C_serial;
94  herk_inplace(C_serial,'L',A,'N',alpha,beta);
95  pherk_inplace(C_parallel,'L',A,'N',alpha,beta);
96  Precision error = frobenius(C_serial-C_parallel) / frobenius(C_serial);
97  myra::out() << " |herk('L','N')-pherk('L','N')| = " << error << std::endl;
98  REQUIRE(error < tolerance);
99  }
100  // Check herk('L','T')
101  {
102  auto C_serial = Matrix<Number>::random(J,J);
103  C_serial = C_serial + hermitian(C_serial);
104  auto C_parallel = C_serial;
105  herk_inplace(C_serial,'L',A,'T',alpha,beta);
106  pherk_inplace(C_parallel,'L',A,'T',alpha,beta);
107  Precision error = frobenius(C_serial-C_parallel) / frobenius(C_serial);
108  myra::out() << " |herk('L','T')-pherk('L','T')| = " << error << std::endl;
109  REQUIRE(error < tolerance);
110  }
111  // Check herk('L','H')
112  {
113  auto C_serial = Matrix<Number>::random(J,J);
114  C_serial = C_serial + hermitian(C_serial);
115  auto C_parallel = C_serial;
116  herk_inplace(C_serial,'L',A,'H',alpha,beta);
117  pherk_inplace(C_parallel,'L',A,'H',alpha,beta);
118  Precision error = frobenius(C_serial-C_parallel) / frobenius(C_serial);
119  myra::out() << " |herk('L','H')-pherk('L','H')| = " << error << std::endl;
120  REQUIRE(error < tolerance);
121  }
122  // Check herk('L','C')
123  {
124  auto C_serial = Matrix<Number>::random(I,I);
125  C_serial = C_serial + hermitian(C_serial);
126  auto C_parallel = C_serial;
127  herk_inplace(C_serial,'L',A,'C',alpha,beta);
128  pherk_inplace(C_parallel,'L',A,'C',alpha,beta);
129  Precision error = frobenius(C_serial-C_parallel) / frobenius(C_serial);
130  myra::out() << " |herk('L','C')-pherk('L','C')| = " << error << std::endl;
131  REQUIRE(error < tolerance);
132  }
133  }
134 
135 } // namespace
136 
137 ADD_TEST("pherk1","[pdense][parallel]")
138  {
139  int I = 512;
140  int J = 256;
141  test<float > (I,J,1.0e-4f);
142  test<double> (I,J,1.0e-8);
143  }
Interface class for representing subranges of dense Matrix&#39;s.
Routines for hermitian rank-k updates, a specialized form of Matrix*Matrix multiplication.
Routines for computing Frobenius norms of various algebraic containers.
Thread parallel version of dense/herk.h, hermitian rank-k updates.
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
Definition: syntax.dox:1
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.


Results: [PASS]

std::complex<float>
|herk('U','N')-pherk('U','N')| = 0
|herk('U','T')-pherk('U','T')| = 0
|herk('U','H')-pherk('U','H')| = 0
|herk('U','C')-pherk('U','C')| = 0
|herk('L','N')-pherk('L','N')| = 0
|herk('L','T')-pherk('L','T')| = 0
|herk('L','H')-pherk('L','H')| = 0
|herk('L','C')-pherk('L','C')| = 0
std::complex<double>
|herk('U','N')-pherk('U','N')| = 0
|herk('U','T')-pherk('U','T')| = 0
|herk('U','H')-pherk('U','H')| = 0
|herk('U','C')-pherk('U','C')| = 0
|herk('L','N')-pherk('L','N')| = 0
|herk('L','T')-pherk('L','T')| = 0
|herk('L','H')-pherk('L','H')| = 0
|herk('L','C')-pherk('L','C')| = 0


Go back to Summary of /test programs.