1 /* =========================================================================
2 Copyright (c) 2010-2016, Institute for Microelectronics,
3 Institute for Analysis and Scientific Computing,
4 TU Wien.
5 Portions of this software are copyright by UChicago Argonne, LLC.
6
7 -----------------
8 ViennaCL - The Vienna Computing Library
9 -----------------
10
11 Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12
13 (A list of authors and contributors can be found in the PDF manual)
14
15 License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17
18
19
20 /** \file tests/src/scheduler_matrix_vector.cpp Tests the scheduler for matrix-vector-operations.
21 * \test Tests the scheduler for matrix-vector-operations.
22 **/
23
24 //
25 // *** System
26 //
27 #include <iostream>
28 #include <vector>
29
30 //
31 // *** ViennaCL
32 //
33 #include "viennacl/scalar.hpp"
34 #include "viennacl/matrix.hpp"
35 #include "viennacl/vector.hpp"
36 #include "viennacl/linalg/prod.hpp"
37 #include "viennacl/linalg/norm_2.hpp"
38 #include "viennacl/linalg/direct_solve.hpp"
39 #include "viennacl/linalg/lu.hpp"
40 #include "viennacl/tools/random.hpp"
41
42 #include "viennacl/scheduler/execute.hpp"
43 #include "viennacl/scheduler/io.hpp"
44
45
46 //
47 // -------------------------------------------------------------
48 //
49 template<typename ScalarType>
diff(ScalarType & s1,viennacl::scalar<ScalarType> & s2)50 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
51 {
52 viennacl::backend::finish();
53 if (s1 != s2)
54 return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
55 return 0;
56 }
57
58 template<typename ScalarType, typename VCLVectorType>
59 ScalarType diff(std::vector<ScalarType> const & v1, VCLVectorType const & v2)
60 {
61 std::vector<ScalarType> v2_cpu(v2.size());
62 viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
63 viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
64
65 ScalarType norm_inf = 0;
66 for (unsigned int i=0;i<v1.size(); ++i)
67 {
68 if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
69 {
70 ScalarType tmp = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
71 if (tmp > norm_inf)
72 norm_inf = tmp;
73 }
74 }
75
76 return norm_inf;
77 }
78
79 template<typename ScalarType, typename VCLMatrixType>
80 ScalarType diff(std::vector<std::vector<ScalarType> > const & mat1, VCLMatrixType const & mat2)
81 {
82 std::vector<std::vector<ScalarType> > mat2_cpu(mat2.size1(), std::vector<ScalarType>(mat2.size2()));
83 viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
84 viennacl::copy(mat2, mat2_cpu);
85 ScalarType ret = 0;
86 ScalarType act = 0;
87
88 for (std::size_t i = 0; i < mat2_cpu.size(); ++i)
89 {
90 for (std::size_t j = 0; j < mat2_cpu[i].size(); ++j)
91 {
92 act = std::fabs(mat2_cpu[i][j] - mat1[i][j]) / std::max( std::fabs(mat2_cpu[i][j]), std::fabs(mat1[i][j]) );
93 if (act > ret)
94 ret = act;
95 }
96 }
97 //std::cout << ret << std::endl;
98 return ret;
99 }
100 //
101 // -------------------------------------------------------------
102 //
103
104 template<typename NumericT, typename Epsilon,
105 typename STLMatrixType, typename STLVectorType,
106 typename VCLMatrixType, typename VCLVectorType1, typename VCLVectorType2>
test_prod_rank1(Epsilon const & epsilon,STLMatrixType & std_m1,STLVectorType & std_v1,STLVectorType & std_v2,VCLMatrixType & vcl_m1,VCLVectorType1 & vcl_v1,VCLVectorType2 & vcl_v2)107 int test_prod_rank1(Epsilon const & epsilon,
108 STLMatrixType & std_m1, STLVectorType & std_v1, STLVectorType & std_v2,
109 VCLMatrixType & vcl_m1, VCLVectorType1 & vcl_v1, VCLVectorType2 & vcl_v2)
110 {
111 int retval = EXIT_SUCCESS;
112
113 // sync data:
114 viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
115 viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
116 viennacl::copy(std_m1, vcl_m1);
117
118 /* TODO: Add rank-1 operations here */
119
120 //reset vcl_matrix:
121 viennacl::copy(std_m1, vcl_m1);
122
123 // --------------------------------------------------------------------------
124 std::cout << "Matrix-Vector product" << std::endl;
125 for (std::size_t i=0; i<std_m1.size(); ++i)
126 {
127 std_v1[i] = 0;
128 for (std::size_t j=0; j<std_m1[i].size(); ++j)
129 std_v1[i] += std_m1[i][j] * std_v2[j];
130 }
131 {
132 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::prod(vcl_m1, vcl_v2));
133 viennacl::scheduler::execute(my_statement);
134 }
135
136 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
137 {
138 std::cout << "# Error at operation: matrix-vector product" << std::endl;
139 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
140 retval = EXIT_FAILURE;
141 }
142
143 std::cout << "Matrix-Vector product with inplace-add" << std::endl;
144 for (std::size_t i=0; i<std_m1.size(); ++i)
145 {
146 for (std::size_t j=0; j<std_m1[i].size(); ++j)
147 std_v1[i] += std_m1[i][j] * std_v2[j];
148 }
149 {
150 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_add(), viennacl::linalg::prod(vcl_m1, vcl_v2));
151 viennacl::scheduler::execute(my_statement);
152 }
153
154 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
155 {
156 std::cout << "# Error at operation: matrix-vector product" << std::endl;
157 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
158 retval = EXIT_FAILURE;
159 }
160
161 std::cout << "Matrix-Vector product with inplace-sub" << std::endl;
162 for (std::size_t i=0; i<std_m1.size(); ++i)
163 {
164 for (std::size_t j=0; j<std_m1[i].size(); ++j)
165 std_v1[i] -= std_m1[i][j] * std_v2[j];
166 }
167 {
168 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_sub(), viennacl::linalg::prod(vcl_m1, vcl_v2));
169 viennacl::scheduler::execute(my_statement);
170 }
171
172 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
173 {
174 std::cout << "# Error at operation: matrix-vector product" << std::endl;
175 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
176 retval = EXIT_FAILURE;
177 }
178
179 // --------------------------------------------------------------------------
180 /*
181 std::cout << "Matrix-Vector product with scaled matrix" << std::endl;
182 std_v1 = viennacl::linalg::prod(NumericT(2.0) * std_m1, std_v2);
183 {
184 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::prod(NumericT(2.0) * vcl_m1, vcl_v2));
185 viennacl::scheduler::execute(my_statement);
186 }
187
188 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
189 {
190 std::cout << "# Error at operation: matrix-vector product" << std::endl;
191 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
192 retval = EXIT_FAILURE;
193 }*/
194
195 // --------------------------------------------------------------------------
196 std::cout << "Matrix-Vector product with scaled vector" << std::endl;
197 /*
198 std_v1 = viennacl::linalg::prod(std_m1, NumericT(2.0) * std_v2);
199 {
200 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::prod(vcl_m1, NumericT(2.0) * vcl_v2));
201 viennacl::scheduler::execute(my_statement);
202 }
203
204 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
205 {
206 std::cout << "# Error at operation: matrix-vector product" << std::endl;
207 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
208 retval = EXIT_FAILURE;
209 }*/
210
211 // --------------------------------------------------------------------------
212 std::cout << "Matrix-Vector product with scaled matrix and scaled vector" << std::endl;
213 /*
214 std_v1 = viennacl::linalg::prod(NumericT(2.0) * std_m1, NumericT(2.0) * std_v2);
215 {
216 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::prod(NumericT(2.0) * vcl_m1, NumericT(2.0) * vcl_v2));
217 viennacl::scheduler::execute(my_statement);
218 }
219
220 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
221 {
222 std::cout << "# Error at operation: matrix-vector product" << std::endl;
223 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
224 retval = EXIT_FAILURE;
225 }*/
226
227
228 // --------------------------------------------------------------------------
229 std::cout << "Matrix-Vector product with scaled add" << std::endl;
230 NumericT alpha = static_cast<NumericT>(2.786);
231 NumericT beta = static_cast<NumericT>(3.1415);
232 viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
233 viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
234
235 for (std::size_t i=0; i<std_m1.size(); ++i)
236 {
237 std_v1[i] = 0;
238 for (std::size_t j=0; j<std_m1[i].size(); ++j)
239 std_v1[i] += std_m1[i][j] * std_v2[j];
240 std_v1[i] = alpha * std_v1[i] - beta * std_v1[i];
241 }
242 {
243 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) - beta * vcl_v1);
244 viennacl::scheduler::execute(my_statement);
245 }
246
247 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
248 {
249 std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
250 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
251 retval = EXIT_FAILURE;
252 }
253
254 std::cout << "Matrix-Vector product with scaled add, inplace-add" << std::endl;
255 viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
256 viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
257
258 for (std::size_t i=0; i<std_m1.size(); ++i)
259 {
260 NumericT tmp = 0;
261 for (std::size_t j=0; j<std_m1[i].size(); ++j)
262 tmp += std_m1[i][j] * std_v2[j];
263 std_v1[i] += alpha * tmp - beta * std_v1[i];
264 }
265 {
266 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_add(), alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) - beta * vcl_v1);
267 viennacl::scheduler::execute(my_statement);
268 }
269
270 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
271 {
272 std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
273 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
274 retval = EXIT_FAILURE;
275 }
276
277 std::cout << "Matrix-Vector product with scaled add, inplace-sub" << std::endl;
278 viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
279 viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
280
281 for (std::size_t i=0; i<std_m1.size(); ++i)
282 {
283 NumericT tmp = 0;
284 for (std::size_t j=0; j<std_m1[i].size(); ++j)
285 tmp += std_m1[i][j] * std_v2[j];
286 std_v1[i] -= alpha * tmp - beta * std_v1[i];
287 }
288 {
289 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_sub(), alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) - beta * vcl_v1);
290 viennacl::scheduler::execute(my_statement);
291 }
292
293 if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
294 {
295 std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
296 std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
297 retval = EXIT_FAILURE;
298 }
299
300 // --------------------------------------------------------------------------
301
302 viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
303 viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
304
305 std::cout << "Transposed Matrix-Vector product" << std::endl;
306 for (std::size_t i=0; i<std_m1[0].size(); ++i)
307 {
308 std_v2[i] = 0;
309 for (std::size_t j=0; j<std_m1.size(); ++j)
310 std_v2[i] += std_m1[j][i] * std_v1[j];
311 }
312 {
313 viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_assign(), viennacl::linalg::prod(trans(vcl_m1), vcl_v1));
314 viennacl::scheduler::execute(my_statement);
315 }
316
317 if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
318 {
319 std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
320 std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
321 retval = EXIT_FAILURE;
322 }
323
324 std::cout << "Transposed Matrix-Vector product, inplace-add" << std::endl;
325 for (std::size_t i=0; i<std_m1[0].size(); ++i)
326 {
327 for (std::size_t j=0; j<std_m1.size(); ++j)
328 std_v2[i] += std_m1[j][i] * std_v1[j];
329 }
330 {
331 viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_inplace_add(), viennacl::linalg::prod(trans(vcl_m1), vcl_v1));
332 viennacl::scheduler::execute(my_statement);
333 }
334
335 if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
336 {
337 std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
338 std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
339 retval = EXIT_FAILURE;
340 }
341
342 std::cout << "Transposed Matrix-Vector product, inplace-sub" << std::endl;
343 for (std::size_t i=0; i<std_m1[0].size(); ++i)
344 {
345 for (std::size_t j=0; j<std_m1.size(); ++j)
346 std_v2[i] -= std_m1[j][i] * std_v1[j];
347 }
348 {
349 viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_inplace_sub(), viennacl::linalg::prod(trans(vcl_m1), vcl_v1));
350 viennacl::scheduler::execute(my_statement);
351 }
352
353 if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
354 {
355 std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
356 std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
357 retval = EXIT_FAILURE;
358 }
359
360 // --------------------------------------------------------------------------
361 std::cout << "Transposed Matrix-Vector product with scaled add" << std::endl;
362 for (std::size_t i=0; i<std_m1[0].size(); ++i)
363 {
364 NumericT tmp = 0;
365 for (std::size_t j=0; j<std_m1.size(); ++j)
366 tmp += std_m1[j][i] * std_v1[j];
367 std_v2[i] = alpha * tmp + beta * std_v2[i];
368 }
369 {
370 viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_assign(), alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2);
371 viennacl::scheduler::execute(my_statement);
372 }
373
374 if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
375 {
376 std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
377 std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
378 retval = EXIT_FAILURE;
379 }
380
381 std::cout << "Transposed Matrix-Vector product with scaled add, inplace-add" << std::endl;
382 for (std::size_t i=0; i<std_m1[0].size(); ++i)
383 {
384 NumericT tmp = 0;
385 for (std::size_t j=0; j<std_m1.size(); ++j)
386 tmp += std_m1[j][i] * std_v1[j];
387 std_v2[i] += alpha * tmp + beta * std_v2[i];
388 }
389 {
390 viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_inplace_add(), alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2);
391 viennacl::scheduler::execute(my_statement);
392 }
393
394 if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
395 {
396 std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
397 std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
398 retval = EXIT_FAILURE;
399 }
400
401 std::cout << "Transposed Matrix-Vector product with scaled add, inplace-sub" << std::endl;
402 for (std::size_t i=0; i<std_m1[0].size(); ++i)
403 {
404 NumericT tmp = 0;
405 for (std::size_t j=0; j<std_m1.size(); ++j)
406 tmp += std_m1[j][i] * std_v1[j];
407 std_v2[i] -= alpha * tmp + beta * std_v2[i];
408 }
409 {
410 viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_inplace_sub(), alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2);
411 viennacl::scheduler::execute(my_statement);
412 }
413
414 if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
415 {
416 std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
417 std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
418 retval = EXIT_FAILURE;
419 }
420
421 // --------------------------------------------------------------------------
422
423 return retval;
424 }
425
426
427 //
428 // -------------------------------------------------------------
429 //
430 template< typename NumericT, typename F, typename Epsilon >
test(Epsilon const & epsilon)431 int test(Epsilon const& epsilon)
432 {
433 int retval = EXIT_SUCCESS;
434
435 viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
436
437 std::size_t num_rows = 141;
438 std::size_t num_cols = 79;
439
440 // --------------------------------------------------------------------------
441 std::vector<NumericT> std_v1(num_rows);
442 for (std::size_t i = 0; i < std_v1.size(); ++i)
443 std_v1[i] = randomNumber();
444 std::vector<NumericT> std_v2 = std::vector<NumericT>(num_cols, NumericT(3.1415));
445
446
447 std::vector<std::vector<NumericT> > std_m1(std_v1.size(), std::vector<NumericT>(std_v2.size()));
448
449 for (std::size_t i = 0; i < std_m1.size(); ++i)
450 for (std::size_t j = 0; j < std_m1[i].size(); ++j)
451 std_m1[i][j] = static_cast<NumericT>(0.1) * randomNumber();
452
453
454 std::vector<std::vector<NumericT> > std_m2(std_v1.size(), std::vector<NumericT>(std_v1.size()));
455
456 for (std::size_t i = 0; i < std_m2.size(); ++i)
457 {
458 for (std::size_t j = 0; j < std_m2[i].size(); ++j)
459 std_m2[i][j] = static_cast<NumericT>(-0.1) * randomNumber();
460 std_m2[i][i] = static_cast<NumericT>(2) + randomNumber();
461 }
462
463
464 viennacl::vector<NumericT> vcl_v1_native(std_v1.size());
465 viennacl::vector<NumericT> vcl_v1_large(4 * std_v1.size());
466 viennacl::vector_range< viennacl::vector<NumericT> > vcl_v1_range(vcl_v1_large, viennacl::range(3, std_v1.size() + 3));
467 viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v1_slice(vcl_v1_large, viennacl::slice(2, 3, std_v1.size()));
468
469 viennacl::vector<NumericT> vcl_v2_native(std_v2.size());
470 viennacl::vector<NumericT> vcl_v2_large(4 * std_v2.size());
471 viennacl::vector_range< viennacl::vector<NumericT> > vcl_v2_range(vcl_v2_large, viennacl::range(8, std_v2.size() + 8));
472 viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v2_slice(vcl_v2_large, viennacl::slice(6, 2, std_v2.size()));
473
474 viennacl::matrix<NumericT, F> vcl_m1_native(std_m1.size(), std_m1[0].size());
475 viennacl::matrix<NumericT, F> vcl_m1_large(4 * std_m1.size(), 4 * std_m1[0].size());
476 viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m1_range(vcl_m1_large,
477 viennacl::range(8, std_m1.size() + 8),
478 viennacl::range(std_m1[0].size(), 2 * std_m1[0].size()) );
479 viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m1_slice(vcl_m1_large,
480 viennacl::slice(6, 2, std_m1.size()),
481 viennacl::slice(std_m1[0].size(), 2, std_m1[0].size()) );
482
483 viennacl::matrix<NumericT, F> vcl_m2_native(std_m2.size(), std_m2[0].size());
484 viennacl::matrix<NumericT, F> vcl_m2_large(4 * std_m2.size(), 4 * std_m2[0].size());
485 viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m2_range(vcl_m2_large,
486 viennacl::range(8, std_m2.size() + 8),
487 viennacl::range(std_m2[0].size(), 2 * std_m2[0].size()) );
488 viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m2_slice(vcl_m2_large,
489 viennacl::slice(6, 2, std_m2.size()),
490 viennacl::slice(std_m2[0].size(), 2, std_m2[0].size()) );
491
492
493 /* std::cout << "Matrix resizing (to larger)" << std::endl;
494 matrix.resize(2*num_rows, 2*num_cols, true);
495 for (unsigned int i = 0; i < matrix.size1(); ++i)
496 {
497 for (unsigned int j = (i<result.size() ? rhs.size() : 0); j < matrix.size2(); ++j)
498 matrix(i,j) = 0;
499 }
500 vcl_matrix.resize(2*num_rows, 2*num_cols, true);
501 viennacl::copy(vcl_matrix, matrix);
502 if ( std::fabs(diff(matrix, vcl_matrix)) > epsilon )
503 {
504 std::cout << "# Error at operation: matrix resize (to larger)" << std::endl;
505 std::cout << " diff: " << std::fabs(diff(matrix, vcl_matrix)) << std::endl;
506 return EXIT_FAILURE;
507 }
508
509 matrix(12, 14) = NumericT(1.9);
510 matrix(19, 16) = NumericT(1.0);
511 matrix (13, 15) = NumericT(-9);
512 vcl_matrix(12, 14) = NumericT(1.9);
513 vcl_matrix(19, 16) = NumericT(1.0);
514 vcl_matrix (13, 15) = NumericT(-9);
515
516 std::cout << "Matrix resizing (to smaller)" << std::endl;
517 matrix.resize(result.size(), rhs.size(), true);
518 vcl_matrix.resize(result.size(), rhs.size(), true);
519 if ( std::fabs(diff(matrix, vcl_matrix)) > epsilon )
520 {
521 std::cout << "# Error at operation: matrix resize (to smaller)" << std::endl;
522 std::cout << " diff: " << std::fabs(diff(matrix, vcl_matrix)) << std::endl;
523 return EXIT_FAILURE;
524 }
525 */
526
527 //
528 // Run a bunch of tests for rank-1-updates, matrix-vector products
529 //
530 std::cout << "------------ Testing rank-1-updates and matrix-vector products ------------------" << std::endl;
531
532 std::cout << "* m = full, v1 = full, v2 = full" << std::endl;
533 retval = test_prod_rank1<NumericT>(epsilon,
534 std_m1, std_v1, std_v2,
535 vcl_m1_native, vcl_v1_native, vcl_v2_native);
536 if (retval == EXIT_FAILURE)
537 {
538 std::cout << " --- FAILED! ---" << std::endl;
539 return retval;
540 }
541 else
542 std::cout << " --- PASSED ---" << std::endl;
543
544
545 std::cout << "* m = full, v1 = full, v2 = range" << std::endl;
546 retval = test_prod_rank1<NumericT>(epsilon,
547 std_m1, std_v1, std_v2,
548 vcl_m1_native, vcl_v1_native, vcl_v2_range);
549 if (retval == EXIT_FAILURE)
550 {
551 std::cout << " --- FAILED! ---" << std::endl;
552 return retval;
553 }
554 else
555 std::cout << " --- PASSED ---" << std::endl;
556
557
558 std::cout << "* m = full, v1 = full, v2 = slice" << std::endl;
559 retval = test_prod_rank1<NumericT>(epsilon,
560 std_m1, std_v1, std_v2,
561 vcl_m1_native, vcl_v1_native, vcl_v2_slice);
562 if (retval == EXIT_FAILURE)
563 {
564 std::cout << " --- FAILED! ---" << std::endl;
565 return retval;
566 }
567 else
568 std::cout << " --- PASSED ---" << std::endl;
569
570
571 // v1 = range
572
573
574 std::cout << "* m = full, v1 = range, v2 = full" << std::endl;
575 retval = test_prod_rank1<NumericT>(epsilon,
576 std_m1, std_v1, std_v2,
577 vcl_m1_native, vcl_v1_range, vcl_v2_native);
578 if (retval == EXIT_FAILURE)
579 {
580 std::cout << " --- FAILED! ---" << std::endl;
581 return retval;
582 }
583 else
584 std::cout << " --- PASSED ---" << std::endl;
585
586
587 std::cout << "* m = full, v1 = range, v2 = range" << std::endl;
588 retval = test_prod_rank1<NumericT>(epsilon,
589 std_m1, std_v1, std_v2,
590 vcl_m1_native, vcl_v1_range, vcl_v2_range);
591 if (retval == EXIT_FAILURE)
592 {
593 std::cout << " --- FAILED! ---" << std::endl;
594 return retval;
595 }
596 else
597 std::cout << " --- PASSED ---" << std::endl;
598
599
600 std::cout << "* m = full, v1 = range, v2 = slice" << std::endl;
601 retval = test_prod_rank1<NumericT>(epsilon,
602 std_m1, std_v1, std_v2,
603 vcl_m1_native, vcl_v1_range, vcl_v2_slice);
604 if (retval == EXIT_FAILURE)
605 {
606 std::cout << " --- FAILED! ---" << std::endl;
607 return retval;
608 }
609 else
610 std::cout << " --- PASSED ---" << std::endl;
611
612
613
614 // v1 = slice
615
616 std::cout << "* m = full, v1 = slice, v2 = full" << std::endl;
617 retval = test_prod_rank1<NumericT>(epsilon,
618 std_m1, std_v1, std_v2,
619 vcl_m1_native, vcl_v1_slice, vcl_v2_native);
620 if (retval == EXIT_FAILURE)
621 {
622 std::cout << " --- FAILED! ---" << std::endl;
623 return retval;
624 }
625 else
626 std::cout << " --- PASSED ---" << std::endl;
627
628
629 std::cout << "* m = full, v1 = slice, v2 = range" << std::endl;
630 retval = test_prod_rank1<NumericT>(epsilon,
631 std_m1, std_v1, std_v2,
632 vcl_m1_native, vcl_v1_slice, vcl_v2_range);
633 if (retval == EXIT_FAILURE)
634 {
635 std::cout << " --- FAILED! ---" << std::endl;
636 return retval;
637 }
638 else
639 std::cout << " --- PASSED ---" << std::endl;
640
641
642 std::cout << "* m = full, v1 = slice, v2 = slice" << std::endl;
643 retval = test_prod_rank1<NumericT>(epsilon,
644 std_m1, std_v1, std_v2,
645 vcl_m1_native, vcl_v1_slice, vcl_v2_slice);
646 if (retval == EXIT_FAILURE)
647 {
648 std::cout << " --- FAILED! ---" << std::endl;
649 return retval;
650 }
651 else
652 std::cout << " --- PASSED ---" << std::endl;
653
654
655 ///////////////////////////// matrix_range
656
657 std::cout << "* m = range, v1 = full, v2 = full" << std::endl;
658 retval = test_prod_rank1<NumericT>(epsilon,
659 std_m1, std_v1, std_v2,
660 vcl_m1_range, vcl_v1_native, vcl_v2_native);
661 if (retval == EXIT_FAILURE)
662 {
663 std::cout << " --- FAILED! ---" << std::endl;
664 return retval;
665 }
666 else
667 std::cout << " --- PASSED ---" << std::endl;
668
669
670 std::cout << "* m = range, v1 = full, v2 = range" << std::endl;
671 retval = test_prod_rank1<NumericT>(epsilon,
672 std_m1, std_v1, std_v2,
673 vcl_m1_range, vcl_v1_native, vcl_v2_range);
674 if (retval == EXIT_FAILURE)
675 {
676 std::cout << " --- FAILED! ---" << std::endl;
677 return retval;
678 }
679 else
680 std::cout << " --- PASSED ---" << std::endl;
681
682
683 std::cout << "* m = range, v1 = full, v2 = slice" << std::endl;
684 retval = test_prod_rank1<NumericT>(epsilon,
685 std_m1, std_v1, std_v2,
686 vcl_m1_range, vcl_v1_native, vcl_v2_slice);
687 if (retval == EXIT_FAILURE)
688 {
689 std::cout << " --- FAILED! ---" << std::endl;
690 return retval;
691 }
692 else
693 std::cout << " --- PASSED ---" << std::endl;
694
695
696 // v1 = range
697
698
699 std::cout << "* m = range, v1 = range, v2 = full" << std::endl;
700 retval = test_prod_rank1<NumericT>(epsilon,
701 std_m1, std_v1, std_v2,
702 vcl_m1_range, vcl_v1_range, vcl_v2_native);
703 if (retval == EXIT_FAILURE)
704 {
705 std::cout << " --- FAILED! ---" << std::endl;
706 return retval;
707 }
708 else
709 std::cout << " --- PASSED ---" << std::endl;
710
711
712 std::cout << "* m = range, v1 = range, v2 = range" << std::endl;
713 retval = test_prod_rank1<NumericT>(epsilon,
714 std_m1, std_v1, std_v2,
715 vcl_m1_range, vcl_v1_range, vcl_v2_range);
716 if (retval == EXIT_FAILURE)
717 {
718 std::cout << " --- FAILED! ---" << std::endl;
719 return retval;
720 }
721 else
722 std::cout << " --- PASSED ---" << std::endl;
723
724
725 std::cout << "* m = range, v1 = range, v2 = slice" << std::endl;
726 retval = test_prod_rank1<NumericT>(epsilon,
727 std_m1, std_v1, std_v2,
728 vcl_m1_range, vcl_v1_range, vcl_v2_slice);
729 if (retval == EXIT_FAILURE)
730 {
731 std::cout << " --- FAILED! ---" << std::endl;
732 return retval;
733 }
734 else
735 std::cout << " --- PASSED ---" << std::endl;
736
737
738
739 // v1 = slice
740
741 std::cout << "* m = range, v1 = slice, v2 = full" << std::endl;
742 retval = test_prod_rank1<NumericT>(epsilon,
743 std_m1, std_v1, std_v2,
744 vcl_m1_range, vcl_v1_slice, vcl_v2_native);
745 if (retval == EXIT_FAILURE)
746 {
747 std::cout << " --- FAILED! ---" << std::endl;
748 return retval;
749 }
750 else
751 std::cout << " --- PASSED ---" << std::endl;
752
753
754 std::cout << "* m = range, v1 = slice, v2 = range" << std::endl;
755 retval = test_prod_rank1<NumericT>(epsilon,
756 std_m1, std_v1, std_v2,
757 vcl_m1_range, vcl_v1_slice, vcl_v2_range);
758 if (retval == EXIT_FAILURE)
759 {
760 std::cout << " --- FAILED! ---" << std::endl;
761 return retval;
762 }
763 else
764 std::cout << " --- PASSED ---" << std::endl;
765
766
767 std::cout << "* m = range, v1 = slice, v2 = slice" << std::endl;
768 retval = test_prod_rank1<NumericT>(epsilon,
769 std_m1, std_v1, std_v2,
770 vcl_m1_range, vcl_v1_slice, vcl_v2_slice);
771 if (retval == EXIT_FAILURE)
772 {
773 std::cout << " --- FAILED! ---" << std::endl;
774 return retval;
775 }
776 else
777 std::cout << " --- PASSED ---" << std::endl;
778
779
780 ///////////////////////////// matrix_slice
781
782 std::cout << "* m = slice, v1 = full, v2 = full" << std::endl;
783 retval = test_prod_rank1<NumericT>(epsilon,
784 std_m1, std_v1, std_v2,
785 vcl_m1_slice, vcl_v1_native, vcl_v2_native);
786 if (retval == EXIT_FAILURE)
787 {
788 std::cout << " --- FAILED! ---" << std::endl;
789 return retval;
790 }
791 else
792 std::cout << " --- PASSED ---" << std::endl;
793
794
795 std::cout << "* m = slice, v1 = full, v2 = range" << std::endl;
796 retval = test_prod_rank1<NumericT>(epsilon,
797 std_m1, std_v1, std_v2,
798 vcl_m1_slice, vcl_v1_native, vcl_v2_range);
799 if (retval == EXIT_FAILURE)
800 {
801 std::cout << " --- FAILED! ---" << std::endl;
802 return retval;
803 }
804 else
805 std::cout << " --- PASSED ---" << std::endl;
806
807
808 std::cout << "* m = slice, v1 = full, v2 = slice" << std::endl;
809 retval = test_prod_rank1<NumericT>(epsilon,
810 std_m1, std_v1, std_v2,
811 vcl_m1_slice, vcl_v1_native, vcl_v2_slice);
812 if (retval == EXIT_FAILURE)
813 {
814 std::cout << " --- FAILED! ---" << std::endl;
815 return retval;
816 }
817 else
818 std::cout << " --- PASSED ---" << std::endl;
819
820
821 // v1 = range
822
823
824 std::cout << "* m = slice, v1 = range, v2 = full" << std::endl;
825 retval = test_prod_rank1<NumericT>(epsilon,
826 std_m1, std_v1, std_v2,
827 vcl_m1_slice, vcl_v1_range, vcl_v2_native);
828 if (retval == EXIT_FAILURE)
829 {
830 std::cout << " --- FAILED! ---" << std::endl;
831 return retval;
832 }
833 else
834 std::cout << " --- PASSED ---" << std::endl;
835
836
837 std::cout << "* m = slice, v1 = range, v2 = range" << std::endl;
838 retval = test_prod_rank1<NumericT>(epsilon,
839 std_m1, std_v1, std_v2,
840 vcl_m1_slice, vcl_v1_range, vcl_v2_range);
841 if (retval == EXIT_FAILURE)
842 {
843 std::cout << " --- FAILED! ---" << std::endl;
844 return retval;
845 }
846 else
847 std::cout << " --- PASSED ---" << std::endl;
848
849
850 std::cout << "* m = slice, v1 = range, v2 = slice" << std::endl;
851 retval = test_prod_rank1<NumericT>(epsilon,
852 std_m1, std_v1, std_v2,
853 vcl_m1_slice, vcl_v1_range, vcl_v2_slice);
854 if (retval == EXIT_FAILURE)
855 {
856 std::cout << " --- FAILED! ---" << std::endl;
857 return retval;
858 }
859 else
860 std::cout << " --- PASSED ---" << std::endl;
861
862
863
864 // v1 = slice
865
866 std::cout << "* m = slice, v1 = slice, v2 = full" << std::endl;
867 retval = test_prod_rank1<NumericT>(epsilon,
868 std_m1, std_v1, std_v2,
869 vcl_m1_slice, vcl_v1_slice, vcl_v2_native);
870 if (retval == EXIT_FAILURE)
871 {
872 std::cout << " --- FAILED! ---" << std::endl;
873 return retval;
874 }
875 else
876 std::cout << " --- PASSED ---" << std::endl;
877
878
879 std::cout << "* m = slice, v1 = slice, v2 = range" << std::endl;
880 retval = test_prod_rank1<NumericT>(epsilon,
881 std_m1, std_v1, std_v2,
882 vcl_m1_slice, vcl_v1_slice, vcl_v2_range);
883 if (retval == EXIT_FAILURE)
884 {
885 std::cout << " --- FAILED! ---" << std::endl;
886 return retval;
887 }
888 else
889 std::cout << " --- PASSED ---" << std::endl;
890
891
892 std::cout << "* m = slice, v1 = slice, v2 = slice" << std::endl;
893 retval = test_prod_rank1<NumericT>(epsilon,
894 std_m1, std_v1, std_v2,
895 vcl_m1_slice, vcl_v1_slice, vcl_v2_slice);
896 if (retval == EXIT_FAILURE)
897 {
898 std::cout << " --- FAILED! ---" << std::endl;
899 return retval;
900 }
901 else
902 std::cout << " --- PASSED ---" << std::endl;
903
904 return retval;
905 }
906 //
907 // -------------------------------------------------------------
908 //
main()909 int main()
910 {
911 std::cout << std::endl;
912 std::cout << "----------------------------------------------" << std::endl;
913 std::cout << "----------------------------------------------" << std::endl;
914 std::cout << "## Test :: Matrix" << std::endl;
915 std::cout << "----------------------------------------------" << std::endl;
916 std::cout << "----------------------------------------------" << std::endl;
917 std::cout << std::endl;
918
919 int retval = EXIT_SUCCESS;
920
921 std::cout << std::endl;
922 std::cout << "----------------------------------------------" << std::endl;
923 std::cout << std::endl;
924 {
925 typedef float NumericT;
926 NumericT epsilon = NumericT(1.0E-3);
927 std::cout << "# Testing setup:" << std::endl;
928 std::cout << " eps: " << epsilon << std::endl;
929 std::cout << " numeric: float" << std::endl;
930 std::cout << " layout: row-major" << std::endl;
931 retval = test<NumericT, viennacl::row_major>(epsilon);
932 if ( retval == EXIT_SUCCESS )
933 std::cout << "# Test passed" << std::endl;
934 else
935 return retval;
936 }
937 std::cout << std::endl;
938 std::cout << "----------------------------------------------" << std::endl;
939 std::cout << std::endl;
940 {
941 typedef float NumericT;
942 NumericT epsilon = NumericT(1.0E-3);
943 std::cout << "# Testing setup:" << std::endl;
944 std::cout << " eps: " << epsilon << std::endl;
945 std::cout << " numeric: float" << std::endl;
946 std::cout << " layout: column-major" << std::endl;
947 retval = test<NumericT, viennacl::column_major>(epsilon);
948 if ( retval == EXIT_SUCCESS )
949 std::cout << "# Test passed" << std::endl;
950 else
951 return retval;
952 }
953 std::cout << std::endl;
954 std::cout << "----------------------------------------------" << std::endl;
955 std::cout << std::endl;
956
957
958 #ifdef VIENNACL_WITH_OPENCL
959 if ( viennacl::ocl::current_device().double_support() )
960 #endif
961 {
962 {
963 typedef double NumericT;
964 NumericT epsilon = 1.0E-11;
965 std::cout << "# Testing setup:" << std::endl;
966 std::cout << " eps: " << epsilon << std::endl;
967 std::cout << " numeric: double" << std::endl;
968 std::cout << " layout: row-major" << std::endl;
969 retval = test<NumericT, viennacl::row_major>(epsilon);
970 if ( retval == EXIT_SUCCESS )
971 std::cout << "# Test passed" << std::endl;
972 else
973 return retval;
974 }
975 std::cout << std::endl;
976 std::cout << "----------------------------------------------" << std::endl;
977 std::cout << std::endl;
978 {
979 typedef double NumericT;
980 NumericT epsilon = 1.0E-11;
981 std::cout << "# Testing setup:" << std::endl;
982 std::cout << " eps: " << epsilon << std::endl;
983 std::cout << " numeric: double" << std::endl;
984 std::cout << " layout: column-major" << std::endl;
985 retval = test<NumericT, viennacl::column_major>(epsilon);
986 if ( retval == EXIT_SUCCESS )
987 std::cout << "# Test passed" << std::endl;
988 else
989 return retval;
990 }
991 std::cout << std::endl;
992 std::cout << "----------------------------------------------" << std::endl;
993 std::cout << std::endl;
994 }
995
996 std::cout << std::endl;
997 std::cout << "------- Test completed --------" << std::endl;
998 std::cout << std::endl;
999
1000
1001 return retval;
1002 }
1003