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_vector.cpp Tests the scheduler for vector-operations.
21 * \test Tests the scheduler for vector-operations.
22 **/
23
24 //
25 // *** System
26 //
27 #include <iostream>
28 #include <iomanip>
29 #include <vector>
30
31 //
32 // *** ViennaCL
33 //
34 #include "viennacl/vector.hpp"
35 #include "viennacl/matrix.hpp"
36 #include "viennacl/vector_proxy.hpp"
37 #include "viennacl/linalg/inner_prod.hpp"
38 #include "viennacl/linalg/norm_1.hpp"
39 #include "viennacl/linalg/norm_2.hpp"
40 #include "viennacl/linalg/norm_inf.hpp"
41
42 #include "viennacl/scheduler/execute.hpp"
43 #include "viennacl/scheduler/io.hpp"
44
45 #include "viennacl/tools/random.hpp"
46
47
48 //
49 // -------------------------------------------------------------
50 //
51 template<typename ScalarType>
diff(ScalarType const & s1,ScalarType const & s2)52 ScalarType diff(ScalarType const & s1, ScalarType const & s2)
53 {
54 viennacl::backend::finish();
55 if (std::fabs(s1 - s2) > 0)
56 return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
57 return 0;
58 }
59 //
60 // -------------------------------------------------------------
61 //
62 template<typename ScalarType>
diff(ScalarType const & s1,viennacl::scalar<ScalarType> const & s2)63 ScalarType diff(ScalarType const & s1, viennacl::scalar<ScalarType> const & s2)
64 {
65 viennacl::backend::finish();
66 if (std::fabs(s1 - s2) > 0)
67 return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
68 return 0;
69 }
70 //
71 // -------------------------------------------------------------
72 //
73 template<typename ScalarType>
diff(ScalarType const & s1,viennacl::entry_proxy<ScalarType> const & s2)74 ScalarType diff(ScalarType const & s1, viennacl::entry_proxy<ScalarType> const & s2)
75 {
76 viennacl::backend::finish();
77 if (std::fabs(s1 - s2) > 0)
78 return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
79 return 0;
80 }
81 //
82 // -------------------------------------------------------------
83 //
84 template<typename ScalarType, typename ViennaCLVectorType>
85 ScalarType diff(std::vector<ScalarType> const & v1, ViennaCLVectorType const & vcl_vec)
86 {
87 std::vector<ScalarType> v2_cpu(vcl_vec.size());
88 viennacl::backend::finish();
89 viennacl::copy(vcl_vec, v2_cpu);
90
91 ScalarType norm_inf_value = 0;
92 for (std::size_t i=0;i<v1.size(); ++i)
93 {
94 ScalarType tmp = 0;
95 if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
96 tmp = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
97
98 norm_inf_value = (tmp > norm_inf_value) ? tmp : norm_inf_value;
99 }
100
101 return norm_inf_value;
102 }
103
104
105 template<typename T1, typename T2>
check(T1 const & t1,T2 const & t2,double epsilon)106 int check(T1 const & t1, T2 const & t2, double epsilon)
107 {
108 int retval = EXIT_SUCCESS;
109
110 double temp = std::fabs(diff(t1, t2));
111 if (temp > epsilon)
112 {
113 std::cout << "# Error! Relative difference: " << temp << std::endl;
114 retval = EXIT_FAILURE;
115 }
116 else
117 std::cout << "PASSED!" << std::endl;
118 return retval;
119 }
120
121
122 //
123 // -------------------------------------------------------------
124 //
125 template< typename NumericT, typename Epsilon, typename STLVectorType, typename ViennaCLVectorType1, typename ViennaCLVectorType2 >
test(Epsilon const & epsilon,STLVectorType & std_v1,STLVectorType & std_v2,ViennaCLVectorType1 & vcl_v1,ViennaCLVectorType2 & vcl_v2)126 int test(Epsilon const& epsilon,
127 STLVectorType & std_v1, STLVectorType & std_v2,
128 ViennaCLVectorType1 & vcl_v1, ViennaCLVectorType2 & vcl_v2)
129 {
130 int retval = EXIT_SUCCESS;
131
132 viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
133
134 NumericT cpu_result = 42.0;
135 viennacl::scalar<NumericT> gpu_result = 43.0;
136 NumericT alpha = NumericT(3.1415);
137 NumericT beta = NumericT(2.7172);
138
139 //
140 // Initializer:
141 //
142 std::cout << "Checking for zero_vector initializer..." << std::endl;
143 std_v1 = std::vector<NumericT>(std_v1.size());
144 vcl_v1 = viennacl::zero_vector<NumericT>(vcl_v1.size());
145 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
146 return EXIT_FAILURE;
147
148 std::cout << "Checking for scalar_vector initializer..." << std::endl;
149 std_v1 = std::vector<NumericT>(std_v1.size(), cpu_result);
150 vcl_v1 = viennacl::scalar_vector<NumericT>(vcl_v1.size(), cpu_result);
151 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
152 return EXIT_FAILURE;
153
154 std_v1 = std::vector<NumericT>(std_v1.size(), gpu_result);
155 vcl_v1 = viennacl::scalar_vector<NumericT>(vcl_v1.size(), gpu_result);
156 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
157 return EXIT_FAILURE;
158
159 std::cout << "Checking for unit_vector initializer..." << std::endl;
160 std_v1 = std::vector<NumericT>(std_v1.size()); std_v1[5] = NumericT(1);
161 vcl_v1 = viennacl::unit_vector<NumericT>(vcl_v1.size(), 5);
162 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
163 return EXIT_FAILURE;
164
165
166 for (std::size_t i=0; i<std_v1.size(); ++i)
167 {
168 std_v1[i] = NumericT(1.0) + randomNumber();
169 std_v2[i] = NumericT(1.0) + randomNumber();
170 }
171
172 viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin()); //resync
173 viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
174
175 std::cout << "Checking for successful copy..." << std::endl;
176 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
177 return EXIT_FAILURE;
178 if (check(std_v2, vcl_v2, epsilon) != EXIT_SUCCESS)
179 return EXIT_FAILURE;
180
181
182 // --------------------------------------------------------------------------
183
184 std::cout << "Testing simple assignments..." << std::endl;
185
186 {
187 for (std::size_t i=0; i<std_v1.size(); ++i)
188 std_v1[i] = std_v2[i];
189 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), vcl_v2); // same as vcl_v1 = vcl_v2;
190 viennacl::scheduler::execute(my_statement);
191
192 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
193 return EXIT_FAILURE;
194 }
195
196 {
197 for (std::size_t i=0; i<std_v1.size(); ++i)
198 std_v1[i] += std_v2[i];
199 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_add(), vcl_v2); // same as vcl_v1 += vcl_v2;
200 viennacl::scheduler::execute(my_statement);
201
202 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
203 return EXIT_FAILURE;
204 }
205
206 {
207 for (std::size_t i=0; i<std_v1.size(); ++i)
208 std_v1[i] -= std_v2[i];
209 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_sub(), vcl_v2); // same as vcl_v1 -= vcl_v2;
210 viennacl::scheduler::execute(my_statement);
211
212 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
213 return EXIT_FAILURE;
214 }
215
216 std::cout << "Testing composite assignments..." << std::endl;
217 {
218 for (std::size_t i=0; i<std_v1.size(); ++i)
219 std_v1[i] = std_v1[i] + std_v2[i];
220 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), vcl_v1 + vcl_v2); // same as vcl_v1 = vcl_v1 + vcl_v2;
221 viennacl::scheduler::execute(my_statement);
222
223 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
224 return EXIT_FAILURE;
225 }
226 {
227 for (std::size_t i=0; i<std_v1.size(); ++i)
228 std_v1[i] += alpha * std_v1[i] - beta * std_v2[i] + std_v1[i] / beta - std_v2[i] / alpha;
229 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_add(), alpha * vcl_v1 - beta * vcl_v2 + vcl_v1 / beta - vcl_v2 / alpha); // same as vcl_v1 += alpha * vcl_v1 - beta * vcl_v2 + beta * vcl_v1 - alpha * vcl_v2;
230 viennacl::scheduler::execute(my_statement);
231
232 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
233 return EXIT_FAILURE;
234 }
235
236 {
237 for (std::size_t i=0; i<std_v1.size(); ++i)
238 std_v1[i] = std_v1[i] - std_v2[i];
239 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), vcl_v1 - vcl_v2); // same as vcl_v1 = vcl_v1 - vcl_v2;
240 viennacl::scheduler::execute(my_statement);
241
242 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
243 return EXIT_FAILURE;
244 }
245
246 std::cout << "--- Testing reductions ---" << std::endl;
247 std::cout << "inner_prod..." << std::endl;
248 {
249 cpu_result = 0;
250 for (std::size_t i=0; i<std_v1.size(); ++i)
251 cpu_result += std_v1[i] * std_v2[i];
252 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::inner_prod(vcl_v1, vcl_v2)); // same as gpu_result = inner_prod(vcl_v1, vcl_v2);
253 viennacl::scheduler::execute(my_statement);
254
255 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
256 return EXIT_FAILURE;
257 }
258
259 {
260 cpu_result = 0;
261 for (std::size_t i=0; i<std_v1.size(); ++i)
262 cpu_result += (std_v1[i] + std_v2[i]) * std_v2[i];
263 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::inner_prod(vcl_v1 + vcl_v2, vcl_v2)); // same as gpu_result = inner_prod(vcl_v1 + vcl_v2, vcl_v2);
264 viennacl::scheduler::execute(my_statement);
265
266 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
267 return EXIT_FAILURE;
268 }
269
270 {
271 cpu_result = 0;
272 for (std::size_t i=0; i<std_v1.size(); ++i)
273 cpu_result += std_v1[i] * (std_v2[i] - std_v1[i]);
274 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::inner_prod(vcl_v1, vcl_v2 - vcl_v1)); // same as gpu_result = inner_prod(vcl_v1, vcl_v2 - vcl_v1);
275 viennacl::scheduler::execute(my_statement);
276
277 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
278 return EXIT_FAILURE;
279 }
280
281 {
282 cpu_result = 0;
283 for (std::size_t i=0; i<std_v1.size(); ++i)
284 cpu_result += (std_v1[i] - std_v2[i]) * (std_v2[i] + std_v1[i]);
285 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::inner_prod(vcl_v1 - vcl_v2, vcl_v2 + vcl_v1)); // same as gpu_result = inner_prod(vcl_v1 - vcl_v2, vcl_v2 + vcl_v1);
286 viennacl::scheduler::execute(my_statement);
287
288 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
289 return EXIT_FAILURE;
290 }
291
292 std::cout << "norm_1..." << std::endl;
293 {
294 cpu_result = 0;
295 for (std::size_t i=0; i<std_v1.size(); ++i)
296 cpu_result += std::fabs(std_v1[i]);
297 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_1(vcl_v1)); // same as gpu_result = norm_1(vcl_v1);
298 viennacl::scheduler::execute(my_statement);
299
300 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
301 return EXIT_FAILURE;
302 }
303
304 {
305 cpu_result = 0;
306 for (std::size_t i=0; i<std_v1.size(); ++i)
307 cpu_result += std::fabs(std_v1[i] + std_v2[i]);
308 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_1(vcl_v1 + vcl_v2)); // same as gpu_result = norm_1(vcl_v1 + vcl_v2);
309 viennacl::scheduler::execute(my_statement);
310
311 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
312 return EXIT_FAILURE;
313 }
314
315 std::cout << "norm_2..." << std::endl;
316 {
317 cpu_result = 0;
318 for (std::size_t i=0; i<std_v1.size(); ++i)
319 cpu_result += std_v1[i] * std_v1[i];
320 cpu_result = std::sqrt(cpu_result);
321 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_2(vcl_v1)); // same as gpu_result = norm_2(vcl_v1);
322 viennacl::scheduler::execute(my_statement);
323
324 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
325 return EXIT_FAILURE;
326 }
327
328 {
329 cpu_result = 0;
330 for (std::size_t i=0; i<std_v1.size(); ++i)
331 cpu_result += (std_v1[i] + std_v2[i]) * (std_v1[i] + std_v2[i]);
332 cpu_result = std::sqrt(cpu_result);
333 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_2(vcl_v1 + vcl_v2)); // same as gpu_result = norm_2(vcl_v1 + vcl_v2);
334 viennacl::scheduler::execute(my_statement);
335
336 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
337 return EXIT_FAILURE;
338 }
339
340 std::cout << "norm_inf..." << std::endl;
341 {
342 cpu_result = 0;
343 for (std::size_t i=0; i<std_v1.size(); ++i)
344 cpu_result = std::max(cpu_result, std::fabs(std_v1[i]));
345 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_inf(vcl_v1)); // same as gpu_result = norm_inf(vcl_v1);
346 viennacl::scheduler::execute(my_statement);
347
348 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
349 return EXIT_FAILURE;
350 }
351
352 {
353 cpu_result = 0;
354 for (std::size_t i=0; i<std_v1.size(); ++i)
355 cpu_result = std::max(cpu_result, std::fabs(std_v1[i] - std_v2[i]));
356 viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_inf(vcl_v1 - vcl_v2)); // same as gpu_result = norm_inf(vcl_v1 - vcl_v2);
357 viennacl::scheduler::execute(my_statement);
358
359 if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
360 return EXIT_FAILURE;
361 }
362
363 std::cout << "--- Testing elementwise operations (binary) ---" << std::endl;
364 std::cout << "x = element_prod(x, y)... ";
365 {
366 for (std::size_t i=0; i<std_v1.size(); ++i)
367 std_v1[i] = std_v1[i] * std_v2[i];
368 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_v1, vcl_v2));
369 viennacl::scheduler::execute(my_statement);
370
371 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
372 return EXIT_FAILURE;
373 }
374
375 std::cout << "x = element_prod(x + y, y)... ";
376 {
377 for (std::size_t i=0; i<std_v1.size(); ++i)
378 std_v1[i] = (std_v1[i] + std_v2[i]) * std_v2[i];
379 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_v1 + vcl_v2, vcl_v2));
380 viennacl::scheduler::execute(my_statement);
381
382 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
383 return EXIT_FAILURE;
384 }
385
386 std::cout << "x = element_prod(x, x + y)... ";
387 {
388 for (std::size_t i=0; i<std_v1.size(); ++i)
389 std_v1[i] = std_v1[i] * (std_v1[i] + std_v2[i]);
390 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_v1, vcl_v2 + vcl_v1));
391 viennacl::scheduler::execute(my_statement);
392
393 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
394 return EXIT_FAILURE;
395 }
396
397 std::cout << "x = element_prod(x - y, y + x)... ";
398 {
399 for (std::size_t i=0; i<std_v1.size(); ++i)
400 std_v1[i] = (std_v1[i] - std_v2[i]) * (std_v2[i] + std_v1[i]);
401 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_v1 - vcl_v2, vcl_v2 + vcl_v1));
402 viennacl::scheduler::execute(my_statement);
403
404 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
405 return EXIT_FAILURE;
406 }
407
408
409
410 std::cout << "x = element_div(x, y)... ";
411 {
412 for (std::size_t i=0; i<std_v1.size(); ++i)
413 std_v1[i] = std_v1[i] / std_v2[i];
414 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_div(vcl_v1, vcl_v2));
415 viennacl::scheduler::execute(my_statement);
416
417 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
418 return EXIT_FAILURE;
419 }
420
421 std::cout << "x = element_div(x + y, y)... ";
422 {
423 for (std::size_t i=0; i<std_v1.size(); ++i)
424 std_v1[i] = (std_v1[i] + std_v2[i]) / std_v2[i];
425 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_div(vcl_v1 + vcl_v2, vcl_v2));
426 viennacl::scheduler::execute(my_statement);
427
428 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
429 return EXIT_FAILURE;
430 }
431
432 std::cout << "x = element_div(x, x + y)... ";
433 {
434 for (std::size_t i=0; i<std_v1.size(); ++i)
435 std_v1[i] = std_v1[i] / (std_v1[i] + std_v2[i]);
436 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_div(vcl_v1, vcl_v2 + vcl_v1));
437 viennacl::scheduler::execute(my_statement);
438
439 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
440 return EXIT_FAILURE;
441 }
442
443 std::cout << "x = element_div(x - y, y + x)... ";
444 {
445 for (std::size_t i=0; i<std_v1.size(); ++i)
446 std_v1[i] = (std_v1[i] - std_v2[i]) / (std_v2[i] + std_v1[i]);
447 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_div(vcl_v1 - vcl_v2, vcl_v2 + vcl_v1));
448 viennacl::scheduler::execute(my_statement);
449
450 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
451 return EXIT_FAILURE;
452 }
453
454
455 std::cout << "x = element_pow(x, y)... ";
456 {
457 for (std::size_t i=0; i<std_v1.size(); ++i)
458 {
459 std_v1[i] = NumericT(2.0) + randomNumber();
460 std_v2[i] = NumericT(1.0) + randomNumber();
461 }
462 viennacl::copy(std_v1, vcl_v1);
463 viennacl::copy(std_v2, vcl_v2);
464
465 for (std::size_t i=0; i<std_v1.size(); ++i)
466 std_v1[i] = std::pow(std_v1[i], std_v2[i]);
467 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_pow(vcl_v1, vcl_v2));
468 viennacl::scheduler::execute(my_statement);
469
470 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
471 return EXIT_FAILURE;
472 }
473
474 std::cout << "x = element_pow(x + y, y)... ";
475 {
476 for (std::size_t i=0; i<std_v1.size(); ++i)
477 {
478 std_v1[i] = NumericT(2.0) + randomNumber();
479 std_v2[i] = NumericT(1.0) + randomNumber();
480 }
481 viennacl::copy(std_v1, vcl_v1);
482 viennacl::copy(std_v2, vcl_v2);
483
484 for (std::size_t i=0; i<std_v1.size(); ++i)
485 std_v1[i] = std::pow(std_v1[i] + std_v2[i], std_v2[i]);
486 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_pow(vcl_v1 + vcl_v2, vcl_v2));
487 viennacl::scheduler::execute(my_statement);
488
489 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
490 return EXIT_FAILURE;
491 }
492
493 std::cout << "x = element_pow(x, x + y)... ";
494 {
495 for (std::size_t i=0; i<std_v1.size(); ++i)
496 {
497 std_v1[i] = NumericT(2.0) + randomNumber();
498 std_v2[i] = NumericT(1.0) + randomNumber();
499 }
500 viennacl::copy(std_v1, vcl_v1);
501 viennacl::copy(std_v2, vcl_v2);
502
503 for (std::size_t i=0; i<std_v1.size(); ++i)
504 std_v1[i] = std::pow(std_v1[i], std_v1[i] + std_v2[i]);
505 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_pow(vcl_v1, vcl_v2 + vcl_v1));
506 viennacl::scheduler::execute(my_statement);
507
508 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
509 return EXIT_FAILURE;
510 }
511
512 std::cout << "x = element_pow(x - y, y + x)... ";
513 {
514 for (std::size_t i=0; i<std_v1.size(); ++i)
515 {
516 std_v1[i] = NumericT(2.0) + randomNumber();
517 std_v2[i] = NumericT(1.0) + randomNumber();
518 }
519 viennacl::copy(std_v1, vcl_v1);
520 viennacl::copy(std_v2, vcl_v2);
521
522 for (std::size_t i=0; i<std_v1.size(); ++i)
523 std_v1[i] = std::pow(std_v1[i] - std_v2[i], std_v2[i] + std_v1[i]);
524 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_pow(vcl_v1 - vcl_v2, vcl_v2 + vcl_v1));
525 viennacl::scheduler::execute(my_statement);
526
527 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
528 return EXIT_FAILURE;
529 }
530
531 std::cout << "--- Testing elementwise operations (unary) ---" << std::endl;
532 #define GENERATE_UNARY_OP_TEST(OPNAME) \
533 std_v1 = std::vector<NumericT>(std_v1.size(), NumericT(0.21)); \
534 for (std::size_t i=0; i<std_v1.size(); ++i) \
535 std_v2[i] = NumericT(3.1415) * std_v1[i]; \
536 viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin()); \
537 viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin()); \
538 { \
539 for (std::size_t i=0; i<std_v1.size(); ++i) \
540 std_v1[i] = std::OPNAME(std_v2[i]); \
541 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_##OPNAME(vcl_v2)); \
542 viennacl::scheduler::execute(my_statement); \
543 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS) \
544 return EXIT_FAILURE; \
545 } \
546 { \
547 for (std::size_t i=0; i<std_v1.size(); ++i) \
548 std_v1[i] = std::OPNAME(std_v2[i] / NumericT(2)); \
549 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_##OPNAME(vcl_v2 / NumericT(2))); \
550 viennacl::scheduler::execute(my_statement); \
551 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS) \
552 return EXIT_FAILURE; \
553 }
554
555 GENERATE_UNARY_OP_TEST(cos);
556 GENERATE_UNARY_OP_TEST(cosh);
557 GENERATE_UNARY_OP_TEST(exp);
558 GENERATE_UNARY_OP_TEST(floor);
559 GENERATE_UNARY_OP_TEST(fabs);
560 GENERATE_UNARY_OP_TEST(log);
561 GENERATE_UNARY_OP_TEST(log10);
562 GENERATE_UNARY_OP_TEST(sin);
563 GENERATE_UNARY_OP_TEST(sinh);
564 GENERATE_UNARY_OP_TEST(fabs);
565 //GENERATE_UNARY_OP_TEST(abs); //OpenCL allows abs on integers only
566 GENERATE_UNARY_OP_TEST(sqrt);
567 GENERATE_UNARY_OP_TEST(tan);
568 GENERATE_UNARY_OP_TEST(tanh);
569
570 #undef GENERATE_UNARY_OP_TEST
571
572 std::cout << "--- Testing complicated composite operations ---" << std::endl;
573 std::cout << "x = inner_prod(x, y) * y..." << std::endl;
574 {
575 cpu_result = 0;
576 for (std::size_t i=0; i<std_v1.size(); ++i)
577 cpu_result += std_v1[i] * std_v2[i];
578 for (std::size_t i=0; i<std_v1.size(); ++i)
579 std_v1[i] = cpu_result * std_v2[i];
580 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::inner_prod(vcl_v1, vcl_v2) * vcl_v2);
581 viennacl::scheduler::execute(my_statement);
582
583 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
584 return EXIT_FAILURE;
585 }
586
587 std::cout << "x = y / norm_1(x)..." << std::endl;
588 {
589 cpu_result = 0;
590 for (std::size_t i=0; i<std_v1.size(); ++i)
591 cpu_result += std::fabs(std_v1[i]);
592 for (std::size_t i=0; i<std_v1.size(); ++i)
593 std_v1[i] = std_v2[i] / cpu_result;
594 viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), vcl_v2 / viennacl::linalg::norm_1(vcl_v1) );
595 viennacl::scheduler::execute(my_statement);
596
597 if (check(std_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
598 return EXIT_FAILURE;
599 }
600
601
602 // --------------------------------------------------------------------------
603 return retval;
604 }
605
606
607 template< typename NumericT, typename Epsilon >
test(Epsilon const & epsilon)608 int test(Epsilon const& epsilon)
609 {
610 int retval = EXIT_SUCCESS;
611 std::size_t size = 24656;
612
613 viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
614
615 std::cout << "Running tests for vector of size " << size << std::endl;
616
617 //
618 // Set up STL objects
619 //
620 std::vector<NumericT> std_full_vec(size);
621 std::vector<NumericT> std_full_vec2(std_full_vec.size());
622
623 for (std::size_t i=0; i<std_full_vec.size(); ++i)
624 {
625 std_full_vec[i] = NumericT(1.0) + randomNumber();
626 std_full_vec2[i] = NumericT(1.0) + randomNumber();
627 }
628
629 std::vector<NumericT> std_range_vec (2 * std_full_vec.size() / 4 - std_full_vec.size() / 4);
630 std::vector<NumericT> std_range_vec2(2 * std_full_vec.size() / 4 - std_full_vec.size() / 4);
631
632 for (std::size_t i=0; i<std_range_vec.size(); ++i)
633 std_range_vec[i] = std_full_vec[i + std_full_vec.size() / 4];
634 for (std::size_t i=0; i<std_range_vec2.size(); ++i)
635 std_range_vec2[i] = std_full_vec2[i + 2 * std_full_vec2.size() / 4];
636
637 std::vector<NumericT> std_slice_vec (std_full_vec.size() / 4);
638 std::vector<NumericT> std_slice_vec2(std_full_vec.size() / 4);
639
640 for (std::size_t i=0; i<std_slice_vec.size(); ++i)
641 std_slice_vec[i] = std_full_vec[3*i + std_full_vec.size() / 4];
642 for (std::size_t i=0; i<std_slice_vec2.size(); ++i)
643 std_slice_vec2[i] = std_full_vec2[2*i + 2 * std_full_vec2.size() / 4];
644
645 //
646 // Set up ViennaCL objects
647 //
648 viennacl::vector<NumericT> vcl_full_vec(std_full_vec.size());
649 viennacl::vector<NumericT> vcl_full_vec2(std_full_vec2.size());
650
651 viennacl::fast_copy(std_full_vec.begin(), std_full_vec.end(), vcl_full_vec.begin());
652 viennacl::copy(std_full_vec2.begin(), std_full_vec2.end(), vcl_full_vec2.begin());
653
654 viennacl::range vcl_r1( vcl_full_vec.size() / 4, 2 * vcl_full_vec.size() / 4);
655 viennacl::range vcl_r2(2 * vcl_full_vec2.size() / 4, 3 * vcl_full_vec2.size() / 4);
656 viennacl::vector_range< viennacl::vector<NumericT> > vcl_range_vec(vcl_full_vec, vcl_r1);
657 viennacl::vector_range< viennacl::vector<NumericT> > vcl_range_vec2(vcl_full_vec2, vcl_r2);
658
659 {
660 viennacl::vector<NumericT> vcl_short_vec(vcl_range_vec);
661 viennacl::vector<NumericT> vcl_short_vec2 = vcl_range_vec2;
662
663 std::vector<NumericT> std_short_vec(std_range_vec);
664 std::vector<NumericT> std_short_vec2(std_range_vec2);
665
666 std::cout << "Testing creation of vectors from range..." << std::endl;
667 if (check(std_short_vec, vcl_short_vec, epsilon) != EXIT_SUCCESS)
668 return EXIT_FAILURE;
669 if (check(std_short_vec2, vcl_short_vec2, epsilon) != EXIT_SUCCESS)
670 return EXIT_FAILURE;
671 }
672
673 viennacl::slice vcl_s1( vcl_full_vec.size() / 4, 3, vcl_full_vec.size() / 4);
674 viennacl::slice vcl_s2(2 * vcl_full_vec2.size() / 4, 2, vcl_full_vec2.size() / 4);
675 viennacl::vector_slice< viennacl::vector<NumericT> > vcl_slice_vec(vcl_full_vec, vcl_s1);
676 viennacl::vector_slice< viennacl::vector<NumericT> > vcl_slice_vec2(vcl_full_vec2, vcl_s2);
677
678 viennacl::vector<NumericT> vcl_short_vec(vcl_slice_vec);
679 viennacl::vector<NumericT> vcl_short_vec2 = vcl_slice_vec2;
680
681 std::vector<NumericT> std_short_vec(std_slice_vec);
682 std::vector<NumericT> std_short_vec2(std_slice_vec2);
683
684 std::cout << "Testing creation of vectors from slice..." << std::endl;
685 if (check(std_short_vec, vcl_short_vec, epsilon) != EXIT_SUCCESS)
686 return EXIT_FAILURE;
687 if (check(std_short_vec2, vcl_short_vec2, epsilon) != EXIT_SUCCESS)
688 return EXIT_FAILURE;
689
690
691 //
692 // Now start running tests for vectors, ranges and slices:
693 //
694
695 std::cout << " ** vcl_v1 = vector, vcl_v2 = vector **" << std::endl;
696 retval = test<NumericT>(epsilon,
697 std_short_vec, std_short_vec2,
698 vcl_short_vec, vcl_short_vec2);
699 if (retval != EXIT_SUCCESS)
700 return EXIT_FAILURE;
701
702 std::cout << " ** vcl_v1 = vector, vcl_v2 = range **" << std::endl;
703 retval = test<NumericT>(epsilon,
704 std_short_vec, std_short_vec2,
705 vcl_short_vec, vcl_range_vec2);
706 if (retval != EXIT_SUCCESS)
707 return EXIT_FAILURE;
708
709 std::cout << " ** vcl_v1 = vector, vcl_v2 = slice **" << std::endl;
710 retval = test<NumericT>(epsilon,
711 std_short_vec, std_short_vec2,
712 vcl_short_vec, vcl_slice_vec2);
713 if (retval != EXIT_SUCCESS)
714 return EXIT_FAILURE;
715
716 ///////
717
718 std::cout << " ** vcl_v1 = range, vcl_v2 = vector **" << std::endl;
719 retval = test<NumericT>(epsilon,
720 std_short_vec, std_short_vec2,
721 vcl_range_vec, vcl_short_vec2);
722 if (retval != EXIT_SUCCESS)
723 return EXIT_FAILURE;
724
725 std::cout << " ** vcl_v1 = range, vcl_v2 = range **" << std::endl;
726 retval = test<NumericT>(epsilon,
727 std_short_vec, std_short_vec2,
728 vcl_range_vec, vcl_range_vec2);
729 if (retval != EXIT_SUCCESS)
730 return EXIT_FAILURE;
731
732 std::cout << " ** vcl_v1 = range, vcl_v2 = slice **" << std::endl;
733 retval = test<NumericT>(epsilon,
734 std_short_vec, std_short_vec2,
735 vcl_range_vec, vcl_slice_vec2);
736 if (retval != EXIT_SUCCESS)
737 return EXIT_FAILURE;
738
739 ///////
740
741 std::cout << " ** vcl_v1 = slice, vcl_v2 = vector **" << std::endl;
742 retval = test<NumericT>(epsilon,
743 std_short_vec, std_short_vec2,
744 vcl_slice_vec, vcl_short_vec2);
745 if (retval != EXIT_SUCCESS)
746 return EXIT_FAILURE;
747
748 std::cout << " ** vcl_v1 = slice, vcl_v2 = range **" << std::endl;
749 retval = test<NumericT>(epsilon,
750 std_short_vec, std_short_vec2,
751 vcl_slice_vec, vcl_range_vec2);
752 if (retval != EXIT_SUCCESS)
753 return EXIT_FAILURE;
754
755 std::cout << " ** vcl_v1 = slice, vcl_v2 = slice **" << std::endl;
756 retval = test<NumericT>(epsilon,
757 std_short_vec, std_short_vec2,
758 vcl_slice_vec, vcl_slice_vec2);
759 if (retval != EXIT_SUCCESS)
760 return EXIT_FAILURE;
761
762 return EXIT_SUCCESS;
763 }
764
765
766
767 //
768 // -------------------------------------------------------------
769 //
main()770 int main()
771 {
772 std::cout << std::endl;
773 std::cout << "----------------------------------------------" << std::endl;
774 std::cout << "----------------------------------------------" << std::endl;
775 std::cout << "## Test :: Vector" << std::endl;
776 std::cout << "----------------------------------------------" << std::endl;
777 std::cout << "----------------------------------------------" << std::endl;
778 std::cout << std::endl;
779
780 int retval = EXIT_SUCCESS;
781
782 std::cout << std::endl;
783 std::cout << "----------------------------------------------" << std::endl;
784 std::cout << std::endl;
785 {
786 typedef float NumericT;
787 NumericT epsilon = static_cast<NumericT>(1.0E-4);
788 std::cout << "# Testing setup:" << std::endl;
789 std::cout << " eps: " << epsilon << std::endl;
790 std::cout << " numeric: float" << std::endl;
791 retval = test<NumericT>(epsilon);
792 if ( retval == EXIT_SUCCESS )
793 std::cout << "# Test passed" << std::endl;
794 else
795 return retval;
796 }
797 std::cout << std::endl;
798 std::cout << "----------------------------------------------" << std::endl;
799 std::cout << std::endl;
800 #ifdef VIENNACL_WITH_OPENCL
801 if ( viennacl::ocl::current_device().double_support() )
802 #endif
803 {
804 {
805 typedef double NumericT;
806 NumericT epsilon = 1.0E-12;
807 std::cout << "# Testing setup:" << std::endl;
808 std::cout << " eps: " << epsilon << std::endl;
809 std::cout << " numeric: double" << std::endl;
810 retval = test<NumericT>(epsilon);
811 if ( retval == EXIT_SUCCESS )
812 std::cout << "# Test passed" << std::endl;
813 else
814 return retval;
815 }
816 std::cout << std::endl;
817 std::cout << "----------------------------------------------" << std::endl;
818 std::cout << std::endl;
819 }
820
821 std::cout << std::endl;
822 std::cout << "------- Test completed --------" << std::endl;
823 std::cout << std::endl;
824
825
826 return retval;
827 }
828