MyraMath
zher2kL


Source: tests/dense/her2k1.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>
14 
15 // Algorithms.
17 #include <myramath/dense/gemm.h>
18 #include <myramath/dense/her2k.h>
23 #include <myramath/dense/tril.h>
24 
25 // Reporting.
26 #include <tests/myratest.h>
27 
28 using namespace myra;
29 
30 namespace {
31 
32 template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
33  {
34  // Useful typedef.
35  typedef typename ReflectPrecision<Number>::type Precision;
36  myra::out() << typestring<Number>() << std::endl;
37  // Make random nonsquare matrices 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  Number one(1);
43  // Compare her2k_inplace('N') to gemm_inplace(A,'N',B,'H'); gemm_inplace(B,'N',A,'H');
44  {
45  auto C1 = Matrix<Number>::random(I,I);
46  C1 = C1 + hermitian(C1);
47  Matrix<Number> C2 = C1;
48  gemm_inplace(C1,A,'N',B,'H',alpha,beta);
49  gemm_inplace(C1,B,'N',A,'H',conjugate(alpha),one);
50  her2k_inplace(C2,'L',A,B,'N',alpha,beta);
51  Precision error = frobenius(tril(C1-C2));
52  myra::out() << " |her2k('N')-gemm('N','H')| = " << error << std::endl;
53  REQUIRE(error < tolerance);
54  }
55  // Compare her2k_inplace('T') to gemm_inplace(A,'T',B,'C'); gemm_inplace(B,'T',A,'C');
56  {
57  auto C1 = Matrix<Number>::random(J,J);
58  C1 = C1 + hermitian(C1);
59  Matrix<Number> C2 = C1;
60  gemm_inplace(C1,A,'T',B,'C',alpha,beta);
61  gemm_inplace(C1,B,'T',A,'C',conjugate(alpha),one);
62  her2k_inplace(C2,'L',A,B,'T',alpha,beta);
63  Precision error = frobenius(tril(C1-C2));
64  myra::out() << " |her2k('T')-gemm('T','C')| = " << error << std::endl;
65  REQUIRE(error < tolerance);
66  }
67  // Compare her2k_inplace('C') to gemm_inplace(A,'C',B,'T'); gemm_inplace(B,'C',A,'T');
68  {
69  auto C1 = Matrix<Number>::random(I,I);
70  C1 = C1 + hermitian(C1);
71  Matrix<Number> C2 = C1;
72  gemm_inplace(C1,A,'C',B,'T',alpha,beta);
73  gemm_inplace(C1,B,'C',A,'T',conjugate(alpha),one);
74  her2k_inplace(C2,'L',A,B,'C',alpha,beta);
75  Precision error = frobenius(tril(C1-C2));
76  myra::out() << " |her2k('C')-gemm('C','H')| = " << error << std::endl;
77  REQUIRE(error < tolerance);
78  }
79  // Compare her2k_inplace('H') to gemm_inplace(A,'H',B,'N'); gemm_inplace(B,'H',A,'N');
80  {
81  auto C1 = Matrix<Number>::random(J,J);
82  C1 = C1 + hermitian(C1);
83  Matrix<Number> C2 = C1;
84  gemm_inplace(C1,A,'H',B,'N',alpha,beta);
85  gemm_inplace(C1,B,'H',A,'N',conjugate(alpha),one);
86  her2k_inplace(C2,'L',A,B,'H',alpha,beta);
87  Precision error = frobenius(tril(C1-C2));
88  myra::out() << " |her2k('H')-gemm('H','N')| = " << error << std::endl;
89  REQUIRE(error < tolerance);
90  }
91  }
92 
93 } // namespace
94 
95 ADD_TEST("sher2kL","[dense][blas]")
96  { test<NumberS>(57,24,1.0e-4f); }
97 
98 ADD_TEST("dher2kL","[dense][blas]")
99  { test<NumberD>(57,24,1.0e-8); }
100 
101 ADD_TEST("cher2kL","[dense][blas]")
102  { test<NumberC>(57,24,1.0e-4f); }
103 
104 ADD_TEST("zher2kL","[dense][blas]")
105  { test<NumberZ>(57,24,1.0e-8); }
106 
Returns a conjugated copy of a Matrix or Vector. Or, conjugate one inplace.
Tabulates an IxJ matrix. Allows random access, has column major layout to be compatible with BLAS/LAP...
Definition: bdsqr.h:20
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
Routines for hermitian rank-2k updates, a specialized form of Matrix*Matrix multiplication.
Definition: syntax.dox:1
Returns a transposed copy of a Matrix. The inplace version only works on a square operand...
Returns the lower triangle of a dense Matrix.
Various utility functions/classes related to scalar Number types.
General purpose dense matrix container, O(i*j) storage.
Reflects Precision trait for a Number, scalar Number types should specialize it.
Definition: Number.h:33
Returns a hermitian copy of a Matrix. The inplace version only works on a square operand.
Simplistic random number functions.
Variety of routines all for dense Matrix*Matrix multiplies. Delegates to the BLAS.


Results: [PASS]

std::complex<double>
|her2k('N')-gemm('N','H')| = 0
|her2k('T')-gemm('T','C')| = 0
|her2k('C')-gemm('C','H')| = 0
|her2k('H')-gemm('H','N')| = 0


Go back to Summary of /test programs.