1 //
2 // Copyright Markus Rickert 2008
3 //
4 // Distributed under the Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7 //
8
9 #include <algorithm>
10 #include <boost/numeric/ublas/io.hpp>
11 #include <boost/numeric/bindings/ublas/matrix.hpp>
12 #include <boost/numeric/bindings/ublas/matrix_proxy.hpp>
13 #include <boost/numeric/bindings/ublas/vector.hpp>
14 #include <boost/numeric/bindings/trans.hpp>
15 #include <boost/numeric/bindings/blas/level3/gemm.hpp>
16 #include <boost/numeric/bindings/blas/level1/axpy.hpp>
17
18 namespace bindings = boost::numeric::bindings;
19
20 int
main(int argc,char ** argv)21 main(int argc, char** argv)
22 {
23 // a * b' = C ; a' * b = d
24 {
25 boost::numeric::ublas::vector<double> a(3);
26 for(std::size_t i = 0; i < a.size(); ++i) a(i) = i;
27 std::cout << "a=" << a << std::endl;
28
29 boost::numeric::ublas::vector<double> b(3);
30 for(std::size_t i = 0; i < b.size(); ++i) b(i) = i;
31 std::cout << "b=" << b << std::endl;
32
33 boost::numeric::ublas::matrix<double, boost::numeric::ublas::column_major> c(3, 3);
34 boost::numeric::bindings::blas::gemm(
35 1.0, a, bindings::trans(b), 0.0, c
36 );
37 std::cout << "C=" << c << std::endl;
38
39 boost::numeric::ublas::vector<double> d(1);
40 boost::numeric::bindings::blas::gemm(
41 1.0, bindings::trans(a), b, 0.0, d
42 );
43 std::cout << "d=" << d << std::endl;
44 }
45
46 std::cout << std::endl;
47
48 // a * b' = C ; a' * b = d
49 {
50 boost::numeric::ublas::bounded_vector<double, 3> a;
51 for(std::size_t i = 0; i < a.size(); ++i) a(i) = i;
52 std::cout << "a=" << a << std::endl;
53
54 boost::numeric::ublas::bounded_vector<double, 3> b;
55 for(std::size_t i = 0; i < b.size(); ++i) b(i) = i;
56 std::cout << "b=" << b << std::endl;
57
58 boost::numeric::ublas::bounded_matrix<double, 3, 3, boost::numeric::ublas::column_major> c;
59 boost::numeric::bindings::blas::gemm(
60 1.0, a, bindings::trans(b), 0.0, c
61 );
62 std::cout << "C=" << c << std::endl;
63
64 boost::numeric::ublas::bounded_vector<double, 1> d;
65 boost::numeric::bindings::blas::gemm(
66 1.0, bindings::trans(a), b, 0.0, d
67 );
68 std::cout << "d=" << d << std::endl;
69 }
70
71 std::cout << std::endl;
72
73 // A * B = C
74 {
75 boost::numeric::ublas::bounded_matrix<double, 4, 3, boost::numeric::ublas::column_major> a;
76 for(std::size_t i = 0; i < a.size1(); ++i) for(std::size_t j = 0; j < a.size2(); ++j) a(i, j) = i * a.size2() + j;
77 std::cout << "A=" << a << std::endl;
78
79 boost::numeric::ublas::bounded_matrix<double, 3, 4, boost::numeric::ublas::column_major> b;
80 for(std::size_t i = 0; i < b.size1(); ++i) for(std::size_t j = 0; j < b.size2(); ++j) b(i, j) = i * b.size2() + j;
81 std::cout << "B=" << b << std::endl;
82
83 boost::numeric::ublas::bounded_matrix<double, 4, 4, boost::numeric::ublas::column_major> c;
84 //boost::numeric::bindings::blas::gemm(a, b, c);
85 boost::numeric::bindings::blas::gemm(
86 1.0, a, b, 0.0, c);
87 std::cout << "C=" << c << std::endl;
88 }
89
90 std::cout << std::endl;
91
92 // A[0:3;0:2] * B[0:2;0:3] = C
93 {
94 boost::numeric::ublas::bounded_matrix<double, 4, 3, boost::numeric::ublas::column_major> a;
95 for(std::size_t i = 0; i < a.size1(); ++i) for(std::size_t j = 0; j < a.size2(); ++j) a(i, j) = i * a.size2() + j;
96 std::cout << "A=" << a << std::endl;
97
98 boost::numeric::ublas::matrix_range<
99 boost::numeric::ublas::bounded_matrix<double, 4, 3, boost::numeric::ublas::column_major>
100 > a2 = boost::numeric::ublas::subrange(a, 0, 3, 0, 2);
101 std::cout << "A2=" << a2 << std::endl;
102
103 boost::numeric::ublas::bounded_matrix<double, 3, 4, boost::numeric::ublas::column_major> b;
104 for(std::size_t i = 0; i < b.size1(); ++i) for(std::size_t j = 0; j < b.size2(); ++j) b(i, j) = i * b.size2() + j;
105 std::cout << "B=" << b << std::endl;
106
107 boost::numeric::ublas::matrix_range<
108 boost::numeric::ublas::bounded_matrix<double, 3, 4, boost::numeric::ublas::column_major>
109 > b2 = boost::numeric::ublas::subrange(b, 0, 2, 0, 3);
110 std::cout << "B2=" << b2 << std::endl;
111
112 boost::numeric::ublas::bounded_matrix<double, 4, 4, boost::numeric::ublas::column_major> c;
113 std::fill(c.data().begin(), c.data().end(), 0.0);
114 boost::numeric::ublas::matrix_range<
115 boost::numeric::ublas::bounded_matrix<double, 4, 4, boost::numeric::ublas::column_major>
116 > c2 = boost::numeric::ublas::subrange(c, 0, 3, 0, 3);
117 //boost::numeric::bindings::blas::gemm(a2, b2, c2);
118 boost::numeric::bindings::blas::gemm(
119 1.0, a2, b2, 0.0, c2);
120 std::cout << "C2=" << c2 << std::endl;
121 std::cout << "C=" << c << std::endl;
122 }
123
124 std::cout << std::endl;
125
126 // a + b = b ; b - a = b
127 {
128 boost::numeric::ublas::bounded_vector<double, 3> a;
129 for(std::size_t i = 0; i < a.size(); ++i) a(i) = i;
130 std::cout << "a=" << a << std::endl;
131
132 boost::numeric::ublas::bounded_vector<double, 3> b;
133 for(std::size_t i = 0; i < b.size(); ++i) b(i) = i;
134 std::cout << "b=" << b << std::endl;
135
136 boost::numeric::bindings::blas::axpy(1.0, a, b);
137 std::cout << "b=" << b << std::endl;
138
139 boost::numeric::bindings::blas::axpy(-1.0, a, b);
140 std::cout << "b=" << b << std::endl;
141 }
142
143 std::cout << std::endl;
144
145 // b + c = c ; c - b = c
146 {
147 boost::numeric::ublas::matrix<double, boost::numeric::ublas::column_major> a(5, 5);
148 for(std::size_t i = 0; i < a.size1(); ++i) for(std::size_t j = 0; j < a.size2(); ++j) a(i, j) = i * a.size2() + j;
149 std::cout << "A=" << a << std::endl;
150
151 boost::numeric::ublas::matrix_vector_range<
152 boost::numeric::ublas::matrix<double, boost::numeric::ublas::column_major>
153 > b(a, boost::numeric::ublas::range(1, 4), boost::numeric::ublas::range(0, 3));
154 std::cout << "b=" << b << std::endl;
155
156 boost::numeric::ublas::matrix_vector_slice<
157 boost::numeric::ublas::matrix<double, boost::numeric::ublas::column_major>
158 > c(a, boost::numeric::ublas::slice(0, 1, 3), boost::numeric::ublas::slice(3, 0, 3));
159 std::cout << "c=" << c << std::endl;
160
161 boost::numeric::bindings::blas::axpy(1.0, b, c);
162 std::cout << "c=" << c << std::endl;
163
164 boost::numeric::bindings::blas::axpy(-1.0, b, c);
165 std::cout << "c=" << c << std::endl;
166 }
167
168 std::cout << std::endl;
169
170 // b + c = c ; c - b = c
171 {
172 boost::numeric::ublas::bounded_matrix<double, 5, 5, boost::numeric::ublas::column_major> a;
173 for(std::size_t i = 0; i < a.size1(); ++i) for(std::size_t j = 0; j < a.size2(); ++j) a(i, j) = i * a.size2() + j;
174 std::cout << "A=" << a << std::endl;
175
176 boost::numeric::ublas::matrix_vector_range<
177 boost::numeric::ublas::bounded_matrix<double, 5, 5, boost::numeric::ublas::column_major>
178 > b(a, boost::numeric::ublas::range(1, 4), boost::numeric::ublas::range(0, 3));
179 std::cout << "b=" << b << std::endl;
180
181 boost::numeric::ublas::matrix_vector_slice<
182 boost::numeric::ublas::bounded_matrix<double, 5, 5, boost::numeric::ublas::column_major>
183 > c(a, boost::numeric::ublas::slice(0, 1, 3), boost::numeric::ublas::slice(3, 0, 3));
184 std::cout << "c=" << c << std::endl;
185
186 boost::numeric::bindings::blas::axpy(1.0, b, c);
187 std::cout << "c=" << c << std::endl;
188
189 boost::numeric::bindings::blas::axpy(-1.0, b, c);
190 std::cout << "c=" << c << std::endl;
191 }
192
193 return 0;
194 }
195