ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
scheduler_vector.cpp
Go to the documentation of this file.
1 /* =========================================================================
2  Copyright (c) 2010-2015, 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 
24 //
25 // *** System
26 //
27 #include <iostream>
28 #include <iomanip>
29 
30 //
31 // *** Boost
32 //
33 #include <boost/numeric/ublas/io.hpp>
34 #include <boost/numeric/ublas/vector.hpp>
35 #include <boost/numeric/ublas/vector_proxy.hpp>
36 
37 //
38 // *** ViennaCL
39 //
40 //#define VIENNACL_DEBUG_ALL
41 #define VIENNACL_WITH_UBLAS 1
42 #include "viennacl/vector.hpp"
43 #include "viennacl/matrix.hpp"
44 #include "viennacl/vector_proxy.hpp"
45 #include "viennacl/linalg/inner_prod.hpp"
46 #include "viennacl/linalg/norm_1.hpp"
47 #include "viennacl/linalg/norm_2.hpp"
48 #include "viennacl/linalg/norm_inf.hpp"
49 
50 #include "viennacl/scheduler/execute.hpp"
51 #include "viennacl/scheduler/io.hpp"
52 
53 #include "viennacl/tools/random.hpp"
54 
55 using namespace boost::numeric;
56 
57 
58 //
59 // -------------------------------------------------------------
60 //
61 template<typename ScalarType>
62 ScalarType diff(ScalarType const & s1, ScalarType const & s2)
63 {
64  viennacl::backend::finish();
65  if (std::fabs(s1 - s2) > 0)
66  return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
67  return 0;
68 }
69 //
70 // -------------------------------------------------------------
71 //
72 template<typename ScalarType>
73 ScalarType diff(ScalarType const & s1, viennacl::scalar<ScalarType> const & s2)
74 {
75  viennacl::backend::finish();
76  if (std::fabs(s1 - s2) > 0)
77  return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
78  return 0;
79 }
80 //
81 // -------------------------------------------------------------
82 //
83 template<typename ScalarType>
84 ScalarType diff(ScalarType const & s1, viennacl::entry_proxy<ScalarType> const & s2)
85 {
86  viennacl::backend::finish();
87  if (std::fabs(s1 - s2) > 0)
88  return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
89  return 0;
90 }
91 //
92 // -------------------------------------------------------------
93 //
94 template<typename ScalarType, typename ViennaCLVectorType>
95 ScalarType diff(ublas::vector<ScalarType> const & v1, ViennaCLVectorType const & vcl_vec)
96 {
97  ublas::vector<ScalarType> v2_cpu(vcl_vec.size());
98  viennacl::backend::finish();
99  viennacl::copy(vcl_vec, v2_cpu);
100 
101  for (unsigned int i=0;i<v1.size(); ++i)
102  {
103  if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
104  v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
105  else
106  v2_cpu[i] = 0.0;
107  }
108 
109  return ublas::norm_inf(v2_cpu);
110 }
111 
112 
113 template<typename T1, typename T2>
114 int check(T1 const & t1, T2 const & t2, double epsilon)
115 {
116  int retval = EXIT_SUCCESS;
117 
118  double temp = std::fabs(diff(t1, t2));
119  if (temp > epsilon)
120  {
121  std::cout << "# Error! Relative difference: " << temp << std::endl;
122  retval = EXIT_FAILURE;
123  }
124  else
125  std::cout << "PASSED!" << std::endl;
126  return retval;
127 }
128 
129 
130 //
131 // -------------------------------------------------------------
132 //
133 template< typename NumericT, typename Epsilon, typename UblasVectorType, typename ViennaCLVectorType1, typename ViennaCLVectorType2 >
134 int test(Epsilon const& epsilon,
135  UblasVectorType & ublas_v1, UblasVectorType & ublas_v2,
136  ViennaCLVectorType1 & vcl_v1, ViennaCLVectorType2 & vcl_v2)
137 {
138  int retval = EXIT_SUCCESS;
139 
140  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
141 
142  NumericT cpu_result = 42.0;
143  viennacl::scalar<NumericT> gpu_result = 43.0;
144  NumericT alpha = NumericT(3.1415);
145  NumericT beta = NumericT(2.7172);
146 
147  //
148  // Initializer:
149  //
150  std::cout << "Checking for zero_vector initializer..." << std::endl;
151  ublas_v1 = ublas::zero_vector<NumericT>(ublas_v1.size());
152  vcl_v1 = viennacl::zero_vector<NumericT>(vcl_v1.size());
153  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
154  return EXIT_FAILURE;
155 
156  std::cout << "Checking for scalar_vector initializer..." << std::endl;
157  ublas_v1 = ublas::scalar_vector<NumericT>(ublas_v1.size(), cpu_result);
158  vcl_v1 = viennacl::scalar_vector<NumericT>(vcl_v1.size(), cpu_result);
159  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
160  return EXIT_FAILURE;
161 
162  ublas_v1 = ublas::scalar_vector<NumericT>(ublas_v1.size(), gpu_result);
163  vcl_v1 = viennacl::scalar_vector<NumericT>(vcl_v1.size(), gpu_result);
164  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
165  return EXIT_FAILURE;
166 
167  std::cout << "Checking for unit_vector initializer..." << std::endl;
168  ublas_v1 = ublas::unit_vector<NumericT>(ublas_v1.size(), 5);
169  vcl_v1 = viennacl::unit_vector<NumericT>(vcl_v1.size(), 5);
170  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
171  return EXIT_FAILURE;
172 
173 
174  for (std::size_t i=0; i<ublas_v1.size(); ++i)
175  {
176  ublas_v1[i] = NumericT(1.0) + randomNumber();
177  ublas_v2[i] = NumericT(1.0) + randomNumber();
178  }
179 
180  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin()); //resync
181  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
182 
183  std::cout << "Checking for successful copy..." << std::endl;
184  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
185  return EXIT_FAILURE;
186  if (check(ublas_v2, vcl_v2, epsilon) != EXIT_SUCCESS)
187  return EXIT_FAILURE;
188 
189 
190  // --------------------------------------------------------------------------
191 
192  std::cout << "Testing simple assignments..." << std::endl;
193 
194  {
195  ublas_v1 = ublas_v2;
196  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), vcl_v2); // same as vcl_v1 = vcl_v2;
197  viennacl::scheduler::execute(my_statement);
198 
199  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
200  return EXIT_FAILURE;
201  }
202 
203  {
204  ublas_v1 += ublas_v2;
205  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_add(), vcl_v2); // same as vcl_v1 += vcl_v2;
206  viennacl::scheduler::execute(my_statement);
207 
208  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
209  return EXIT_FAILURE;
210  }
211 
212  {
213  ublas_v1 -= ublas_v2;
214  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_sub(), vcl_v2); // same as vcl_v1 -= vcl_v2;
215  viennacl::scheduler::execute(my_statement);
216 
217  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
218  return EXIT_FAILURE;
219  }
220 
221  std::cout << "Testing composite assignments..." << std::endl;
222  {
223  ublas_v1 = ublas_v1 + ublas_v2;
224  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), vcl_v1 + vcl_v2); // same as vcl_v1 = vcl_v1 + vcl_v2;
225  viennacl::scheduler::execute(my_statement);
226 
227  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
228  return EXIT_FAILURE;
229  }
230  {
231  ublas_v1 += alpha * ublas_v1 - beta * ublas_v2 + ublas_v1 / beta - ublas_v2 / alpha;
232  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;
233  viennacl::scheduler::execute(my_statement);
234 
235  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
236  return EXIT_FAILURE;
237  }
238 
239  {
240  ublas_v1 = ublas_v1 - ublas_v2;
241  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), vcl_v1 - vcl_v2); // same as vcl_v1 = vcl_v1 - vcl_v2;
242  viennacl::scheduler::execute(my_statement);
243 
244  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
245  return EXIT_FAILURE;
246  }
247 
248  std::cout << "--- Testing reductions ---" << std::endl;
249  std::cout << "inner_prod..." << std::endl;
250  {
251  cpu_result = inner_prod(ublas_v1, ublas_v2);
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 = inner_prod(ublas_v1 + ublas_v2, ublas_v2);
261  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);
262  viennacl::scheduler::execute(my_statement);
263 
264  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
265  return EXIT_FAILURE;
266  }
267 
268  {
269  cpu_result = inner_prod(ublas_v1, ublas_v2 - ublas_v1);
270  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);
271  viennacl::scheduler::execute(my_statement);
272 
273  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
274  return EXIT_FAILURE;
275  }
276 
277  {
278  cpu_result = inner_prod(ublas_v1 - ublas_v2, ublas_v2 + ublas_v1);
279  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);
280  viennacl::scheduler::execute(my_statement);
281 
282  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
283  return EXIT_FAILURE;
284  }
285 
286  std::cout << "norm_1..." << std::endl;
287  {
288  cpu_result = norm_1(ublas_v1);
289  viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_1(vcl_v1)); // same as gpu_result = norm_1(vcl_v1);
290  viennacl::scheduler::execute(my_statement);
291 
292  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
293  return EXIT_FAILURE;
294  }
295 
296  {
297  cpu_result = norm_1(ublas_v1 + ublas_v2);
298  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);
299  viennacl::scheduler::execute(my_statement);
300 
301  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
302  return EXIT_FAILURE;
303  }
304 
305  std::cout << "norm_2..." << std::endl;
306  {
307  cpu_result = norm_2(ublas_v1);
308  viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_2(vcl_v1)); // same as gpu_result = norm_2(vcl_v1);
309  viennacl::scheduler::execute(my_statement);
310 
311  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
312  return EXIT_FAILURE;
313  }
314 
315  {
316  cpu_result = norm_2(ublas_v1 + ublas_v2);
317  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);
318  viennacl::scheduler::execute(my_statement);
319 
320  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
321  return EXIT_FAILURE;
322  }
323 
324  std::cout << "norm_inf..." << std::endl;
325  {
326  cpu_result = norm_inf(ublas_v1);
327  viennacl::scheduler::statement my_statement(gpu_result, viennacl::op_assign(), viennacl::linalg::norm_inf(vcl_v1)); // same as gpu_result = norm_inf(vcl_v1);
328  viennacl::scheduler::execute(my_statement);
329 
330  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
331  return EXIT_FAILURE;
332  }
333 
334  {
335  cpu_result = norm_inf(ublas_v1 - ublas_v2);
336  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);
337  viennacl::scheduler::execute(my_statement);
338 
339  if (check(cpu_result, gpu_result, epsilon) != EXIT_SUCCESS)
340  return EXIT_FAILURE;
341  }
342 
343  std::cout << "--- Testing elementwise operations (binary) ---" << std::endl;
344  std::cout << "x = element_prod(x, y)... ";
345  {
346  ublas_v1 = element_prod(ublas_v1, ublas_v2);
347  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_v1, vcl_v2));
348  viennacl::scheduler::execute(my_statement);
349 
350  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
351  return EXIT_FAILURE;
352  }
353 
354  std::cout << "x = element_prod(x + y, y)... ";
355  {
356  ublas_v1 = element_prod(ublas_v1 + ublas_v2, ublas_v2);
357  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_v1 + vcl_v2, vcl_v2));
358  viennacl::scheduler::execute(my_statement);
359 
360  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
361  return EXIT_FAILURE;
362  }
363 
364  std::cout << "x = element_prod(x, x + y)... ";
365  {
366  ublas_v1 = element_prod(ublas_v1, ublas_v1 + ublas_v2);
367  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_v1, vcl_v2 + vcl_v1));
368  viennacl::scheduler::execute(my_statement);
369 
370  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
371  return EXIT_FAILURE;
372  }
373 
374  std::cout << "x = element_prod(x - y, y + x)... ";
375  {
376  ublas_v1 = element_prod(ublas_v1 - ublas_v2, ublas_v2 + ublas_v1);
377  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_v1 - vcl_v2, vcl_v2 + vcl_v1));
378  viennacl::scheduler::execute(my_statement);
379 
380  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
381  return EXIT_FAILURE;
382  }
383 
384 
385 
386  std::cout << "x = element_div(x, y)... ";
387  {
388  ublas_v1 = element_div(ublas_v1, ublas_v2);
389  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_div(vcl_v1, vcl_v2));
390  viennacl::scheduler::execute(my_statement);
391 
392  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
393  return EXIT_FAILURE;
394  }
395 
396  std::cout << "x = element_div(x + y, y)... ";
397  {
398  ublas_v1 = element_div(ublas_v1 + ublas_v2, ublas_v2);
399  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_div(vcl_v1 + vcl_v2, vcl_v2));
400  viennacl::scheduler::execute(my_statement);
401 
402  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
403  return EXIT_FAILURE;
404  }
405 
406  std::cout << "x = element_div(x, x + y)... ";
407  {
408  ublas_v1 = element_div(ublas_v1, ublas_v1 + ublas_v2);
409  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_div(vcl_v1, vcl_v2 + vcl_v1));
410  viennacl::scheduler::execute(my_statement);
411 
412  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
413  return EXIT_FAILURE;
414  }
415 
416  std::cout << "x = element_div(x - y, y + x)... ";
417  {
418  ublas_v1 = element_div(ublas_v1 - ublas_v2, ublas_v2 + ublas_v1);
419  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_div(vcl_v1 - vcl_v2, vcl_v2 + vcl_v1));
420  viennacl::scheduler::execute(my_statement);
421 
422  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
423  return EXIT_FAILURE;
424  }
425 
426 
427  std::cout << "x = element_pow(x, y)... ";
428  {
429  for (std::size_t i=0; i<ublas_v1.size(); ++i)
430  {
431  ublas_v1[i] = NumericT(2.0) + randomNumber();
432  ublas_v2[i] = NumericT(1.0) + randomNumber();
433  }
434  viennacl::copy(ublas_v1, vcl_v1);
435  viennacl::copy(ublas_v2, vcl_v2);
436 
437  for (std::size_t i=0; i<ublas_v1.size(); ++i)
438  ublas_v1[i] = std::pow(ublas_v1[i], ublas_v2[i]);
439  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_pow(vcl_v1, vcl_v2));
440  viennacl::scheduler::execute(my_statement);
441 
442  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
443  return EXIT_FAILURE;
444  }
445 
446  std::cout << "x = element_pow(x + y, y)... ";
447  {
448  for (std::size_t i=0; i<ublas_v1.size(); ++i)
449  {
450  ublas_v1[i] = NumericT(2.0) + randomNumber();
451  ublas_v2[i] = NumericT(1.0) + randomNumber();
452  }
453  viennacl::copy(ublas_v1, vcl_v1);
454  viennacl::copy(ublas_v2, vcl_v2);
455 
456  for (std::size_t i=0; i<ublas_v1.size(); ++i)
457  ublas_v1[i] = std::pow(ublas_v1[i] + ublas_v2[i], ublas_v2[i]);
458  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_pow(vcl_v1 + vcl_v2, vcl_v2));
459  viennacl::scheduler::execute(my_statement);
460 
461  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
462  return EXIT_FAILURE;
463  }
464 
465  std::cout << "x = element_pow(x, x + y)... ";
466  {
467  for (std::size_t i=0; i<ublas_v1.size(); ++i)
468  {
469  ublas_v1[i] = NumericT(2.0) + randomNumber();
470  ublas_v2[i] = NumericT(1.0) + randomNumber();
471  }
472  viennacl::copy(ublas_v1, vcl_v1);
473  viennacl::copy(ublas_v2, vcl_v2);
474 
475  for (std::size_t i=0; i<ublas_v1.size(); ++i)
476  ublas_v1[i] = std::pow(ublas_v1[i], ublas_v1[i] + ublas_v2[i]);
477  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_pow(vcl_v1, vcl_v2 + vcl_v1));
478  viennacl::scheduler::execute(my_statement);
479 
480  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
481  return EXIT_FAILURE;
482  }
483 
484  std::cout << "x = element_pow(x - y, y + x)... ";
485  {
486  for (std::size_t i=0; i<ublas_v1.size(); ++i)
487  {
488  ublas_v1[i] = NumericT(2.0) + randomNumber();
489  ublas_v2[i] = NumericT(1.0) + randomNumber();
490  }
491  viennacl::copy(ublas_v1, vcl_v1);
492  viennacl::copy(ublas_v2, vcl_v2);
493 
494  for (std::size_t i=0; i<ublas_v1.size(); ++i)
495  ublas_v1[i] = std::pow(ublas_v1[i] - ublas_v2[i], ublas_v2[i] + ublas_v1[i]);
496  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_pow(vcl_v1 - vcl_v2, vcl_v2 + vcl_v1));
497  viennacl::scheduler::execute(my_statement);
498 
499  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
500  return EXIT_FAILURE;
501  }
502 
503  std::cout << "--- Testing elementwise operations (unary) ---" << std::endl;
504 #define GENERATE_UNARY_OP_TEST(OPNAME) \
505  ublas_v1 = ublas::scalar_vector<NumericT>(ublas_v1.size(), NumericT(0.21)); \
506  ublas_v2 = NumericT(3.1415) * ublas_v1; \
507  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin()); \
508  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin()); \
509  { \
510  for (std::size_t i=0; i<ublas_v1.size(); ++i) \
511  ublas_v1[i] = std::OPNAME(ublas_v2[i]); \
512  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_##OPNAME(vcl_v2)); \
513  viennacl::scheduler::execute(my_statement); \
514  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS) \
515  return EXIT_FAILURE; \
516  } \
517  { \
518  for (std::size_t i=0; i<ublas_v1.size(); ++i) \
519  ublas_v1[i] = std::OPNAME(ublas_v2[i] / NumericT(2)); \
520  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::element_##OPNAME(vcl_v2 / NumericT(2))); \
521  viennacl::scheduler::execute(my_statement); \
522  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS) \
523  return EXIT_FAILURE; \
524  }
525 
526  GENERATE_UNARY_OP_TEST(cos);
527  GENERATE_UNARY_OP_TEST(cosh);
528  GENERATE_UNARY_OP_TEST(exp);
529  GENERATE_UNARY_OP_TEST(floor);
530  GENERATE_UNARY_OP_TEST(fabs);
531  GENERATE_UNARY_OP_TEST(log);
532  GENERATE_UNARY_OP_TEST(log10);
533  GENERATE_UNARY_OP_TEST(sin);
534  GENERATE_UNARY_OP_TEST(sinh);
535  GENERATE_UNARY_OP_TEST(fabs);
536  //GENERATE_UNARY_OP_TEST(abs); //OpenCL allows abs on integers only
537  GENERATE_UNARY_OP_TEST(sqrt);
538  GENERATE_UNARY_OP_TEST(tan);
539  GENERATE_UNARY_OP_TEST(tanh);
540 
541 #undef GENERATE_UNARY_OP_TEST
542 
543  std::cout << "--- Testing complicated composite operations ---" << std::endl;
544  std::cout << "x = inner_prod(x, y) * y..." << std::endl;
545  {
546  ublas_v1 = inner_prod(ublas_v1, ublas_v2) * ublas_v2;
547  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::inner_prod(vcl_v1, vcl_v2) * vcl_v2);
548  viennacl::scheduler::execute(my_statement);
549 
550  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
551  return EXIT_FAILURE;
552  }
553 
554  std::cout << "x = y / norm_1(x)..." << std::endl;
555  {
556  ublas_v1 = ublas_v2 / norm_1(ublas_v1);
557  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), vcl_v2 / viennacl::linalg::norm_1(vcl_v1) );
558  viennacl::scheduler::execute(my_statement);
559 
560  if (check(ublas_v1, vcl_v1, epsilon) != EXIT_SUCCESS)
561  return EXIT_FAILURE;
562  }
563 
564 
565  // --------------------------------------------------------------------------
566  return retval;
567 }
568 
569 
570 template< typename NumericT, typename Epsilon >
571 int test(Epsilon const& epsilon)
572 {
573  int retval = EXIT_SUCCESS;
574  std::size_t size = 24656;
575 
576  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
577 
578  std::cout << "Running tests for vector of size " << size << std::endl;
579 
580  //
581  // Set up UBLAS objects
582  //
583  ublas::vector<NumericT> ublas_full_vec(size);
584  ublas::vector<NumericT> ublas_full_vec2(ublas_full_vec.size());
585 
586  for (std::size_t i=0; i<ublas_full_vec.size(); ++i)
587  {
588  ublas_full_vec[i] = NumericT(1.0) + randomNumber();
589  ublas_full_vec2[i] = NumericT(1.0) + randomNumber();
590  }
591 
592  ublas::range r1( ublas_full_vec.size() / 4, 2 * ublas_full_vec.size() / 4);
593  ublas::range r2(2 * ublas_full_vec2.size() / 4, 3 * ublas_full_vec2.size() / 4);
594  ublas::vector_range< ublas::vector<NumericT> > ublas_range_vec(ublas_full_vec, r1);
595  ublas::vector_range< ublas::vector<NumericT> > ublas_range_vec2(ublas_full_vec2, r2);
596 
597  ublas::slice s1( ublas_full_vec.size() / 4, 3, ublas_full_vec.size() / 4);
598  ublas::slice s2(2 * ublas_full_vec2.size() / 4, 2, ublas_full_vec2.size() / 4);
599  ublas::vector_slice< ublas::vector<NumericT> > ublas_slice_vec(ublas_full_vec, s1);
600  ublas::vector_slice< ublas::vector<NumericT> > ublas_slice_vec2(ublas_full_vec2, s2);
601 
602  //
603  // Set up ViennaCL objects
604  //
605  viennacl::vector<NumericT> vcl_full_vec(ublas_full_vec.size());
606  viennacl::vector<NumericT> vcl_full_vec2(ublas_full_vec2.size());
607 
608  viennacl::fast_copy(ublas_full_vec.begin(), ublas_full_vec.end(), vcl_full_vec.begin());
609  viennacl::copy(ublas_full_vec2.begin(), ublas_full_vec2.end(), vcl_full_vec2.begin());
610 
611  viennacl::range vcl_r1( vcl_full_vec.size() / 4, 2 * vcl_full_vec.size() / 4);
612  viennacl::range vcl_r2(2 * vcl_full_vec2.size() / 4, 3 * vcl_full_vec2.size() / 4);
613  viennacl::vector_range< viennacl::vector<NumericT> > vcl_range_vec(vcl_full_vec, vcl_r1);
614  viennacl::vector_range< viennacl::vector<NumericT> > vcl_range_vec2(vcl_full_vec2, vcl_r2);
615 
616  {
617  viennacl::vector<NumericT> vcl_short_vec(vcl_range_vec);
618  viennacl::vector<NumericT> vcl_short_vec2 = vcl_range_vec2;
619 
620  ublas::vector<NumericT> ublas_short_vec(ublas_range_vec);
621  ublas::vector<NumericT> ublas_short_vec2(ublas_range_vec2);
622 
623  std::cout << "Testing creation of vectors from range..." << std::endl;
624  if (check(ublas_short_vec, vcl_short_vec, epsilon) != EXIT_SUCCESS)
625  return EXIT_FAILURE;
626  if (check(ublas_short_vec2, vcl_short_vec2, epsilon) != EXIT_SUCCESS)
627  return EXIT_FAILURE;
628  }
629 
630  viennacl::slice vcl_s1( vcl_full_vec.size() / 4, 3, vcl_full_vec.size() / 4);
631  viennacl::slice vcl_s2(2 * vcl_full_vec2.size() / 4, 2, vcl_full_vec2.size() / 4);
632  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_slice_vec(vcl_full_vec, vcl_s1);
633  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_slice_vec2(vcl_full_vec2, vcl_s2);
634 
635  viennacl::vector<NumericT> vcl_short_vec(vcl_slice_vec);
636  viennacl::vector<NumericT> vcl_short_vec2 = vcl_slice_vec2;
637 
638  ublas::vector<NumericT> ublas_short_vec(ublas_slice_vec);
639  ublas::vector<NumericT> ublas_short_vec2(ublas_slice_vec2);
640 
641  std::cout << "Testing creation of vectors from slice..." << std::endl;
642  if (check(ublas_short_vec, vcl_short_vec, epsilon) != EXIT_SUCCESS)
643  return EXIT_FAILURE;
644  if (check(ublas_short_vec2, vcl_short_vec2, epsilon) != EXIT_SUCCESS)
645  return EXIT_FAILURE;
646 
647 
648  //
649  // Now start running tests for vectors, ranges and slices:
650  //
651 
652  std::cout << " ** vcl_v1 = vector, vcl_v2 = vector **" << std::endl;
653  retval = test<NumericT>(epsilon,
654  ublas_short_vec, ublas_short_vec2,
655  vcl_short_vec, vcl_short_vec2);
656  if (retval != EXIT_SUCCESS)
657  return EXIT_FAILURE;
658 
659  std::cout << " ** vcl_v1 = vector, vcl_v2 = range **" << std::endl;
660  retval = test<NumericT>(epsilon,
661  ublas_short_vec, ublas_short_vec2,
662  vcl_short_vec, vcl_range_vec2);
663  if (retval != EXIT_SUCCESS)
664  return EXIT_FAILURE;
665 
666  std::cout << " ** vcl_v1 = vector, vcl_v2 = slice **" << std::endl;
667  retval = test<NumericT>(epsilon,
668  ublas_short_vec, ublas_short_vec2,
669  vcl_short_vec, vcl_slice_vec2);
670  if (retval != EXIT_SUCCESS)
671  return EXIT_FAILURE;
672 
674 
675  std::cout << " ** vcl_v1 = range, vcl_v2 = vector **" << std::endl;
676  retval = test<NumericT>(epsilon,
677  ublas_short_vec, ublas_short_vec2,
678  vcl_range_vec, vcl_short_vec2);
679  if (retval != EXIT_SUCCESS)
680  return EXIT_FAILURE;
681 
682  std::cout << " ** vcl_v1 = range, vcl_v2 = range **" << std::endl;
683  retval = test<NumericT>(epsilon,
684  ublas_short_vec, ublas_short_vec2,
685  vcl_range_vec, vcl_range_vec2);
686  if (retval != EXIT_SUCCESS)
687  return EXIT_FAILURE;
688 
689  std::cout << " ** vcl_v1 = range, vcl_v2 = slice **" << std::endl;
690  retval = test<NumericT>(epsilon,
691  ublas_short_vec, ublas_short_vec2,
692  vcl_range_vec, vcl_slice_vec2);
693  if (retval != EXIT_SUCCESS)
694  return EXIT_FAILURE;
695 
697 
698  std::cout << " ** vcl_v1 = slice, vcl_v2 = vector **" << std::endl;
699  retval = test<NumericT>(epsilon,
700  ublas_short_vec, ublas_short_vec2,
701  vcl_slice_vec, vcl_short_vec2);
702  if (retval != EXIT_SUCCESS)
703  return EXIT_FAILURE;
704 
705  std::cout << " ** vcl_v1 = slice, vcl_v2 = range **" << std::endl;
706  retval = test<NumericT>(epsilon,
707  ublas_short_vec, ublas_short_vec2,
708  vcl_slice_vec, vcl_range_vec2);
709  if (retval != EXIT_SUCCESS)
710  return EXIT_FAILURE;
711 
712  std::cout << " ** vcl_v1 = slice, vcl_v2 = slice **" << std::endl;
713  retval = test<NumericT>(epsilon,
714  ublas_short_vec, ublas_short_vec2,
715  vcl_slice_vec, vcl_slice_vec2);
716  if (retval != EXIT_SUCCESS)
717  return EXIT_FAILURE;
718 
719  return EXIT_SUCCESS;
720 }
721 
722 
723 
724 //
725 // -------------------------------------------------------------
726 //
727 int main()
728 {
729  std::cout << std::endl;
730  std::cout << "----------------------------------------------" << std::endl;
731  std::cout << "----------------------------------------------" << std::endl;
732  std::cout << "## Test :: Vector" << std::endl;
733  std::cout << "----------------------------------------------" << std::endl;
734  std::cout << "----------------------------------------------" << std::endl;
735  std::cout << std::endl;
736 
737  int retval = EXIT_SUCCESS;
738 
739  std::cout << std::endl;
740  std::cout << "----------------------------------------------" << std::endl;
741  std::cout << std::endl;
742  {
743  typedef float NumericT;
744  NumericT epsilon = static_cast<NumericT>(1.0E-4);
745  std::cout << "# Testing setup:" << std::endl;
746  std::cout << " eps: " << epsilon << std::endl;
747  std::cout << " numeric: float" << std::endl;
748  retval = test<NumericT>(epsilon);
749  if ( retval == EXIT_SUCCESS )
750  std::cout << "# Test passed" << std::endl;
751  else
752  return retval;
753  }
754  std::cout << std::endl;
755  std::cout << "----------------------------------------------" << std::endl;
756  std::cout << std::endl;
757 #ifdef VIENNACL_WITH_OPENCL
758  if ( viennacl::ocl::current_device().double_support() )
759 #endif
760  {
761  {
762  typedef double NumericT;
763  NumericT epsilon = 1.0E-12;
764  std::cout << "# Testing setup:" << std::endl;
765  std::cout << " eps: " << epsilon << std::endl;
766  std::cout << " numeric: double" << std::endl;
767  retval = test<NumericT>(epsilon);
768  if ( retval == EXIT_SUCCESS )
769  std::cout << "# Test passed" << std::endl;
770  else
771  return retval;
772  }
773  std::cout << std::endl;
774  std::cout << "----------------------------------------------" << std::endl;
775  std::cout << std::endl;
776  }
777 
778  std::cout << std::endl;
779  std::cout << "------- Test completed --------" << std::endl;
780  std::cout << std::endl;
781 
782 
783  return retval;
784 }
viennacl::linalg::element_div
viennacl::vector_expression< const vector_base< T >, const vector_base< T >, op_element_binary< op_div > > element_div(vector_base< T > const &v1, vector_base< T > const &v2)
viennacl::linalg::norm_2
T norm_2(std::vector< T, A > const &v1)
Definition: norm_2.hpp:96
viennacl::scalar
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
norm_2.hpp
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
test
int test(Epsilon const &epsilon, UblasVectorType &ublas_v1, UblasVectorType &ublas_v2, ViennaCLVectorType1 &vcl_v1, ViennaCLVectorType2 &vcl_v2)
Definition: scheduler_vector.cpp:134
matrix.hpp
Implementation of the dense matrix class.
io.hpp
Some helper routines for reading/writing/printing scheduler expressions.
viennacl::op_assign
A tag class representing assignment.
Definition: forwards.h:81
viennacl::backend::finish
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
check
int check(T1 const &t1, T2 const &t2, double epsilon)
Definition: scheduler_vector.cpp:114
viennacl::zero_vector
Definition: vector_def.hpp:93
viennacl::linalg::inner_prod
viennacl::enable_if< viennacl::is_stl< typename viennacl::traits::tag_of< VectorT1 >::type >::value, typename VectorT1::value_type >::type inner_prod(VectorT1 const &v1, VectorT2 const &v2)
Definition: inner_prod.hpp:100
viennacl::scheduler::execute
void execute(statement const &s)
Definition: execute.hpp:279
s2
viennacl::scalar< int > s2
Definition: global_variables.cpp:58
s1
viennacl::scalar< float > s1
Definition: global_variables.cpp:57
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
viennacl::ocl::current_device
viennacl::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Definition: backend.hpp:351
inner_prod.hpp
Generic interface for the computation of inner products. See viennacl/linalg/vector_operations.hpp for implementations.
main
int main()
Definition: scheduler_vector.cpp:727
viennacl::op_inplace_add
A tag class representing inplace addition.
Definition: forwards.h:83
norm_1.hpp
Generic interface for the l^1-norm. See viennacl/linalg/vector_operations.hpp for implementations...
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::range
basic_range range
Definition: forwards.h:424
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
viennacl::traits::size
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:235
viennacl::tools::uniform_random_numbers
Random number generator for returning uniformly distributed values in the closed interval [0...
Definition: random.hpp:44
viennacl::vector_range
Class for representing non-strided subvectors of a bigger vector x.
Definition: forwards.h:434
viennacl::vector_slice
Class for representing strided subvectors of a bigger vector x.
Definition: forwards.h:437
viennacl::ocl::device::double_support
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Definition: device.hpp:956
GENERATE_UNARY_OP_TEST
#define GENERATE_UNARY_OP_TEST(OPNAME)
vector_proxy.hpp
Proxy classes for vectors.
viennacl::op_inplace_sub
A tag class representing inplace subtraction.
Definition: forwards.h:85
viennacl::unit_vector
Represents a vector consisting of 1 at a given index and zeros otherwise.
Definition: vector_def.hpp:76
viennacl::slice
basic_slice slice
Definition: forwards.h:429
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::scalar_vector
Represents a vector consisting of scalars 's' only, i.e. v[i] = s for all i. To be used as an initial...
Definition: vector_def.hpp:87
viennacl::linalg::norm_inf
T norm_inf(std::vector< T, A > const &v1)
Definition: norm_inf.hpp:60
viennacl::copy
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
random.hpp
A small collection of sequential random number generators.
viennacl::linalg::norm_1
T norm_1(std::vector< T, A > const &v1)
Definition: norm_1.hpp:61
viennacl::basic_range
A range class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:424
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
diff
ScalarType diff(ScalarType const &s1, ScalarType const &s2)
Definition: scheduler_vector.cpp:62
viennacl::linalg::element_prod
viennacl::vector_expression< const vector_base< T >, const vector_base< T >, op_element_binary< op_prod > > element_prod(vector_base< T > const &v1, vector_base< T > const &v2)
execute.hpp
Provides the datastructures for dealing with a single statement such as 'x = y + z;'.
viennacl::scheduler::statement
The main class for representing a statement such as x = inner_prod(y,z); at runtime.
Definition: forwards.h:502
viennacl::basic_slice
A slice class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:429
viennacl::entry_proxy
A proxy class for a single element of a vector or matrix. This proxy should not be noticed by end-use...
Definition: forwards.h:233
norm_inf.hpp
Generic interface for the l^infty-norm. See viennacl/linalg/vector_operations.hpp for implementations...
viennacl::fast_copy
void fast_copy(const const_vector_iterator< SCALARTYPE, ALIGNMENT > &gpu_begin, const const_vector_iterator< SCALARTYPE, ALIGNMENT > &gpu_end, CPU_ITERATOR cpu_begin)