ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
matrix_float_double.hpp
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 #define VIENNACL_WITH_UBLAS
19 //#define NDEBUG
20 //#define VIENNACL_BUILD_INFO
21 
22 // We don't need debug mode in UBLAS:
23 #define BOOST_UBLAS_NDEBUG
24 
25 #include <utility>
26 #include <iostream>
27 #include <fstream>
28 #include <string>
29 #include <cmath>
30 #include <algorithm>
31 #include <cstdio>
32 #include <ctime>
33 
34 #include "viennacl/scalar.hpp"
35 #include "viennacl/matrix.hpp"
36 #include "viennacl/linalg/prod.hpp"
37 #include "viennacl/linalg/norm_1.hpp"
38 #include "viennacl/linalg/norm_inf.hpp"
39 #include "viennacl/linalg/norm_frobenius.hpp"
40 #include "viennacl/matrix_proxy.hpp"
41 #include "viennacl/vector_proxy.hpp"
42 #include "viennacl/linalg/norm_1.hpp"
43 #include "viennacl/linalg/norm_2.hpp"
44 #include "viennacl/linalg/norm_inf.hpp"
45 #include "viennacl/linalg/norm_frobenius.hpp"
46 
47 #include "boost/numeric/ublas/vector.hpp"
48 #include "boost/numeric/ublas/matrix.hpp"
49 #include "boost/numeric/ublas/matrix_proxy.hpp"
50 #include "boost/numeric/ublas/vector_proxy.hpp"
51 #include "boost/numeric/ublas/io.hpp"
52 
53 template<typename MatrixType, typename VCLMatrixType>
54 bool check_for_equality(MatrixType const & ublas_A, VCLMatrixType const & vcl_A, double epsilon)
55 {
56  typedef typename MatrixType::value_type value_type;
57 
58  boost::numeric::ublas::matrix<value_type> vcl_A_cpu(vcl_A.size1(), vcl_A.size2());
59  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
60  viennacl::copy(vcl_A, vcl_A_cpu);
61 
62  for (std::size_t i=0; i<ublas_A.size1(); ++i)
63  {
64  for (std::size_t j=0; j<ublas_A.size2(); ++j)
65  {
66  if (std::fabs(ublas_A(i,j) - vcl_A_cpu(i,j)) > 0)
67  {
68  if ( (std::abs(ublas_A(i,j) - vcl_A_cpu(i,j)) / std::max(std::fabs(ublas_A(i,j)), std::fabs(vcl_A_cpu(i,j))) > epsilon) || std::fabs(vcl_A_cpu(i,j) - vcl_A_cpu(i,j)) > 0 )
69  {
70  std::cout << "Error at index (" << i << ", " << j << "): " << ublas_A(i,j) << " vs " << vcl_A_cpu(i,j) << std::endl;
71  std::cout << std::endl << "TEST failed!" << std::endl;
72  return false;
73  }
74  }
75  }
76  }
77 
78  std::cout << "PASSED!" << std::endl;
79  return true;
80 }
81 
82 
83 
84 
85 template<typename UBLASMatrixType,
86  typename ViennaCLMatrixType1, typename ViennaCLMatrixType2, typename ViennaCLMatrixType3>
87 int run_test(double epsilon,
88  UBLASMatrixType & ublas_A, UBLASMatrixType & ublas_B, UBLASMatrixType & ublas_C,
89  ViennaCLMatrixType1 & vcl_A, ViennaCLMatrixType2 & vcl_B, ViennaCLMatrixType3 vcl_C)
90 {
91 
92  typedef typename viennacl::result_of::cpu_value_type<typename ViennaCLMatrixType1::value_type>::type cpu_value_type;
93 
94  cpu_value_type alpha = cpu_value_type(3.1415);
95  viennacl::scalar<cpu_value_type> gpu_alpha = alpha;
96 
97  cpu_value_type beta = cpu_value_type(2.7182);
98  viennacl::scalar<cpu_value_type> gpu_beta = beta;
99 
100 
101  //
102  // Initializer:
103  //
104  std::cout << "Checking for zero_matrix initializer..." << std::endl;
105  ublas_A = boost::numeric::ublas::zero_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2());
106  vcl_A = viennacl::zero_matrix<cpu_value_type>(vcl_A.size1(), vcl_A.size2());
107  if (!check_for_equality(ublas_A, vcl_A, epsilon))
108  return EXIT_FAILURE;
109 
110  std::cout << "Checking for scalar_matrix initializer..." << std::endl;
111  ublas_A = boost::numeric::ublas::scalar_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2(), alpha);
112  vcl_A = viennacl::scalar_matrix<cpu_value_type>(vcl_A.size1(), vcl_A.size2(), alpha);
113  if (!check_for_equality(ublas_A, vcl_A, epsilon))
114  return EXIT_FAILURE;
115 
116  ublas_A = boost::numeric::ublas::scalar_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2(), gpu_beta);
117  vcl_A = viennacl::scalar_matrix<cpu_value_type>( vcl_A.size1(), vcl_A.size2(), gpu_beta);
118  if (!check_for_equality(ublas_A, vcl_A, epsilon))
119  return EXIT_FAILURE;
120 
121  /*std::cout << "Checking for identity initializer..." << std::endl;
122  ublas_A = boost::numeric::ublas::identity_matrix<cpu_value_type>(ublas_A.size1());
123  vcl_A = viennacl::identity_matrix<cpu_value_type>(vcl_A.size1());
124  if (!check_for_equality(ublas_A, vcl_A, epsilon))
125  return EXIT_FAILURE;*/
126 
127 
128  std::cout << std::endl;
129  //std::cout << "//" << std::endl;
130  //std::cout << "////////// Test: Assignments //////////" << std::endl;
131  //std::cout << "//" << std::endl;
132 
133  if (!check_for_equality(ublas_B, vcl_B, epsilon))
134  return EXIT_FAILURE;
135 
136  std::cout << "Testing matrix assignment... ";
137  //std::cout << ublas_B(0,0) << " vs. " << vcl_B(0,0) << std::endl;
138  ublas_A = ublas_B;
139  vcl_A = vcl_B;
140  if (!check_for_equality(ublas_A, vcl_A, epsilon))
141  return EXIT_FAILURE;
142 
143  if (ublas_A.size1() == ublas_A.size2())
144  {
145  std::cout << "Testing matrix assignment (transposed)... ";
146  ublas_A = trans(ublas_B);
147  vcl_A = viennacl::trans(vcl_B);
148 
149  if (!check_for_equality(ublas_A, vcl_A, epsilon))
150  return EXIT_FAILURE;
151  }
152 
153 
154 
155  //std::cout << std::endl;
156  //std::cout << "//" << std::endl;
157  //std::cout << "////////// Test 1: Copy to GPU //////////" << std::endl;
158  //std::cout << "//" << std::endl;
159 
160  ublas_A = ublas_B;
161  viennacl::copy(ublas_B, vcl_A);
162  std::cout << "Testing upper left copy to GPU... ";
163  if (!check_for_equality(ublas_A, vcl_A, epsilon))
164  return EXIT_FAILURE;
165 
166 
167  ublas_C = ublas_B;
168  viennacl::copy(ublas_B, vcl_C);
169  std::cout << "Testing lower right copy to GPU... ";
170  if (!check_for_equality(ublas_C, vcl_C, epsilon))
171  return EXIT_FAILURE;
172 
173 
174  //std::cout << std::endl;
175  //std::cout << "//" << std::endl;
176  //std::cout << "////////// Test 2: Copy from GPU //////////" << std::endl;
177  //std::cout << "//" << std::endl;
178 
179  std::cout << "Testing upper left copy to A... ";
180  if (!check_for_equality(ublas_A, vcl_A, epsilon))
181  return EXIT_FAILURE;
182 
183  std::cout << "Testing lower right copy to C... ";
184  if (!check_for_equality(ublas_C, vcl_C, epsilon))
185  return EXIT_FAILURE;
186 
187 
188 
189  //std::cout << "//" << std::endl;
190  //std::cout << "////////// Test 3: Addition //////////" << std::endl;
191  //std::cout << "//" << std::endl;
192  viennacl::copy(ublas_C, vcl_C);
193 
194  std::cout << "Inplace add: ";
195  ublas_C += ublas_C;
196  vcl_C += vcl_C;
197 
198  if (!check_for_equality(ublas_C, vcl_C, epsilon))
199  return EXIT_FAILURE;
200 
201  if (ublas_C.size1() == ublas_C.size2())
202  {
203  std::cout << "Inplace add (transposed): ";
204  ublas_C += trans(ublas_C);
205  vcl_C += viennacl::trans(vcl_C);
206 
207  if (!check_for_equality(ublas_C, vcl_C, epsilon))
208  return EXIT_FAILURE;
209  }
210 
211  std::cout << "Scaled inplace add: ";
212  ublas_C += beta * ublas_A;
213  vcl_C += gpu_beta * vcl_A;
214 
215  if (!check_for_equality(ublas_C, vcl_C, epsilon))
216  return EXIT_FAILURE;
217 
218  std::cout << "Add: ";
219  ublas_C = ublas_A + ublas_B;
220  vcl_C = vcl_A + vcl_B;
221 
222  if (!check_for_equality(ublas_C, vcl_C, epsilon))
223  return EXIT_FAILURE;
224 
225  std::cout << "Add with flipsign: ";
226  ublas_C = - ublas_A + ublas_B;
227  vcl_C = - vcl_A + vcl_B;
228 
229  if (!check_for_equality(ublas_C, vcl_C, epsilon))
230  return EXIT_FAILURE;
231 
232 
233  std::cout << "Scaled add (left): ";
234  ublas_C = long(alpha) * ublas_A + ublas_B;
235  vcl_C = long(alpha) * vcl_A + vcl_B;
236 
237  if (!check_for_equality(ublas_C, vcl_C, epsilon))
238  return EXIT_FAILURE;
239 
240  ublas_C = float(alpha) * ublas_A + ublas_B;
241  vcl_C = float(alpha) * vcl_A + vcl_B;
242 
243  if (!check_for_equality(ublas_C, vcl_C, epsilon))
244  return EXIT_FAILURE;
245 
246  ublas_C = double(alpha) * ublas_A + ublas_B;
247  vcl_C = double(alpha) * vcl_A + vcl_B;
248 
249  if (!check_for_equality(ublas_C, vcl_C, epsilon))
250  return EXIT_FAILURE;
251 
252  std::cout << "Scaled add (left): ";
253  vcl_C = gpu_alpha * vcl_A + vcl_B;
254  if (!check_for_equality(ublas_C, vcl_C, epsilon))
255  return EXIT_FAILURE;
256 
257 
258  std::cout << "Scaled add (right): ";
259  ublas_C = ublas_A + ublas_B * long(beta);
260  vcl_C = vcl_A + vcl_B * long(beta);
261 
262  if (!check_for_equality(ublas_C, vcl_C, epsilon))
263  return EXIT_FAILURE;
264 
265  ublas_C = ublas_A + ublas_B * float(beta);
266  vcl_C = vcl_A + vcl_B * float(beta);
267 
268  if (!check_for_equality(ublas_C, vcl_C, epsilon))
269  return EXIT_FAILURE;
270 
271  ublas_C = ublas_A + ublas_B * double(beta);
272  vcl_C = vcl_A + vcl_B * double(beta);
273 
274  if (!check_for_equality(ublas_C, vcl_C, epsilon))
275  return EXIT_FAILURE;
276 
277  std::cout << "Scaled add (right): ";
278  vcl_C = vcl_A + gpu_beta * vcl_B;
279  if (!check_for_equality(ublas_C, vcl_C, epsilon))
280  return EXIT_FAILURE;
281 
282  std::cout << "Scaled add (right, with division): ";
283  ublas_C = ublas_A + ublas_B / long(beta);
284  vcl_C = vcl_A + vcl_B / long(beta);
285 
286  if (!check_for_equality(ublas_C, vcl_C, epsilon))
287  return EXIT_FAILURE;
288 
289  ublas_C = ublas_A + ublas_B / float(beta);
290  vcl_C = vcl_A + vcl_B / float(beta);
291 
292  if (!check_for_equality(ublas_C, vcl_C, epsilon))
293  return EXIT_FAILURE;
294 
295  ublas_C = ublas_A + ublas_B / double(beta);
296  vcl_C = vcl_A + vcl_B / double(beta);
297 
298  if (!check_for_equality(ublas_C, vcl_C, epsilon))
299  return EXIT_FAILURE;
300 
301 
302  std::cout << "Scaled add (both): ";
303  ublas_C = alpha * ublas_A + beta * ublas_B;
304  vcl_C = alpha * vcl_A + beta * vcl_B;
305 
306  if (!check_for_equality(ublas_C, vcl_C, epsilon))
307  return EXIT_FAILURE;
308 
309  std::cout << "Scaled add (both): ";
310  vcl_C = gpu_alpha * vcl_A + gpu_beta * vcl_B;
311  if (!check_for_equality(ublas_C, vcl_C, epsilon))
312  return EXIT_FAILURE;
313 
314  //std::cout << "//" << std::endl;
315  //std::cout << "////////// Test 4: Subtraction //////////" << std::endl;
316  //std::cout << "//" << std::endl;
317  viennacl::copy(ublas_C, vcl_C);
318 
319  std::cout << "Inplace sub: ";
320  ublas_C -= ublas_B;
321  vcl_C -= vcl_B;
322 
323  if (!check_for_equality(ublas_C, vcl_C, epsilon))
324  return EXIT_FAILURE;
325 
326  if (ublas_C.size1() == ublas_C.size2())
327  {
328  std::cout << "Inplace add (transposed): ";
329  ublas_C -= trans(ublas_C);
330  vcl_C -= viennacl::trans(vcl_C);
331  }
332 
333  if (!check_for_equality(ublas_C, vcl_C, epsilon))
334  return EXIT_FAILURE;
335 
336  std::cout << "Scaled Inplace sub: ";
337  ublas_C -= alpha * ublas_B;
338  vcl_C -= alpha * vcl_B;
339 
340  if (!check_for_equality(ublas_C, vcl_C, epsilon))
341  return EXIT_FAILURE;
342 
343 
344 
345 
346  std::cout << "Sub: ";
347  ublas_C = ublas_A - ublas_B;
348  vcl_C = vcl_A - vcl_B;
349 
350  if (!check_for_equality(ublas_C, vcl_C, epsilon))
351  return EXIT_FAILURE;
352 
353  std::cout << "Scaled sub (left): ";
354  ublas_B = alpha * ublas_A - ublas_C;
355  vcl_B = alpha * vcl_A - vcl_C;
356 
357  if (!check_for_equality(ublas_B, vcl_B, epsilon))
358  return EXIT_FAILURE;
359 
360  std::cout << "Scaled sub (left): ";
361  vcl_B = gpu_alpha * vcl_A - vcl_C;
362  if (!check_for_equality(ublas_B, vcl_B, epsilon))
363  return EXIT_FAILURE;
364 
365 
366  std::cout << "Scaled sub (right): ";
367  ublas_B = ublas_A - beta * ublas_C;
368  vcl_B = vcl_A - vcl_C * beta;
369 
370  if (!check_for_equality(ublas_B, vcl_B, epsilon))
371  return EXIT_FAILURE;
372 
373  std::cout << "Scaled sub (right): ";
374  vcl_B = vcl_A - vcl_C * gpu_beta;
375  if (!check_for_equality(ublas_B, vcl_B, epsilon))
376  return EXIT_FAILURE;
377 
378 
379  std::cout << "Scaled sub (both): ";
380  ublas_B = alpha * ublas_A - beta * ublas_C;
381  vcl_B = alpha * vcl_A - vcl_C * beta;
382 
383  if (!check_for_equality(ublas_B, vcl_B, epsilon))
384  return EXIT_FAILURE;
385 
386  std::cout << "Scaled sub (both): ";
387  vcl_B = gpu_alpha * vcl_A - vcl_C * gpu_beta;
388  if (!check_for_equality(ublas_B, vcl_B, epsilon))
389  return EXIT_FAILURE;
390 
391 
392  std::cout << "Unary operator-: ";
393  ublas_C = - ublas_A;
394  vcl_C = - vcl_A;
395 
396  if (!check_for_equality(ublas_C, vcl_C, epsilon))
397  return EXIT_FAILURE;
398 
399 
400 
401  //std::cout << "//" << std::endl;
402  //std::cout << "////////// Test 5: Scaling //////////" << std::endl;
403  //std::cout << "//" << std::endl;
404  viennacl::copy(ublas_A, vcl_A);
405 
406  std::cout << "Multiplication with CPU scalar: ";
407  ublas_A *= long(alpha);
408  vcl_A *= long(alpha);
409 
410  if (!check_for_equality(ublas_A, vcl_A, epsilon))
411  return EXIT_FAILURE;
412 
413  ublas_A *= float(alpha);
414  vcl_A *= float(alpha);
415 
416  if (!check_for_equality(ublas_A, vcl_A, epsilon))
417  return EXIT_FAILURE;
418 
419  ublas_A *= double(alpha);
420  vcl_A *= double(alpha);
421 
422  if (!check_for_equality(ublas_A, vcl_A, epsilon))
423  return EXIT_FAILURE;
424 
425  std::cout << "Multiplication with GPU scalar: ";
426  ublas_A *= beta;
427  vcl_A *= gpu_beta;
428 
429  if (!check_for_equality(ublas_A, vcl_A, epsilon))
430  return EXIT_FAILURE;
431 
432 
433  std::cout << "Division with CPU scalar: ";
434  ublas_A /= long(alpha);
435  vcl_A /= long(alpha);
436 
437  if (!check_for_equality(ublas_A, vcl_A, epsilon))
438  return EXIT_FAILURE;
439 
440  ublas_A /= float(alpha);
441  vcl_A /= float(alpha);
442 
443  if (!check_for_equality(ublas_A, vcl_A, epsilon))
444  return EXIT_FAILURE;
445 
446  ublas_A /= double(alpha);
447  vcl_A /= double(alpha);
448 
449  if (!check_for_equality(ublas_A, vcl_A, epsilon))
450  return EXIT_FAILURE;
451 
452  std::cout << "Division with GPU scalar: ";
453  ublas_A /= beta;
454  vcl_A /= gpu_beta;
455 
456  if (!check_for_equality(ublas_A, vcl_A, epsilon))
457  return EXIT_FAILURE;
458 
459 
460 
461  std::cout << "Testing elementwise multiplication..." << std::endl;
462  ublas_B = boost::numeric::ublas::scalar_matrix<cpu_value_type>(ublas_B.size1(), ublas_B.size2(), cpu_value_type(1.4142));
463  ublas_A = cpu_value_type(3.1415) * ublas_B;
464  viennacl::copy(ublas_A, vcl_A);
465  viennacl::copy(ublas_B, vcl_B);
466  viennacl::copy(ublas_B, vcl_B);
467  ublas_A = boost::numeric::ublas::element_prod(ublas_A, ublas_B);
468  vcl_A = viennacl::linalg::element_prod(vcl_A, vcl_B);
469 
470  if (!check_for_equality(ublas_A, vcl_A, epsilon))
471  return EXIT_FAILURE;
472 
473  ublas_A += boost::numeric::ublas::element_prod(ublas_A, ublas_B);
474  vcl_A += viennacl::linalg::element_prod(vcl_A, vcl_B);
475 
476  if (!check_for_equality(ublas_A, vcl_A, epsilon))
477  return EXIT_FAILURE;
478 
479  ublas_A -= boost::numeric::ublas::element_prod(ublas_A, ublas_B);
480  vcl_A -= viennacl::linalg::element_prod(vcl_A, vcl_B);
481 
482  if (!check_for_equality(ublas_A, vcl_A, epsilon))
483  return EXIT_FAILURE;
484 
486  ublas_A = boost::numeric::ublas::element_prod(ublas_A + ublas_B, ublas_B);
487  vcl_A = viennacl::linalg::element_prod(vcl_A + vcl_B, vcl_B);
488 
489  if (!check_for_equality(ublas_A, vcl_A, epsilon))
490  return EXIT_FAILURE;
491 
492  ublas_A += boost::numeric::ublas::element_prod(ublas_A + ublas_B, ublas_B);
493  vcl_A += viennacl::linalg::element_prod(vcl_A + vcl_B, vcl_B);
494 
495  if (!check_for_equality(ublas_A, vcl_A, epsilon))
496  return EXIT_FAILURE;
497 
498  ublas_A -= boost::numeric::ublas::element_prod(ublas_A + ublas_B, ublas_B);
499  vcl_A -= viennacl::linalg::element_prod(vcl_A + vcl_B, vcl_B);
500 
501  if (!check_for_equality(ublas_A, vcl_A, epsilon))
502  return EXIT_FAILURE;
503 
505  ublas_A = boost::numeric::ublas::element_prod(ublas_A, ublas_B + ublas_A);
506  vcl_A = viennacl::linalg::element_prod(vcl_A, vcl_B + vcl_A);
507 
508  if (!check_for_equality(ublas_A, vcl_A, epsilon))
509  return EXIT_FAILURE;
510 
511  ublas_A += boost::numeric::ublas::element_prod(ublas_A, ublas_B + ublas_A);
512  vcl_A += viennacl::linalg::element_prod(vcl_A, vcl_B + vcl_A);
513 
514  if (!check_for_equality(ublas_A, vcl_A, epsilon))
515  return EXIT_FAILURE;
516 
517  ublas_A -= boost::numeric::ublas::element_prod(ublas_A, ublas_B + ublas_A);
518  vcl_A -= viennacl::linalg::element_prod(vcl_A, vcl_B + vcl_A);
519 
520  if (!check_for_equality(ublas_A, vcl_A, epsilon))
521  return EXIT_FAILURE;
522 
524  ublas_A = boost::numeric::ublas::element_prod(ublas_A + ublas_B, ublas_B + ublas_A);
525  vcl_A = viennacl::linalg::element_prod(vcl_A + vcl_B, vcl_B + vcl_A);
526 
527  if (!check_for_equality(ublas_A, vcl_A, epsilon))
528  return EXIT_FAILURE;
529 
530  ublas_A += boost::numeric::ublas::element_prod(ublas_A + ublas_B, ublas_B + ublas_A);
531  vcl_A += viennacl::linalg::element_prod(vcl_A + vcl_B, vcl_B + vcl_A);
532 
533  if (!check_for_equality(ublas_A, vcl_A, epsilon))
534  return EXIT_FAILURE;
535 
536  ublas_A -= boost::numeric::ublas::element_prod(ublas_A + ublas_B, ublas_B + ublas_A);
537  vcl_A -= viennacl::linalg::element_prod(vcl_A + vcl_B, vcl_B + vcl_A);
538 
539  if (!check_for_equality(ublas_A, vcl_A, epsilon))
540  return EXIT_FAILURE;
541 
542 
543  ublas_B = boost::numeric::ublas::scalar_matrix<cpu_value_type>(ublas_B.size1(), ublas_B.size2(), cpu_value_type(1.4142));
544  ublas_A = cpu_value_type(3.1415) * ublas_B;
545  viennacl::copy(ublas_A, vcl_A);
546  viennacl::copy(ublas_B, vcl_B);
547 
548  ublas_A = boost::numeric::ublas::element_div(ublas_A, ublas_B);
549  vcl_A = viennacl::linalg::element_div(vcl_A, vcl_B);
550 
551  if (!check_for_equality(ublas_A, vcl_A, epsilon))
552  return EXIT_FAILURE;
553 
554  ublas_A += boost::numeric::ublas::element_div(ublas_A, ublas_B);
555  vcl_A += viennacl::linalg::element_div(vcl_A, vcl_B);
556 
557  if (!check_for_equality(ublas_A, vcl_A, epsilon))
558  return EXIT_FAILURE;
559 
560  ublas_A -= boost::numeric::ublas::element_div(ublas_A, ublas_B);
561  vcl_A -= viennacl::linalg::element_div(vcl_A, vcl_B);
562 
563  if (!check_for_equality(ublas_A, vcl_A, epsilon))
564  return EXIT_FAILURE;
565 
567  ublas_A = boost::numeric::ublas::element_div(ublas_A + ublas_B, ublas_B);
568  vcl_A = viennacl::linalg::element_div(vcl_A + vcl_B, vcl_B);
569 
570  if (!check_for_equality(ublas_A, vcl_A, epsilon))
571  return EXIT_FAILURE;
572 
573  ublas_A += boost::numeric::ublas::element_div(ublas_A + ublas_B, ublas_B);
574  vcl_A += viennacl::linalg::element_div(vcl_A + vcl_B, vcl_B);
575 
576  if (!check_for_equality(ublas_A, vcl_A, epsilon))
577  return EXIT_FAILURE;
578 
579  ublas_A -= boost::numeric::ublas::element_div(ublas_A + ublas_B, ublas_B);
580  vcl_A -= viennacl::linalg::element_div(vcl_A + vcl_B, vcl_B);
581 
582  if (!check_for_equality(ublas_A, vcl_A, epsilon))
583  return EXIT_FAILURE;
584 
586  ublas_A = boost::numeric::ublas::element_div(ublas_A, ublas_B + ublas_A);
587  vcl_A = viennacl::linalg::element_div(vcl_A, vcl_B + vcl_A);
588 
589  if (!check_for_equality(ublas_A, vcl_A, epsilon))
590  return EXIT_FAILURE;
591 
592  ublas_A += boost::numeric::ublas::element_div(ublas_A, ublas_B + ublas_A);
593  vcl_A += viennacl::linalg::element_div(vcl_A, vcl_B + vcl_A);
594 
595  if (!check_for_equality(ublas_A, vcl_A, epsilon))
596  return EXIT_FAILURE;
597 
598  ublas_A -= boost::numeric::ublas::element_div(ublas_A, ublas_B + ublas_A);
599  vcl_A -= viennacl::linalg::element_div(vcl_A, vcl_B + vcl_A);
600 
601  if (!check_for_equality(ublas_A, vcl_A, epsilon))
602  return EXIT_FAILURE;
603 
605  ublas_A = boost::numeric::ublas::element_div(ublas_A + ublas_B, ublas_B + ublas_A);
606  vcl_A = viennacl::linalg::element_div(vcl_A + vcl_B, vcl_B + vcl_A);
607 
608  if (!check_for_equality(ublas_A, vcl_A, epsilon))
609  return EXIT_FAILURE;
610 
611  ublas_A += boost::numeric::ublas::element_div(ublas_A + ublas_B, ublas_B + ublas_A);
612  vcl_A += viennacl::linalg::element_div(vcl_A + vcl_B, vcl_B + vcl_A);
613 
614  if (!check_for_equality(ublas_A, vcl_A, epsilon))
615  return EXIT_FAILURE;
616 
617  ublas_A -= boost::numeric::ublas::element_div(ublas_A + ublas_B, ublas_B + ublas_A);
618  vcl_A -= viennacl::linalg::element_div(vcl_A + vcl_B, vcl_B + vcl_A);
619 
620  if (!check_for_equality(ublas_A, vcl_A, epsilon))
621  return EXIT_FAILURE;
622 
623  // element_pow
624  std::cout << "Testing unary element_pow()..." << std::endl;
625 
626  ublas_B = boost::numeric::ublas::scalar_matrix<cpu_value_type>(ublas_B.size1(), ublas_B.size2(), cpu_value_type(1.4142));
627  ublas_A = cpu_value_type(3.1415) * ublas_B;
628  viennacl::copy(ublas_A, vcl_A);
629  viennacl::copy(ublas_B, vcl_B);
630 
631  for (std::size_t i=0; i<ublas_C.size1(); i++)
632  for (std::size_t j=0; j<ublas_C.size2(); ++j)
633  ublas_C(i,j) = std::pow(ublas_A(i,j), ublas_B(i,j));
634  vcl_C = viennacl::linalg::element_pow(vcl_A, vcl_B);
635 
636  if (!check_for_equality(ublas_C, vcl_C, epsilon))
637  return EXIT_FAILURE;
638 
639  for (std::size_t i=0; i<ublas_C.size1(); i++)
640  for (std::size_t j=0; j<ublas_C.size2(); ++j)
641  ublas_C(i,j) += std::pow(ublas_A(i,j), ublas_B(i,j));
642  vcl_C += viennacl::linalg::element_pow(vcl_A, vcl_B);
643 
644  if (!check_for_equality(ublas_C, vcl_C, epsilon))
645  return EXIT_FAILURE;
646 
647  for (std::size_t i=0; i<ublas_C.size1(); i++)
648  for (std::size_t j=0; j<ublas_C.size2(); ++j)
649  ublas_C(i,j) -= std::pow(ublas_A(i,j), ublas_B(i,j));
650  vcl_C -= viennacl::linalg::element_pow(vcl_A, vcl_B);
651 
652  if (!check_for_equality(ublas_C, vcl_C, epsilon))
653  return EXIT_FAILURE;
654 
656  for (std::size_t i=0; i<ublas_C.size1(); i++)
657  for (std::size_t j=0; j<ublas_C.size2(); ++j)
658  ublas_C(i,j) = std::pow(ublas_A(i,j) + ublas_B(i,j), ublas_B(i,j));
659  vcl_C = viennacl::linalg::element_pow(vcl_A + vcl_B, vcl_B);
660 
661  if (!check_for_equality(ublas_C, vcl_C, epsilon))
662  return EXIT_FAILURE;
663 
664  for (std::size_t i=0; i<ublas_C.size1(); i++)
665  for (std::size_t j=0; j<ublas_C.size2(); ++j)
666  ublas_C(i,j) += std::pow(ublas_A(i,j) + ublas_B(i,j), ublas_B(i,j));
667  vcl_C += viennacl::linalg::element_pow(vcl_A + vcl_B, vcl_B);
668 
669  if (!check_for_equality(ublas_C, vcl_C, epsilon))
670  return EXIT_FAILURE;
671 
672  for (std::size_t i=0; i<ublas_C.size1(); i++)
673  for (std::size_t j=0; j<ublas_C.size2(); ++j)
674  ublas_C(i,j) -= std::pow(ublas_A(i,j) + ublas_B(i,j), ublas_B(i,j));
675  vcl_C -= viennacl::linalg::element_pow(vcl_A + vcl_B, vcl_B);
676 
677  if (!check_for_equality(ublas_C, vcl_C, epsilon))
678  return EXIT_FAILURE;
679 
681  for (std::size_t i=0; i<ublas_C.size1(); i++)
682  for (std::size_t j=0; j<ublas_C.size2(); ++j)
683  ublas_C(i,j) = std::pow(ublas_A(i,j), ublas_B(i,j) + ublas_A(i,j));
684  vcl_C = viennacl::linalg::element_pow(vcl_A, vcl_B + vcl_A);
685 
686  if (!check_for_equality(ublas_C, vcl_C, epsilon))
687  return EXIT_FAILURE;
688 
689  for (std::size_t i=0; i<ublas_C.size1(); i++)
690  for (std::size_t j=0; j<ublas_C.size2(); ++j)
691  ublas_C(i,j) += std::pow(ublas_A(i,j), ublas_B(i,j) + ublas_A(i,j));
692  vcl_C += viennacl::linalg::element_pow(vcl_A, vcl_B + vcl_A);
693 
694  if (!check_for_equality(ublas_C, vcl_C, epsilon))
695  return EXIT_FAILURE;
696 
697  for (std::size_t i=0; i<ublas_C.size1(); i++)
698  for (std::size_t j=0; j<ublas_C.size2(); ++j)
699  ublas_C(i,j) -= std::pow(ublas_A(i,j), ublas_B(i,j) + ublas_A(i,j));
700  vcl_C -= viennacl::linalg::element_pow(vcl_A, vcl_B + vcl_A);
701 
702  if (!check_for_equality(ublas_C, vcl_C, epsilon))
703  return EXIT_FAILURE;
704 
706  for (std::size_t i=0; i<ublas_C.size1(); i++)
707  for (std::size_t j=0; j<ublas_C.size2(); ++j)
708  ublas_C(i,j) = std::pow(ublas_A(i,j) + ublas_B(i,j), ublas_B(i,j) + ublas_A(i,j));
709  vcl_C = viennacl::linalg::element_pow(vcl_A + vcl_B, vcl_B + vcl_A);
710 
711  if (!check_for_equality(ublas_C, vcl_C, epsilon))
712  return EXIT_FAILURE;
713 
714  for (std::size_t i=0; i<ublas_C.size1(); i++)
715  for (std::size_t j=0; j<ublas_C.size2(); ++j)
716  ublas_C(i,j) += std::pow(ublas_A(i,j) + ublas_B(i,j), ublas_B(i,j) + ublas_A(i,j));
717  vcl_C += viennacl::linalg::element_pow(vcl_A + vcl_B, vcl_B + vcl_A);
718 
719  if (!check_for_equality(ublas_C, vcl_C, epsilon))
720  return EXIT_FAILURE;
721 
722  for (std::size_t i=0; i<ublas_C.size1(); i++)
723  for (std::size_t j=0; j<ublas_C.size2(); ++j)
724  ublas_C(i,j) -= std::pow(ublas_A(i,j) + ublas_B(i,j), ublas_B(i,j) + ublas_A(i,j));
725  vcl_C -= viennacl::linalg::element_pow(vcl_A + vcl_B, vcl_B + vcl_A);
726 
727  if (!check_for_equality(ublas_C, vcl_C, epsilon))
728  return EXIT_FAILURE;
729 
730 
731  std::cout << "Testing unary elementwise operations..." << std::endl;
732 
733 #define GENERATE_UNARY_OP_TEST(FUNCNAME) \
734  ublas_B = boost::numeric::ublas::scalar_matrix<cpu_value_type>(ublas_B.size1(), ublas_B.size2(), cpu_value_type(1.4142)); \
735  ublas_A = cpu_value_type(3.1415) * ublas_B; \
736  ublas_C = cpu_value_type(2.7172) * ublas_A; \
737  viennacl::copy(ublas_A, vcl_A); \
738  viennacl::copy(ublas_B, vcl_B); \
739  viennacl::copy(ublas_C, vcl_C); \
740  viennacl::copy(ublas_B, vcl_B); \
741  \
742  for (std::size_t i=0; i<ublas_C.size1(); ++i) \
743  for (std::size_t j=0; j<ublas_C.size2(); ++j) \
744  ublas_C(i,j) = std::FUNCNAME(ublas_A(i,j)); \
745  vcl_C = viennacl::linalg::element_##FUNCNAME(vcl_A); \
746  \
747  if (!check_for_equality(ublas_C, vcl_C, epsilon)) \
748  { \
749  std::cout << "Failure at C = " << #FUNCNAME << "(A)" << std::endl; \
750  return EXIT_FAILURE; \
751  } \
752  \
753  for (std::size_t i=0; i<ublas_C.size1(); ++i) \
754  for (std::size_t j=0; j<ublas_C.size2(); ++j) \
755  ublas_C(i,j) = std::FUNCNAME(ublas_A(i,j) + ublas_B(i,j)); \
756  vcl_C = viennacl::linalg::element_##FUNCNAME(vcl_A + vcl_B); \
757  \
758  if (!check_for_equality(ublas_C, vcl_C, epsilon)) \
759  { \
760  std::cout << "Failure at C = " << #FUNCNAME << "(A + B)" << std::endl; \
761  return EXIT_FAILURE; \
762  } \
763  \
764  for (std::size_t i=0; i<ublas_C.size1(); ++i) \
765  for (std::size_t j=0; j<ublas_C.size2(); ++j) \
766  ublas_C(i,j) += std::FUNCNAME(ublas_A(i,j)); \
767  vcl_C += viennacl::linalg::element_##FUNCNAME(vcl_A); \
768  \
769  if (!check_for_equality(ublas_C, vcl_C, epsilon)) \
770  { \
771  std::cout << "Failure at C += " << #FUNCNAME << "(A)" << std::endl; \
772  return EXIT_FAILURE; \
773  } \
774  \
775  for (std::size_t i=0; i<ublas_C.size1(); ++i) \
776  for (std::size_t j=0; j<ublas_C.size2(); ++j) \
777  ublas_C(i,j) += std::FUNCNAME(ublas_A(i,j) + ublas_B(i,j)); \
778  vcl_C += viennacl::linalg::element_##FUNCNAME(vcl_A + vcl_B); \
779  \
780  if (!check_for_equality(ublas_C, vcl_C, epsilon)) \
781  { \
782  std::cout << "Failure at C += " << #FUNCNAME << "(A + B)" << std::endl; \
783  return EXIT_FAILURE; \
784  } \
785  \
786  for (std::size_t i=0; i<ublas_C.size1(); ++i) \
787  for (std::size_t j=0; j<ublas_C.size2(); ++j) \
788  ublas_C(i,j) -= std::FUNCNAME(ublas_A(i,j)); \
789  vcl_C -= viennacl::linalg::element_##FUNCNAME(vcl_A); \
790  \
791  if (!check_for_equality(ublas_C, vcl_C, epsilon)) \
792  { \
793  std::cout << "Failure at C -= " << #FUNCNAME << "(A)" << std::endl; \
794  return EXIT_FAILURE; \
795  } \
796  \
797  for (std::size_t i=0; i<ublas_C.size1(); ++i) \
798  for (std::size_t j=0; j<ublas_C.size2(); ++j) \
799  ublas_C(i,j) -= std::FUNCNAME(ublas_A(i,j) + ublas_B(i,j)); \
800  vcl_C -= viennacl::linalg::element_##FUNCNAME(vcl_A + vcl_B); \
801  \
802  if (!check_for_equality(ublas_C, vcl_C, epsilon)) \
803  { \
804  std::cout << "Failure at C -= " << #FUNCNAME << "(A + B)" << std::endl; \
805  return EXIT_FAILURE; \
806  } \
807  \
808 
809  GENERATE_UNARY_OP_TEST(cos);
810  GENERATE_UNARY_OP_TEST(cosh);
811  GENERATE_UNARY_OP_TEST(exp);
812  GENERATE_UNARY_OP_TEST(floor);
813  GENERATE_UNARY_OP_TEST(fabs);
814  GENERATE_UNARY_OP_TEST(log);
815  GENERATE_UNARY_OP_TEST(log10);
816  GENERATE_UNARY_OP_TEST(sin);
817  GENERATE_UNARY_OP_TEST(sinh);
818  GENERATE_UNARY_OP_TEST(fabs);
819  GENERATE_UNARY_OP_TEST(sqrt);
820  GENERATE_UNARY_OP_TEST(tan);
821  GENERATE_UNARY_OP_TEST(tanh);
822 
823  std::cout << "Complicated expressions: ";
824  //std::cout << "ublas_A: " << ublas_A << std::endl;
825  //std::cout << "ublas_B: " << ublas_B << std::endl;
826  //std::cout << "ublas_C: " << ublas_C << std::endl;
827  ublas_B += alpha * (- ublas_A - beta * ublas_C + ublas_A);
828  vcl_B += gpu_alpha * (- vcl_A - vcl_C * beta + vcl_A);
829 
830  if (!check_for_equality(ublas_B, vcl_B, epsilon))
831  return EXIT_FAILURE;
832 
833  ublas_B += (- ublas_A - beta * ublas_C + ublas_A * beta) / gpu_alpha;
834  vcl_B += (- vcl_A - vcl_C * beta + gpu_beta * vcl_A) / gpu_alpha;
835 
836  if (!check_for_equality(ublas_B, vcl_B, epsilon))
837  return EXIT_FAILURE;
838 
839 
840  ublas_B -= alpha * (- ublas_A - beta * ublas_C - ublas_A);
841  vcl_B -= gpu_alpha * (- vcl_A - vcl_C * beta - vcl_A);
842 
843  if (!check_for_equality(ublas_B, vcl_B, epsilon))
844  return EXIT_FAILURE;
845 
846  ublas_B -= (- ublas_A - beta * ublas_C - ublas_A * beta) / alpha;
847  vcl_B -= (- vcl_A - vcl_C * beta - gpu_beta * vcl_A) / gpu_alpha;
848 
849  if (!check_for_equality(ublas_B, vcl_B, epsilon))
850  return EXIT_FAILURE;
851 
852  std::cout << std::endl;
853  std::cout << "----------------------------------------------" << std::endl;
854  std::cout << std::endl;
855 
856 
857  return EXIT_SUCCESS;
858 }
859 
860 
861 
862 
863 template<typename T, typename ScalarType>
864 int run_test(double epsilon)
865 {
866  //typedef float ScalarType;
867  typedef boost::numeric::ublas::matrix<ScalarType> MatrixType;
868 
869  typedef viennacl::matrix<ScalarType, T> VCLMatrixType;
870 
871  std::size_t dim_rows = 131;
872  std::size_t dim_cols = 33;
873  //std::size_t dim_rows = 5;
874  //std::size_t dim_cols = 3;
875 
876  //setup ublas objects:
877  MatrixType ublas_A(dim_rows, dim_cols);
878  MatrixType ublas_B(dim_rows, dim_cols);
879  MatrixType ublas_C(dim_rows, dim_cols);
880 
881  for (std::size_t i=0; i<ublas_A.size1(); ++i)
882  for (std::size_t j=0; j<ublas_A.size2(); ++j)
883  {
884  ublas_A(i,j) = ScalarType((i+2) + (j+1)*(i+2));
885  ublas_B(i,j) = ScalarType((j+2) + (j+1)*(j+2));
886  ublas_C(i,j) = ScalarType((i+1) + (i+1)*(i+2));
887  }
888 
889  MatrixType ublas_A_large(4 * dim_rows, 4 * dim_cols);
890  for (std::size_t i=0; i<ublas_A_large.size1(); ++i)
891  for (std::size_t j=0; j<ublas_A_large.size2(); ++j)
892  ublas_A_large(i,j) = ScalarType(i * ublas_A_large.size2() + j);
893 
894  //Setup ViennaCL objects
895  VCLMatrixType vcl_A_full(4 * dim_rows, 4 * dim_cols);
896  VCLMatrixType vcl_B_full(4 * dim_rows, 4 * dim_cols);
897  VCLMatrixType vcl_C_full(4 * dim_rows, 4 * dim_cols);
898 
899  viennacl::copy(ublas_A_large, vcl_A_full);
900  viennacl::copy(ublas_A_large, vcl_B_full);
901  viennacl::copy(ublas_A_large, vcl_C_full);
902 
903  //
904  // Create A
905  //
906  VCLMatrixType vcl_A(dim_rows, dim_cols);
907 
908  viennacl::range vcl_A_r1(2 * dim_rows, 3 * dim_rows);
909  viennacl::range vcl_A_r2(dim_cols, 2 * dim_cols);
910  viennacl::matrix_range<VCLMatrixType> vcl_range_A(vcl_A_full, vcl_A_r1, vcl_A_r2);
911 
912  viennacl::slice vcl_A_s1(2, 3, dim_rows);
913  viennacl::slice vcl_A_s2(2 * dim_cols, 2, dim_cols);
914  viennacl::matrix_slice<VCLMatrixType> vcl_slice_A(vcl_A_full, vcl_A_s1, vcl_A_s2);
915 
916 
917  //
918  // Create B
919  //
920  VCLMatrixType vcl_B(dim_rows, dim_cols);
921 
922  viennacl::range vcl_B_r1(dim_rows, 2 * dim_rows);
923  viennacl::range vcl_B_r2(2 * dim_cols, 3 * dim_cols);
924  viennacl::matrix_range<VCLMatrixType> vcl_range_B(vcl_B_full, vcl_B_r1, vcl_B_r2);
925 
926  viennacl::slice vcl_B_s1(2 * dim_rows, 2, dim_rows);
927  viennacl::slice vcl_B_s2(dim_cols, 3, dim_cols);
928  viennacl::matrix_slice<VCLMatrixType> vcl_slice_B(vcl_B_full, vcl_B_s1, vcl_B_s2);
929 
930 
931  //
932  // Create C
933  //
934  VCLMatrixType vcl_C(dim_rows, dim_cols);
935 
936  viennacl::range vcl_C_r1(2 * dim_rows, 3 * dim_rows);
937  viennacl::range vcl_C_r2(3 * dim_cols, 4 * dim_cols);
938  viennacl::matrix_range<VCLMatrixType> vcl_range_C(vcl_C_full, vcl_C_r1, vcl_C_r2);
939 
940  viennacl::slice vcl_C_s1(dim_rows, 2, dim_rows);
941  viennacl::slice vcl_C_s2(0, 3, dim_cols);
942  viennacl::matrix_slice<VCLMatrixType> vcl_slice_C(vcl_C_full, vcl_C_s1, vcl_C_s2);
943 
944  viennacl::copy(ublas_A, vcl_A);
945  viennacl::copy(ublas_A, vcl_range_A);
946  viennacl::copy(ublas_A, vcl_slice_A);
947 
948  viennacl::copy(ublas_B, vcl_B);
949  viennacl::copy(ublas_B, vcl_range_B);
950  viennacl::copy(ublas_B, vcl_slice_B);
951 
952  viennacl::copy(ublas_C, vcl_C);
953  viennacl::copy(ublas_C, vcl_range_C);
954  viennacl::copy(ublas_C, vcl_slice_C);
955 
956 
957  std::cout << std::endl;
958  std::cout << "//" << std::endl;
959  std::cout << "////////// Test: Copy CTOR //////////" << std::endl;
960  std::cout << "//" << std::endl;
961 
962  {
963  std::cout << "Testing matrix created from range... ";
964  VCLMatrixType vcl_temp = vcl_range_A;
965  if (check_for_equality(ublas_A, vcl_temp, epsilon))
966  std::cout << "PASSED!" << std::endl;
967  else
968  {
969  std::cout << "ublas_A: " << ublas_A << std::endl;
970  std::cout << "vcl_temp: " << vcl_temp << std::endl;
971  std::cout << "vcl_range_A: " << vcl_range_A << std::endl;
972  std::cout << "vcl_A: " << vcl_A << std::endl;
973  std::cout << std::endl << "TEST failed!" << std::endl;
974  return EXIT_FAILURE;
975  }
976 
977  std::cout << "Testing matrix created from slice... ";
978  VCLMatrixType vcl_temp2 = vcl_range_B;
979  if (check_for_equality(ublas_B, vcl_temp2, epsilon))
980  std::cout << "PASSED!" << std::endl;
981  else
982  {
983  std::cout << std::endl << "TEST failed!" << std::endl;
984  return EXIT_FAILURE;
985  }
986  }
987 
988  std::cout << "//" << std::endl;
989  std::cout << "////////// Test: Initializer for matrix type //////////" << std::endl;
990  std::cout << "//" << std::endl;
991 
992  {
993  boost::numeric::ublas::matrix<ScalarType> ublas_dummy1 = boost::numeric::ublas::identity_matrix<ScalarType>(ublas_A.size1());
994  boost::numeric::ublas::matrix<ScalarType> ublas_dummy2 = boost::numeric::ublas::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
995  boost::numeric::ublas::matrix<ScalarType> ublas_dummy3 = boost::numeric::ublas::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
996 
997  viennacl::matrix<ScalarType> vcl_dummy1 = viennacl::identity_matrix<ScalarType>(ublas_A.size1());
998  viennacl::matrix<ScalarType> vcl_dummy2 = viennacl::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
999  viennacl::matrix<ScalarType> vcl_dummy3 = viennacl::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
1000 
1001  std::cout << "Testing initializer CTOR... ";
1002  if ( check_for_equality(ublas_dummy1, vcl_dummy1, epsilon)
1003  && check_for_equality(ublas_dummy2, vcl_dummy2, epsilon)
1004  && check_for_equality(ublas_dummy3, vcl_dummy3, epsilon)
1005  )
1006  std::cout << "PASSED!" << std::endl;
1007  else
1008  {
1009  std::cout << std::endl << "TEST failed!" << std::endl;
1010  return EXIT_FAILURE;
1011  }
1012 
1013  ublas_dummy1 = boost::numeric::ublas::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
1014  ublas_dummy2 = boost::numeric::ublas::identity_matrix<ScalarType>(ublas_A.size1());
1015  ublas_dummy3 = boost::numeric::ublas::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
1016 
1017  vcl_dummy1 = viennacl::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
1018  vcl_dummy2 = viennacl::identity_matrix<ScalarType>(ublas_A.size1());
1019  vcl_dummy3 = viennacl::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
1020 
1021  std::cout << "Testing initializer assignment... ";
1022  if ( check_for_equality(ublas_dummy1, vcl_dummy1, epsilon)
1023  && check_for_equality(ublas_dummy2, vcl_dummy2, epsilon)
1024  && check_for_equality(ublas_dummy3, vcl_dummy3, epsilon)
1025  )
1026  std::cout << "PASSED!" << std::endl;
1027  else
1028  {
1029  std::cout << std::endl << "TEST failed!" << std::endl;
1030  return EXIT_FAILURE;
1031  }
1032  }
1033 
1034  std::cout << "//" << std::endl;
1035  std::cout << "////////// Test: Norms //////////" << std::endl;
1036  std::cout << "//" << std::endl;
1037 
1038  /*ScalarType ublas_norm_1 = viennacl::linalg::norm_1(ublas_C);
1039  ScalarType vcl_norm_1 = viennacl::linalg::norm_1(vcl_C);
1040  if ( std::fabs(ublas_norm_1 - vcl_norm_1) / ublas_norm_1 > epsilon)
1041  {
1042  std::cerr << "Failure at norm_1(): " << std::fabs(ublas_norm_1 - vcl_norm_1) / ublas_norm_1 << std::endl;
1043  return EXIT_FAILURE;
1044  }
1045 
1046  ScalarType ublas_norm_inf = boost::numeric::ublas::norm_inf(ublas_C);
1047  ScalarType vcl_norm_inf = viennacl::linalg::norm_inf(vcl_C);
1048  if ( std::fabs(ublas_norm_inf - vcl_norm_inf) / ublas_norm_inf > epsilon)
1049  {
1050  std::cerr << "Failure at norm_inf(): " << std::fabs(ublas_norm_inf - vcl_norm_inf) / ublas_norm_inf << std::endl;
1051  return EXIT_FAILURE;
1052  }*/
1053 
1054  ScalarType ublas_norm_frobenius = viennacl::linalg::norm_frobenius(ublas_C);
1055  ScalarType vcl_norm_frobenius = viennacl::linalg::norm_frobenius(vcl_C);
1056  if ( std::fabs(ublas_norm_frobenius - vcl_norm_frobenius) / ublas_norm_frobenius > epsilon)
1057  {
1058  std::cerr << "Failure at norm_frobenius()" << std::endl;
1059  return EXIT_FAILURE;
1060  }
1061 
1062  viennacl::scalar<ScalarType> device_ublas_norm_frobenius = viennacl::linalg::norm_frobenius(ublas_C);
1063  viennacl::scalar<ScalarType> device_vcl_norm_frobenius = viennacl::linalg::norm_frobenius(vcl_C);
1064  if ( std::fabs(device_ublas_norm_frobenius - device_vcl_norm_frobenius) / device_ublas_norm_frobenius > epsilon)
1065  {
1066  std::cerr << "Failure at norm_frobenius()" << std::endl;
1067  return EXIT_FAILURE;
1068  }
1069 
1070  vcl_norm_frobenius = viennacl::linalg::norm_frobenius(vcl_range_C);
1071  if ( std::fabs(ublas_norm_frobenius - vcl_norm_frobenius) / ublas_norm_frobenius > epsilon)
1072  {
1073  std::cerr << "Failure at norm_frobenius() with range" << std::endl;
1074  return EXIT_FAILURE;
1075  }
1076 
1077  device_vcl_norm_frobenius = viennacl::linalg::norm_frobenius(vcl_range_C);
1078  if ( std::fabs(device_ublas_norm_frobenius - device_vcl_norm_frobenius) / device_ublas_norm_frobenius > epsilon)
1079  {
1080  std::cerr << "Failure at norm_frobenius() with range" << std::endl;
1081  return EXIT_FAILURE;
1082  }
1083 
1084  vcl_norm_frobenius = viennacl::linalg::norm_frobenius(vcl_slice_C);
1085  if ( std::fabs(ublas_norm_frobenius - vcl_norm_frobenius) / ublas_norm_frobenius > epsilon)
1086  {
1087  std::cerr << "Failure at norm_frobenius() with slice" << std::endl;
1088  return EXIT_FAILURE;
1089  }
1090 
1091  device_vcl_norm_frobenius = viennacl::linalg::norm_frobenius(vcl_slice_C);
1092  if ( std::fabs(device_ublas_norm_frobenius - device_vcl_norm_frobenius) / device_ublas_norm_frobenius > epsilon)
1093  {
1094  std::cerr << "Failure at norm_frobenius() with slice" << std::endl;
1095  return EXIT_FAILURE;
1096  }
1097 
1098  std::cout << "PASSED!" << std::endl;
1099 
1100 
1101  //
1102  // run operation tests:
1103  //
1104 
1105  std::cout << "//" << std::endl;
1106  std::cout << "////////// Test: Operations //////////" << std::endl;
1107  std::cout << "//" << std::endl;
1108 
1109 
1111  std::cout << "Testing A=matrix, B=matrix, C=matrix ..." << std::endl;
1112  viennacl::copy(ublas_A, vcl_A);
1113  viennacl::copy(ublas_B, vcl_B);
1114  viennacl::copy(ublas_C, vcl_C);
1115  if (run_test(epsilon,
1116  ublas_A, ublas_B, ublas_C,
1117  vcl_A, vcl_B, vcl_C) != EXIT_SUCCESS)
1118  {
1119  return EXIT_FAILURE;
1120  }
1121 
1122  std::cout << "Testing A=matrix, B=matrix, C=range ..." << std::endl;
1123  viennacl::copy(ublas_A, vcl_A);
1124  viennacl::copy(ublas_B, vcl_B);
1125  viennacl::copy(ublas_C, vcl_range_C);
1126  if (run_test(epsilon,
1127  ublas_A, ublas_B, ublas_C,
1128  vcl_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
1129  {
1130  return EXIT_FAILURE;
1131  }
1132 
1133  std::cout << "Testing A=matrix, B=matrix, C=slice ..." << std::endl;
1134  viennacl::copy(ublas_A, vcl_A);
1135  viennacl::copy(ublas_B, vcl_B);
1136  viennacl::copy(ublas_C, vcl_slice_C);
1137  if (run_test(epsilon,
1138  ublas_A, ublas_B, ublas_C,
1139  vcl_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
1140  {
1141  return EXIT_FAILURE;
1142  }
1143 
1144  std::cout << "Testing A=matrix, B=range, C=matrix ..." << std::endl;
1145  viennacl::copy(ublas_A, vcl_A);
1146  viennacl::copy(ublas_B, vcl_range_B);
1147  viennacl::copy(ublas_C, vcl_C);
1148  if (run_test(epsilon,
1149  ublas_A, ublas_B, ublas_C,
1150  vcl_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
1151  {
1152  return EXIT_FAILURE;
1153  }
1154 
1155  std::cout << "Testing A=matrix, B=range, C=range ..." << std::endl;
1156  viennacl::copy(ublas_A, vcl_A);
1157  viennacl::copy(ublas_B, vcl_range_B);
1158  viennacl::copy(ublas_C, vcl_range_C);
1159  if (run_test(epsilon,
1160  ublas_A, ublas_B, ublas_C,
1161  vcl_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
1162  {
1163  return EXIT_FAILURE;
1164  }
1165 
1166  std::cout << "Testing A=matrix, B=range, C=slice ..." << std::endl;
1167  viennacl::copy(ublas_A, vcl_A);
1168  viennacl::copy(ublas_B, vcl_range_B);
1169  viennacl::copy(ublas_C, vcl_slice_C);
1170  if (run_test(epsilon,
1171  ublas_A, ublas_B, ublas_C,
1172  vcl_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
1173  {
1174  return EXIT_FAILURE;
1175  }
1176 
1177 
1178  std::cout << "Testing A=matrix, B=slice, C=matrix ..." << std::endl;
1179  viennacl::copy(ublas_A, vcl_A);
1180  viennacl::copy(ublas_B, vcl_slice_B);
1181  viennacl::copy(ublas_C, vcl_C);
1182  if (run_test(epsilon,
1183  ublas_A, ublas_B, ublas_C,
1184  vcl_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
1185  {
1186  return EXIT_FAILURE;
1187  }
1188 
1189  std::cout << "Testing A=matrix, B=slice, C=range ..." << std::endl;
1190  viennacl::copy(ublas_A, vcl_A);
1191  viennacl::copy(ublas_B, vcl_slice_B);
1192  viennacl::copy(ublas_C, vcl_range_C);
1193  if (run_test(epsilon,
1194  ublas_A, ublas_B, ublas_C,
1195  vcl_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
1196  {
1197  return EXIT_FAILURE;
1198  }
1199 
1200  std::cout << "Testing A=matrix, B=slice, C=slice ..." << std::endl;
1201  viennacl::copy(ublas_A, vcl_A);
1202  viennacl::copy(ublas_B, vcl_slice_B);
1203  viennacl::copy(ublas_C, vcl_slice_C);
1204  if (run_test(epsilon,
1205  ublas_A, ublas_B, ublas_C,
1206  vcl_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
1207  {
1208  return EXIT_FAILURE;
1209  }
1210 
1211 
1212 
1214  std::cout << "Testing A=range, B=matrix, C=matrix ..." << std::endl;
1215  viennacl::copy(ublas_A, vcl_range_A);
1216  viennacl::copy(ublas_B, vcl_B);
1217  viennacl::copy(ublas_C, vcl_C);
1218  if (run_test(epsilon,
1219  ublas_A, ublas_B, ublas_C,
1220  vcl_range_A, vcl_B, vcl_C) != EXIT_SUCCESS)
1221  {
1222  return EXIT_FAILURE;
1223  }
1224 
1225  std::cout << "Testing A=range, B=matrix, C=range ..." << std::endl;
1226  viennacl::copy(ublas_A, vcl_range_A);
1227  viennacl::copy(ublas_B, vcl_B);
1228  viennacl::copy(ublas_C, vcl_range_C);
1229  if (run_test(epsilon,
1230  ublas_A, ublas_B, ublas_C,
1231  vcl_range_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
1232  {
1233  return EXIT_FAILURE;
1234  }
1235 
1236  std::cout << "Testing A=range, B=matrix, C=slice ..." << std::endl;
1237  viennacl::copy(ublas_A, vcl_range_A);
1238  viennacl::copy(ublas_B, vcl_B);
1239  viennacl::copy(ublas_C, vcl_slice_C);
1240  if (run_test(epsilon,
1241  ublas_A, ublas_B, ublas_C,
1242  vcl_range_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
1243  {
1244  return EXIT_FAILURE;
1245  }
1246 
1247 
1248 
1249  std::cout << "Testing A=range, B=range, C=matrix ..." << std::endl;
1250  viennacl::copy(ublas_A, vcl_range_A);
1251  viennacl::copy(ublas_B, vcl_range_B);
1252  viennacl::copy(ublas_C, vcl_C);
1253  if (run_test(epsilon,
1254  ublas_A, ublas_B, ublas_C,
1255  vcl_range_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
1256  {
1257  return EXIT_FAILURE;
1258  }
1259 
1260  std::cout << "Testing A=range, B=range, C=range ..." << std::endl;
1261  viennacl::copy(ublas_A, vcl_range_A);
1262  viennacl::copy(ublas_B, vcl_range_B);
1263  viennacl::copy(ublas_C, vcl_range_C);
1264  if (run_test(epsilon,
1265  ublas_A, ublas_B, ublas_C,
1266  vcl_range_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
1267  {
1268  return EXIT_FAILURE;
1269  }
1270 
1271  std::cout << "Testing A=range, B=range, C=slice ..." << std::endl;
1272  viennacl::copy(ublas_A, vcl_range_A);
1273  viennacl::copy(ublas_B, vcl_range_B);
1274  viennacl::copy(ublas_C, vcl_slice_C);
1275  if (run_test(epsilon,
1276  ublas_A, ublas_B, ublas_C,
1277  vcl_range_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
1278  {
1279  return EXIT_FAILURE;
1280  }
1281 
1282 
1283 
1284  std::cout << "Testing A=range, B=slice, C=matrix ..." << std::endl;
1285  viennacl::copy(ublas_A, vcl_range_A);
1286  viennacl::copy(ublas_B, vcl_slice_B);
1287  viennacl::copy(ublas_C, vcl_C);
1288  if (run_test(epsilon,
1289  ublas_A, ublas_B, ublas_C,
1290  vcl_range_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
1291  {
1292  return EXIT_FAILURE;
1293  }
1294 
1295  std::cout << "Testing A=range, B=slice, C=range ..." << std::endl;
1296  viennacl::copy(ublas_A, vcl_range_A);
1297  viennacl::copy(ublas_B, vcl_slice_B);
1298  viennacl::copy(ublas_C, vcl_range_C);
1299  if (run_test(epsilon,
1300  ublas_A, ublas_B, ublas_C,
1301  vcl_range_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
1302  {
1303  return EXIT_FAILURE;
1304  }
1305 
1306  std::cout << "Testing A=range, B=slice, C=slice ..." << std::endl;
1307  viennacl::copy(ublas_A, vcl_range_A);
1308  viennacl::copy(ublas_B, vcl_slice_B);
1309  viennacl::copy(ublas_C, vcl_slice_C);
1310  if (run_test(epsilon,
1311  ublas_A, ublas_B, ublas_C,
1312  vcl_range_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
1313  {
1314  return EXIT_FAILURE;
1315  }
1316 
1317 
1319  std::cout << "Testing A=slice, B=matrix, C=matrix ..." << std::endl;
1320  viennacl::copy(ublas_A, vcl_slice_A);
1321  viennacl::copy(ublas_B, vcl_B);
1322  viennacl::copy(ublas_C, vcl_C);
1323  if (run_test(epsilon,
1324  ublas_A, ublas_B, ublas_C,
1325  vcl_slice_A, vcl_B, vcl_C) != EXIT_SUCCESS)
1326  {
1327  return EXIT_FAILURE;
1328  }
1329 
1330  std::cout << "Testing A=slice, B=matrix, C=range ..." << std::endl;
1331  viennacl::copy(ublas_A, vcl_slice_A);
1332  viennacl::copy(ublas_B, vcl_B);
1333  viennacl::copy(ublas_C, vcl_range_C);
1334  if (run_test(epsilon,
1335  ublas_A, ublas_B, ublas_C,
1336  vcl_slice_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
1337  {
1338  return EXIT_FAILURE;
1339  }
1340 
1341  std::cout << "Testing A=slice, B=matrix, C=slice ..." << std::endl;
1342  viennacl::copy(ublas_A, vcl_slice_A);
1343  viennacl::copy(ublas_B, vcl_B);
1344  viennacl::copy(ublas_C, vcl_slice_C);
1345  if (run_test(epsilon,
1346  ublas_A, ublas_B, ublas_C,
1347  vcl_slice_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
1348  {
1349  return EXIT_FAILURE;
1350  }
1351 
1352 
1353 
1354  std::cout << "Testing A=slice, B=range, C=matrix ..." << std::endl;
1355  viennacl::copy(ublas_A, vcl_slice_A);
1356  viennacl::copy(ublas_B, vcl_range_B);
1357  viennacl::copy(ublas_C, vcl_C);
1358  if (run_test(epsilon,
1359  ublas_A, ublas_B, ublas_C,
1360  vcl_slice_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
1361  {
1362  return EXIT_FAILURE;
1363  }
1364 
1365  std::cout << "Testing A=slice, B=range, C=range ..." << std::endl;
1366  viennacl::copy(ublas_A, vcl_slice_A);
1367  viennacl::copy(ublas_B, vcl_range_B);
1368  viennacl::copy(ublas_C, vcl_range_C);
1369  if (run_test(epsilon,
1370  ublas_A, ublas_B, ublas_C,
1371  vcl_slice_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
1372  {
1373  return EXIT_FAILURE;
1374  }
1375 
1376  std::cout << "Testing A=slice, B=range, C=slice ..." << std::endl;
1377  viennacl::copy(ublas_A, vcl_slice_A);
1378  viennacl::copy(ublas_B, vcl_range_B);
1379  viennacl::copy(ublas_C, vcl_slice_C);
1380  if (run_test(epsilon,
1381  ublas_A, ublas_B, ublas_C,
1382  vcl_slice_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
1383  {
1384  return EXIT_FAILURE;
1385  }
1386 
1387 
1388 
1389  std::cout << "Testing A=slice, B=slice, C=matrix ..." << std::endl;
1390  viennacl::copy(ublas_A, vcl_slice_A);
1391  viennacl::copy(ublas_B, vcl_slice_B);
1392  viennacl::copy(ublas_C, vcl_C);
1393  if (run_test(epsilon,
1394  ublas_A, ublas_B, ublas_C,
1395  vcl_slice_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
1396  {
1397  return EXIT_FAILURE;
1398  }
1399 
1400  std::cout << "Testing A=slice, B=slice, C=range ..." << std::endl;
1401  viennacl::copy(ublas_A, vcl_slice_A);
1402  viennacl::copy(ublas_B, vcl_slice_B);
1403  viennacl::copy(ublas_C, vcl_range_C);
1404  if (run_test(epsilon,
1405  ublas_A, ublas_B, ublas_C,
1406  vcl_slice_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
1407  {
1408  return EXIT_FAILURE;
1409  }
1410 
1411  std::cout << "Testing A=slice, B=slice, C=slice ..." << std::endl;
1412  viennacl::copy(ublas_A, vcl_slice_A);
1413  viennacl::copy(ublas_B, vcl_slice_B);
1414  viennacl::copy(ublas_C, vcl_slice_C);
1415  if (run_test(epsilon,
1416  ublas_A, ublas_B, ublas_C,
1417  vcl_slice_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
1418  {
1419  return EXIT_FAILURE;
1420  }
1421 
1422 
1423  return EXIT_SUCCESS;
1424 }
1425 
1426 
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::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...
viennacl::trans
viennacl::enable_if< viennacl::is_any_sparse_matrix< M1 >::value, matrix_expression< const M1, const M1, op_trans > >::type trans(const M1 &mat)
Returns an expression template class representing a transposed matrix.
Definition: sparse_matrix_operations.hpp:374
viennacl::matrix_slice
Class for representing strided submatrices of a bigger matrix A.
Definition: forwards.h:443
viennacl::linalg::cuda::trans
void trans(matrix_expression< const matrix_base< NumericT, SizeT, DistanceT >, const matrix_base< NumericT, SizeT, DistanceT >, op_trans > const &proxy, matrix_base< NumericT > &temp_trans)
Definition: matrix_operations.hpp:94
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
matrix.hpp
Implementation of the dense matrix class.
viennacl::backend::finish
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
viennacl::identity_matrix
Represents a vector consisting of 1 at a given index and zeros otherwise. To be used as an initialize...
Definition: matrix_def.hpp:69
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
norm_1.hpp
Generic interface for the l^1-norm. See viennacl/linalg/vector_operations.hpp for implementations...
norm_frobenius.hpp
Generic interface for the Frobenius norm.
viennacl::scalar_matrix
Represents a vector consisting of scalars 's' only, i.e. v[i] = s for all i. To be used as an initial...
Definition: matrix_def.hpp:93
GENERATE_UNARY_OP_TEST
#define GENERATE_UNARY_OP_TEST(FUNCNAME)
viennacl::zero_matrix
Represents a vector consisting of zeros only. To be used as an initializer for viennacl::vector, vector_range, or vector_slize only.
Definition: matrix_def.hpp:81
viennacl::result_of::cpu_value_type::type
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:271
vector_proxy.hpp
Proxy classes for vectors.
check_for_equality
bool check_for_equality(MatrixType const &ublas_A, VCLMatrixType const &vcl_A, double epsilon)
Definition: matrix_float_double.hpp:54
matrix_proxy.hpp
Proxy classes for matrices.
viennacl::linalg::norm_frobenius
scalar_expression< const matrix_base< NumericT >, const matrix_base< NumericT >, op_norm_frobenius > norm_frobenius(const matrix_base< NumericT > &A)
Definition: norm_frobenius.hpp:61
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
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
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)
viennacl::matrix_range
Class for representing non-strided submatrices of a bigger matrix A.
Definition: forwards.h:440
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
scalar.hpp
Implementation of the ViennaCL scalar class.
run_test
int run_test(double epsilon, UBLASMatrixType &ublas_A, UBLASMatrixType &ublas_B, UBLASMatrixType &ublas_C, ViennaCLMatrixType1 &vcl_A, ViennaCLMatrixType2 &vcl_B, ViennaCLMatrixType3 vcl_C)
Definition: matrix_float_double.hpp:87
norm_inf.hpp
Generic interface for the l^infty-norm. See viennacl/linalg/vector_operations.hpp for implementations...