ViennaCL - The Vienna Computing Library  1.7.1
Free open-source GPU-accelerated linear algebra and solver library.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
scheduler_vector.cpp
Go to the documentation of this file.
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 
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>
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>
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>
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>
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 >
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 >
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 
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 
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 //
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 }
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...
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:106
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:770
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
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
test
int test(Epsilon const &epsilon, STLVectorType &std_v1, STLVectorType &std_v2, ViennaCLVectorType1 &vcl_v1, ViennaCLVectorType2 &vcl_v2)
Definition: scheduler_vector.cpp:126
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:239
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
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:52
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)