MyraMath
pher2k1


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


Results: [PASS]

std::complex<float>
|her2k('U','N')-pher2k('U','N')| = 0
|her2k('U','T')-pher2k('U','T')| = 0
|her2k('U','H')-pher2k('U','H')| = 0
|her2k('U','C')-pher2k('U','C')| = 0
|her2k('L','N')-pher2k('L','N')| = 3.78802e-08
|her2k('L','T')-pher2k('L','T')| = 0
|her2k('L','H')-pher2k('L','H')| = 0
|her2k('L','C')-pher2k('L','C')| = 3.76193e-08
std::complex<double>
|her2k('U','N')-pher2k('U','N')| = 1.11281e-16
|her2k('U','T')-pher2k('U','T')| = 1.16702e-16
|her2k('U','H')-pher2k('U','H')| = 1.17732e-16
|her2k('U','C')-pher2k('U','C')| = 1.11603e-16
|her2k('L','N')-pher2k('L','N')| = 1.08926e-16
|her2k('L','T')-pher2k('L','T')| = 0
|her2k('L','H')-pher2k('L','H')| = 0
|her2k('L','C')-pher2k('L','C')| = 1.08844e-16


Go back to Summary of /test programs.