1
2 // BLAS level 3
3 // hermitian matrices, herk
4
5 #include <stddef.h>
6 #include <iostream>
7 #include <boost/numeric/bindings/blas/level3.hpp>
8 #include <boost/numeric/bindings/ublas/hermitian.hpp>
9 #include <boost/numeric/bindings/upper.hpp>
10 #include <boost/numeric/bindings/lower.hpp>
11 #include <boost/numeric/bindings/conj.hpp>
12 #include "utils.h"
13
14 namespace ublas = boost::numeric::ublas;
15 namespace blas = boost::numeric::bindings::blas;
16 namespace bindings = boost::numeric::bindings;
17
18 using std::cout;
19 using std::cin;
20 using std::endl;
21
22 typedef double real_t;
23 typedef std::complex<real_t> cmplx_t;
24
25 typedef ublas::matrix<cmplx_t, ublas::column_major> cm_t;
26 typedef ublas::matrix<cmplx_t, ublas::row_major> rm_t;
27 typedef ublas::hermitian_adaptor<cm_t, ublas::upper> ucha_t;
28 typedef ublas::hermitian_adaptor<cm_t, ublas::lower> lcha_t;
29 typedef ublas::hermitian_adaptor<rm_t, ublas::upper> urha_t;
30 typedef ublas::hermitian_adaptor<rm_t, ublas::lower> lrha_t;
31
32 #define N 3
33 #define K 4
34
main()35 int main() {
36
37 cm_t ac (N, K);
38 rm_t ar (N, K);
39
40 ac(0,0) = ar(0,0) = cmplx_t (1., 1.);
41 ac(1,0) = ar(1,0) = cmplx_t (2., 1.);
42 ac(2,0) = ar(2,0) = cmplx_t (3., 1.);
43 #if (K > 1)
44 ac(0,1) = ar(0,1) = cmplx_t (1., 1.);
45 ac(1,1) = ar(1,1) = cmplx_t (2., 1.);
46 ac(2,1) = ar(2,1) = cmplx_t (3., 1.);
47 #if (K > 2)
48 ac(0,2) = ar(0,2) = cmplx_t (1., 1.);
49 ac(1,2) = ar(1,2) = cmplx_t (2., 1.);
50 ac(2,2) = ar(2,2) = cmplx_t (3., 1.);
51 #if (K > 3)
52 ac(0,3) = ar(0,3) = cmplx_t (1., 1.);
53 ac(1,3) = ar(1,3) = cmplx_t (2., 1.);
54 ac(2,3) = ar(2,3) = cmplx_t (3., 1.);
55 #endif
56 #endif
57 #endif
58
59 print_m (ac, "ac");
60 cout << endl;
61 print_m (ar, "ar");
62 cout << endl << endl;
63
64 cm_t cmu (N, N);
65 cm_t cml (N, N);
66 rm_t rmu (N, N);
67 rm_t rml (N, N);
68 ucha_t ucha (cmu);
69 lcha_t lcha (cml);
70 urha_t urha (rmu);
71 lrha_t lrha (rml);
72
73 blas::herk (1.0, ac, 0.0, ucha);
74 blas::herk (1.0, ac, 0.0, lcha);
75 blas::herk (1.0, ar, 0.0, urha);
76 blas::herk (1.0, ar, 0.0, lrha);
77
78 print_m (ucha, "ucha");
79 cout << endl;
80 print_m (lcha, "lcha");
81 cout << endl;
82 print_m (urha, "urha");
83 cout << endl;
84 print_m (lrha, "lrha");
85 cout << endl << endl;
86
87 // part 2
88
89 cm_t act (ublas::herm (ac));
90 rm_t art (ublas::herm (ar));
91 print_m (act, "act");
92 cout << endl;
93 print_m (art, "art");
94 cout << endl << endl;
95
96 init_m (cmu, const_val<cmplx_t> (cmplx_t (0, 0)));
97 init_m (cml, const_val<cmplx_t> (cmplx_t (0, 0)));
98 init_m (rmu, const_val<cmplx_t> (cmplx_t (0, 0)));
99 init_m (rml, const_val<cmplx_t> (cmplx_t (0, 0)));
100
101 blas::herk ( 1.0, bindings::conj(act), 0.0, bindings::upper(cmu));
102 blas::herk ( 1.0, bindings::conj(act), 0.0, bindings::lower(cml));
103 blas::herk ( 1.0, bindings::conj(art), 0.0, bindings::upper(rmu));
104 blas::herk ( 1.0, bindings::conj(art), 0.0, bindings::lower(rml));
105
106 print_m (cmu, "cmu");
107 cout << endl;
108 print_m (cml, "cml");
109 cout << endl;
110 print_m (rmu, "rmu");
111 cout << endl;
112 print_m (rml, "rml");
113 cout << endl;
114 }
115
116