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
matrix_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 
23 //
24 // *** System
25 //
26 #include <iostream>
27 #include <vector>
28 
29 //
30 // *** ViennaCL
31 //
32 //#define VIENNACL_DEBUG_ALL
33 #include "viennacl/scalar.hpp"
34 #include "viennacl/matrix.hpp"
35 #include "viennacl/vector.hpp"
36 #include "viennacl/linalg/prod.hpp"
37 #include "viennacl/linalg/norm_2.hpp"
38 #include "viennacl/linalg/direct_solve.hpp"
39 #include "viennacl/linalg/lu.hpp"
40 #include "viennacl/linalg/sum.hpp"
41 #include "viennacl/tools/random.hpp"
42 
43 //
44 // -------------------------------------------------------------
45 //
46 template<typename ScalarType>
47 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
48 {
49  viennacl::backend::finish();
50  if (s1 != s2)
51  return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
52  return 0;
53 }
54 
55 template<typename ScalarType, typename VCLVectorType>
56 ScalarType diff(std::vector<ScalarType> const & v1, VCLVectorType const & v2)
57 {
58  std::vector<ScalarType> v2_cpu(v2.size());
59  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
60  viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
61 
62  ScalarType norm_inf = 0;
63  for (unsigned int i=0;i<v1.size(); ++i)
64  {
65  if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
66  {
67  ScalarType tmp = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
68  if (tmp > norm_inf)
69  norm_inf = tmp;
70  }
71  }
72 
73  return norm_inf;
74 }
75 
76 template<typename ScalarType, typename VCLMatrixType>
77 ScalarType diff(std::vector<std::vector<ScalarType> > const & mat1, VCLMatrixType const & mat2)
78 {
79  std::vector<std::vector<ScalarType> > mat2_cpu(mat2.size1(), std::vector<ScalarType>(mat2.size2()));
80  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
81  viennacl::copy(mat2, mat2_cpu);
82  ScalarType ret = 0;
83  ScalarType act = 0;
84 
85  for (std::size_t i = 0; i < mat2_cpu.size(); ++i)
86  {
87  for (std::size_t j = 0; j < mat2_cpu[i].size(); ++j)
88  {
89  act = std::fabs(mat2_cpu[i][j] - mat1[i][j]) / std::max( std::fabs(mat2_cpu[i][j]), std::fabs(mat1[i][j]) );
90  if (act > ret)
91  ret = act;
92  }
93  }
94  //std::cout << ret << std::endl;
95  return ret;
96 }
97 //
98 // -------------------------------------------------------------
99 //
100 
101 template<typename NumericT, typename Epsilon,
102  typename STLMatrixType, typename STLVectorType,
103  typename VCLMatrixType, typename VCLVectorType1, typename VCLVectorType2>
104 int test_prod_rank1(Epsilon const & epsilon,
105  STLMatrixType & std_m1, STLVectorType & std_v1, STLVectorType & std_v2, STLMatrixType & std_m2,
106  VCLMatrixType & vcl_m1, VCLVectorType1 & vcl_v1, VCLVectorType2 & vcl_v2, VCLMatrixType & vcl_m2)
107 {
108  int retval = EXIT_SUCCESS;
109 
110  // sync data:
111  std_v1 = std::vector<NumericT>(std_v1.size(), NumericT(0.1234));
112  std_v2 = std::vector<NumericT>(std_v2.size(), NumericT(0.4321));
113  viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
114  viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
115  viennacl::copy(std_m1, vcl_m1);
116 
117  // --------------------------------------------------------------------------
118  std::cout << "Rank 1 update" << std::endl;
119 
120  for (std::size_t i=0; i<std_m1.size(); ++i)
121  for (std::size_t j=0; j<std_m1[i].size(); ++j)
122  std_m1[i][j] += std_v1[i] * std_v2[j];
123  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
124  if ( std::fabs(diff(std_m1, vcl_m1)) > epsilon )
125  {
126  std::cout << "# Error at operation: rank 1 update" << std::endl;
127  std::cout << " diff: " << std::fabs(diff(std_m1, vcl_m1)) << std::endl;
128  return EXIT_FAILURE;
129  }
130 
131 
132 
133  // --------------------------------------------------------------------------
134  std::cout << "Scaled rank 1 update - CPU Scalar" << std::endl;
135  for (std::size_t i=0; i<std_m1.size(); ++i)
136  for (std::size_t j=0; j<std_m1[i].size(); ++j)
137  std_m1[i][j] += NumericT(4.2) * std_v1[i] * std_v2[j];
138  vcl_m1 += NumericT(2.1) * viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
139  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2) * NumericT(2.1); //check proper compilation
140  if ( std::fabs(diff(std_m1, vcl_m1)) > epsilon )
141  {
142  std::cout << "# Error at operation: scaled rank 1 update - CPU Scalar" << std::endl;
143  std::cout << " diff: " << std::fabs(diff(std_m1, vcl_m1)) << std::endl;
144  return EXIT_FAILURE;
145  }
146 
147  // --------------------------------------------------------------------------
148  std::cout << "Scaled rank 1 update - GPU Scalar" << std::endl;
149  for (std::size_t i=0; i<std_m1.size(); ++i)
150  for (std::size_t j=0; j<std_m1[i].size(); ++j)
151  std_m1[i][j] += NumericT(4.2) * std_v1[i] * std_v2[j];
152  vcl_m1 += viennacl::scalar<NumericT>(NumericT(2.1)) * viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
153  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2) * viennacl::scalar<NumericT>(NumericT(2.1)); //check proper compilation
154  if ( std::fabs(diff(std_m1, vcl_m1)) > epsilon )
155  {
156  std::cout << "# Error at operation: scaled rank 1 update - GPU Scalar" << std::endl;
157  std::cout << " diff: " << std::fabs(diff(std_m1, vcl_m1)) << std::endl;
158  return EXIT_FAILURE;
159  }
160 
161  //reset vcl_matrix:
162  viennacl::copy(std_m1, vcl_m1);
163 
164  // --------------------------------------------------------------------------
165  std::cout << "Matrix-Vector product" << std::endl;
166  for (std::size_t i=0; i<std_m1.size(); ++i)
167  {
168  std_v1[i] = 0;
169  for (std::size_t j=0; j<std_m1[i].size(); ++j)
170  std_v1[i] += std_m1[i][j] * std_v2[j];
171  }
172  vcl_v1 = viennacl::linalg::prod(vcl_m1, vcl_v2);
173 
174  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
175  {
176  std::cout << "# Error at operation: matrix-vector product" << std::endl;
177  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
178  retval = EXIT_FAILURE;
179  }
180  // --------------------------------------------------------------------------
181  std::cout << "Matrix-Vector product with scaled add" << std::endl;
182  NumericT alpha = static_cast<NumericT>(2.786);
183  NumericT beta = static_cast<NumericT>(1.432);
184  viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
185  viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
186 
187  for (std::size_t i=0; i<std_m1.size(); ++i)
188  {
189  NumericT tmp = 0;
190  for (std::size_t j=0; j<std_m1[i].size(); ++j)
191  tmp += std_m1[i][j] * std_v2[j];
192  std_v1[i] = alpha * tmp + beta * std_v1[i];
193  }
194  vcl_v1 = alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) + beta * vcl_v1;
195 
196  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
197  {
198  std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
199  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
200  retval = EXIT_FAILURE;
201  }
202  // --------------------------------------------------------------------------
203  std::cout << "Matrix-Vector product with matrix expression" << std::endl;
204  for (std::size_t i=0; i<std_m1.size(); ++i)
205  {
206  std_v1[i] = 0;
207  for (std::size_t j=0; j<std_m1[i].size(); ++j)
208  std_v1[i] += (std_m1[i][j] + std_m1[i][j]) * std_v2[j];
209  }
210  vcl_v1 = viennacl::linalg::prod(vcl_m1 + vcl_m1, vcl_v2);
211 
212  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
213  {
214  std::cout << "# Error at operation: matrix-vector product" << std::endl;
215  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
216  retval = EXIT_FAILURE;
217  }
218  // --------------------------------------------------------------------------
219  std::cout << "Matrix-Vector product with vector expression" << std::endl;
220  for (std::size_t i=0; i<std_m1.size(); ++i)
221  {
222  std_v1[i] = 0;
223  for (std::size_t j=0; j<std_m1[i].size(); ++j)
224  std_v1[i] += std_m1[i][j] * NumericT(3) * std_v2[j];
225  }
226  vcl_v1 = viennacl::linalg::prod(vcl_m1, NumericT(3) * vcl_v2);
227 
228  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
229  {
230  std::cout << "# Error at operation: matrix-vector product" << std::endl;
231  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
232  retval = EXIT_FAILURE;
233  }
234  // --------------------------------------------------------------------------
235  std::cout << "Matrix-Vector product with matrix and vector expression" << std::endl;
236  for (std::size_t i=0; i<std_m1.size(); ++i)
237  {
238  std_v1[i] = 0;
239  for (std::size_t j=0; j<std_m1[i].size(); ++j)
240  std_v1[i] += (std_m1[i][j] + std_m1[i][j]) * (std_v2[j] + std_v2[j]);
241  }
242  vcl_v1 = viennacl::linalg::prod(vcl_m1 + vcl_m1, vcl_v2 + vcl_v2);
243 
244  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
245  {
246  std::cout << "# Error at operation: matrix-vector product" << std::endl;
247  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
248  retval = EXIT_FAILURE;
249  }
250  // --------------------------------------------------------------------------
251 
252  viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
253  viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
254 
255  std::cout << "Transposed Matrix-Vector product" << std::endl;
256  for (std::size_t i=0; i<std_m1[0].size(); ++i)
257  {
258  std_v2[i] = 0;
259  for (std::size_t j=0; j<std_m1.size(); ++j)
260  std_v2[i] += alpha * std_m1[j][i] * std_v1[j];
261  }
262  vcl_v2 = alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1);
263 
264  if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
265  {
266  std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
267  std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
268  retval = EXIT_FAILURE;
269  }
270 
271  std::cout << "Transposed Matrix-Vector product with scaled add" << std::endl;
272  for (std::size_t i=0; i<std_m1[0].size(); ++i)
273  {
274  NumericT tmp = 0;
275  for (std::size_t j=0; j<std_m1.size(); ++j)
276  tmp += std_m1[j][i] * std_v1[j];
277  std_v2[i] = alpha * tmp + beta * std_v2[i];
278  }
279  vcl_v2 = alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2;
280 
281  if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
282  {
283  std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
284  std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
285  retval = EXIT_FAILURE;
286  }
287  // --------------------------------------------------------------------------
288 
289  std::cout << "Row sum with matrix" << std::endl;
290  for (std::size_t i=0; i<std_m1.size(); ++i)
291  {
292  std_v1[i] = 0;
293  for (std::size_t j=0; j<std_m1[i].size(); ++j)
294  std_v1[i] += std_m1[i][j];
295  }
296  vcl_v1 = viennacl::linalg::row_sum(vcl_m1);
297 
298  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
299  {
300  std::cout << "# Error at operation: row sum" << std::endl;
301  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
302  retval = EXIT_FAILURE;
303  }
304  // --------------------------------------------------------------------------
305 
306  std::cout << "Row sum with matrix expression" << std::endl;
307  for (std::size_t i=0; i<std_m1.size(); ++i)
308  {
309  std_v1[i] = 0;
310  for (std::size_t j=0; j<std_m1[i].size(); ++j)
311  std_v1[i] += std_m1[i][j] + std_m1[i][j];
312  }
313  vcl_v1 = viennacl::linalg::row_sum(vcl_m1 + vcl_m1);
314 
315  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
316  {
317  std::cout << "# Error at operation: row sum (with expression)" << std::endl;
318  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
319  retval = EXIT_FAILURE;
320  }
321  // --------------------------------------------------------------------------
322 
323  std::cout << "Column sum with matrix" << std::endl;
324  for (std::size_t i=0; i<std_m1[0].size(); ++i)
325  {
326  std_v2[i] = 0;
327  for (std::size_t j=0; j<std_m1.size(); ++j)
328  std_v2[i] += std_m1[j][i];
329  }
330  vcl_v2 = viennacl::linalg::column_sum(vcl_m1);
331 
332  if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
333  {
334  std::cout << "# Error at operation: column sum" << std::endl;
335  std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
336  retval = EXIT_FAILURE;
337  }
338  // --------------------------------------------------------------------------
339 
340  std::cout << "Column sum with matrix expression" << std::endl;
341  for (std::size_t i=0; i<std_m1[0].size(); ++i)
342  {
343  std_v2[i] = 0;
344  for (std::size_t j=0; j<std_m1.size(); ++j)
345  std_v2[i] += std_m1[j][i] + std_m1[j][i];
346  }
347  vcl_v2 = viennacl::linalg::column_sum(vcl_m1 + vcl_m1);
348 
349  if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
350  {
351  std::cout << "# Error at operation: column sum (with expression)" << std::endl;
352  std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
353  retval = EXIT_FAILURE;
354  }
355  // --------------------------------------------------------------------------
356 
357 
358  viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
359  viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
360 
361  std::cout << "Row extraction from matrix" << std::endl;
362  for (std::size_t j=0; j<std_m1[7].size(); ++j)
363  std_v2[j] = std_m1[7][j];
364  vcl_v2 = row(vcl_m1, std::size_t(7));
365 
366  if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
367  {
368  std::cout << "# Error at operation: diagonal extraction from matrix" << std::endl;
369  std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
370  retval = EXIT_FAILURE;
371  }
372 
373  std::cout << "Column extraction from matrix" << std::endl;
374  for (std::size_t i=0; i<std_m1.size(); ++i)
375  std_v1[i] = std_m1[i][7];
376  vcl_v1 = column(vcl_m1, std::size_t(7));
377 
378  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
379  {
380  std::cout << "# Error at operation: diagonal extraction from matrix" << std::endl;
381  std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
382  retval = EXIT_FAILURE;
383  }
384 
385  // --------------------------------------------------------------------------
386 
387  viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
388  viennacl::copy(std_v2.begin(), std_v2.end(), vcl_v2.begin());
389  viennacl::copy(std_m2, vcl_m2);
390  STLMatrixType A = std_m2;
391 
392  std::cout << "Diagonal extraction from matrix" << std::endl;
393  for (std::size_t i=0; i<std_m1[0].size(); ++i)
394  std_v2[i] = std_m1[i + 3][i];
395  vcl_v2 = diag(vcl_m1, static_cast<int>(-3));
396 
397  if ( std::fabs(diff(std_v2, vcl_v2)) > epsilon )
398  {
399  std::cout << "# Error at operation: diagonal extraction from matrix" << std::endl;
400  std::cout << " diff: " << std::fabs(diff(std_v2, vcl_v2)) << std::endl;
401  retval = EXIT_FAILURE;
402  }
403 
404  std::cout << "Matrix diagonal assignment from vector" << std::endl;
405  A = std::vector<std::vector<NumericT> >(A.size(), std::vector<NumericT>(A[0].size()));
406  for (std::size_t i=0; i<std_m1[0].size(); ++i)
407  A[i + (A.size() - std_m1[i].size())][i] = std_v2[i];
408  vcl_m2 = diag(vcl_v2, static_cast<int>(std_m1[0].size()) - static_cast<int>(A.size()));
409 
410  if ( std::fabs(diff(A, vcl_m2)) > epsilon )
411  {
412  std::cout << "# Error at operation: Matrix assignment from diagonal" << std::endl;
413  std::cout << " diff: " << std::fabs(diff(A, vcl_m2)) << std::endl;
414  retval = EXIT_FAILURE;
415  }
416  // --------------------------------------------------------------------------
417 
418  return retval;
419 }
420 
421 template<typename NumericT>
422 void inplace_solve_upper(std::vector<std::vector<NumericT> > const & A, std::vector<NumericT> & b, bool unit_diagonal)
423 {
424  for (std::size_t i2=0; i2<A.size(); ++i2)
425  {
426  std::size_t i = A.size() - i2 - 1;
427  for (std::size_t j = i+1; j < A.size(); ++j)
428  b[i] -= A[i][j] * b[j];
429  b[i] = unit_diagonal ? b[i] : b[i] / A[i][i];
430  }
431 }
432 
433 template<typename NumericT>
434 void inplace_solve(std::vector<std::vector<NumericT> > const & A, std::vector<NumericT> & b, viennacl::linalg::upper_tag)
435 {
436  inplace_solve_upper(A, b, false);
437 }
438 
439 template<typename NumericT>
440 void inplace_solve(std::vector<std::vector<NumericT> > const & A, std::vector<NumericT> & b, viennacl::linalg::unit_upper_tag)
441 {
442  inplace_solve_upper(A, b, true);
443 }
444 
445 
446 template<typename NumericT>
447 void inplace_solve_lower(std::vector<std::vector<NumericT> > const & A, std::vector<NumericT> & b, bool unit_diagonal)
448 {
449  for (std::size_t i=0; i<A.size(); ++i)
450  {
451  for (std::size_t j = 0; j < i; ++j)
452  b[i] -= A[i][j] * b[j];
453  b[i] = unit_diagonal ? b[i] : b[i] / A[i][i];
454  }
455 }
456 
457 template<typename NumericT>
458 void inplace_solve(std::vector<std::vector<NumericT> > const & A, std::vector<NumericT> & b, viennacl::linalg::lower_tag)
459 {
460  inplace_solve_lower(A, b, false);
461 }
462 
463 template<typename NumericT>
464 void inplace_solve(std::vector<std::vector<NumericT> > const & A, std::vector<NumericT> & b, viennacl::linalg::unit_lower_tag)
465 {
466  inplace_solve_lower(A, b, true);
467 }
468 
469 
470 template<typename NumericT, typename TagT>
471 std::vector<NumericT> solve(std::vector<std::vector<NumericT> > const & A, std::vector<NumericT> const & b, TagT)
472 {
473  std::vector<NumericT> ret(b);
474  inplace_solve(A, ret, TagT());
475  return ret;
476 }
477 
478 template<typename NumericT, typename Epsilon,
479  typename STLMatrixType, typename STLVectorType,
480  typename VCLMatrixType, typename VCLVectorType1>
481 int test_solve(Epsilon const & epsilon,
482  STLMatrixType & std_m1, STLVectorType & std_v1,
483  VCLMatrixType & vcl_m1, VCLVectorType1 & vcl_v1)
484 {
485  int retval = EXIT_SUCCESS;
486 
487  // sync data:
488  //viennacl::copy(std_v1.begin(), std_v1.end(), vcl_v1.begin());
489  viennacl::copy(std_v1, vcl_v1);
490  viennacl::copy(std_m1, vcl_m1);
491 
493 
494  //upper triangular:
495  std::cout << "Upper triangular solver" << std::endl;
496  std_v1 = solve(std_m1, std_v1, viennacl::linalg::upper_tag());
497  vcl_v1 = viennacl::linalg::solve(vcl_m1, vcl_v1, viennacl::linalg::upper_tag());
498  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
499  {
500  std::cout << "# Error at operation: upper triangular solver" << std::endl;
501  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
502  retval = EXIT_FAILURE;
503  }
504 
505  //upper unit triangular:
506  std::cout << "Upper unit triangular solver" << std::endl;
507  viennacl::copy(std_v1, vcl_v1);
508  std_v1 = solve(std_m1, std_v1, viennacl::linalg::unit_upper_tag());
509  vcl_v1 = viennacl::linalg::solve(vcl_m1, vcl_v1, viennacl::linalg::unit_upper_tag());
510  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
511  {
512  std::cout << "# Error at operation: unit upper triangular solver" << std::endl;
513  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
514  retval = EXIT_FAILURE;
515  }
516 
517  //lower triangular:
518  std::cout << "Lower triangular solver" << std::endl;
519  viennacl::copy(std_v1, vcl_v1);
520  std_v1 = solve(std_m1, std_v1, viennacl::linalg::lower_tag());
521  vcl_v1 = viennacl::linalg::solve(vcl_m1, vcl_v1, viennacl::linalg::lower_tag());
522  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
523  {
524  std::cout << "# Error at operation: lower triangular solver" << std::endl;
525  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
526  retval = EXIT_FAILURE;
527  }
528 
529  //lower unit triangular:
530  std::cout << "Lower unit triangular solver" << std::endl;
531  viennacl::copy(std_v1, vcl_v1);
532  std_v1 = solve(std_m1, std_v1, viennacl::linalg::unit_lower_tag());
533  vcl_v1 = viennacl::linalg::solve(vcl_m1, vcl_v1, viennacl::linalg::unit_lower_tag());
534  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
535  {
536  std::cout << "# Error at operation: unit lower triangular solver" << std::endl;
537  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
538  retval = EXIT_FAILURE;
539  }
540 
541 
542 
543 
544  STLMatrixType std_m1_trans(std_m1[0].size(), std::vector<NumericT>(std_m1.size()));
545  for (std::size_t i=0; i<std_m1.size(); ++i)
546  for (std::size_t j=0; j<std_m1[i].size(); ++j)
547  std_m1_trans[j][i] = std_m1[i][j];
548 
549 
550  //transposed upper triangular:
551  std::cout << "Transposed upper triangular solver" << std::endl;
552  viennacl::copy(std_v1, vcl_v1);
553  std_v1 = solve(std_m1_trans, std_v1, viennacl::linalg::upper_tag());
554  vcl_v1 = viennacl::linalg::solve(trans(vcl_m1), vcl_v1, viennacl::linalg::upper_tag());
555  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
556  {
557  std::cout << "# Error at operation: upper triangular solver" << std::endl;
558  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
559  retval = EXIT_FAILURE;
560  }
561 
562  //transposed upper unit triangular:
563  std::cout << "Transposed unit upper triangular solver" << std::endl;
564  viennacl::copy(std_v1, vcl_v1);
565  std_v1 = solve(std_m1_trans, std_v1, viennacl::linalg::unit_upper_tag());
566  vcl_v1 = viennacl::linalg::solve(trans(vcl_m1), vcl_v1, viennacl::linalg::unit_upper_tag());
567  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
568  {
569  std::cout << "# Error at operation: unit upper triangular solver" << std::endl;
570  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
571  retval = EXIT_FAILURE;
572  }
573 
574  //transposed lower triangular:
575  std::cout << "Transposed lower triangular solver" << std::endl;
576  viennacl::copy(std_v1, vcl_v1);
577  std_v1 = solve(std_m1_trans, std_v1, viennacl::linalg::lower_tag());
578  vcl_v1 = viennacl::linalg::solve(trans(vcl_m1), vcl_v1, viennacl::linalg::lower_tag());
579  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
580  {
581  std::cout << "# Error at operation: lower triangular solver" << std::endl;
582  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
583  retval = EXIT_FAILURE;
584  }
585 
586  //transposed lower unit triangular:
587  std::cout << "Transposed unit lower triangular solver" << std::endl;
588  viennacl::copy(std_v1, vcl_v1);
589  std_v1 = solve(std_m1_trans, std_v1, viennacl::linalg::unit_lower_tag());
590  vcl_v1 = viennacl::linalg::solve(trans(vcl_m1), vcl_v1, viennacl::linalg::unit_lower_tag());
591  if ( std::fabs(diff(std_v1, vcl_v1)) > epsilon )
592  {
593  std::cout << "# Error at operation: unit lower triangular solver" << std::endl;
594  std::cout << " diff: " << std::fabs(diff(std_v1, vcl_v1)) << std::endl;
595  retval = EXIT_FAILURE;
596  }
597 
598  return retval;
599 }
600 
601 
602 //
603 // -------------------------------------------------------------
604 //
605 template< typename NumericT, typename F, typename Epsilon >
606 int test(Epsilon const& epsilon)
607 {
608  int retval = EXIT_SUCCESS;
609 
610  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
611 
612  std::size_t num_rows = 141; //note: use num_rows > num_cols + 3 for diag() tests to work
613  std::size_t num_cols = 103;
614 
615  // --------------------------------------------------------------------------
616  std::vector<NumericT> std_v1(num_rows);
617  for (std::size_t i = 0; i < std_v1.size(); ++i)
618  std_v1[i] = randomNumber();
619  std::vector<NumericT> std_v2 = std::vector<NumericT>(num_cols, NumericT(3.1415));
620 
621 
622  std::vector<std::vector<NumericT> > std_m1(std_v1.size(), std::vector<NumericT>(std_v2.size()));
623 
624  for (std::size_t i = 0; i < std_m1.size(); ++i)
625  for (std::size_t j = 0; j < std_m1[i].size(); ++j)
626  std_m1[i][j] = static_cast<NumericT>(0.1) * randomNumber();
627 
628 
629  std::vector<std::vector<NumericT> > std_m2(std_v1.size(), std::vector<NumericT>(std_v1.size()));
630 
631  for (std::size_t i = 0; i < std_m2.size(); ++i)
632  {
633  for (std::size_t j = 0; j < std_m2[i].size(); ++j)
634  std_m2[i][j] = static_cast<NumericT>(-0.1) * randomNumber();
635  std_m2[i][i] = static_cast<NumericT>(2) + randomNumber();
636  }
637 
638 
639  viennacl::vector<NumericT> vcl_v1_native(std_v1.size());
640  viennacl::vector<NumericT> vcl_v1_large(4 * std_v1.size());
641  viennacl::vector_range< viennacl::vector<NumericT> > vcl_v1_range(vcl_v1_large, viennacl::range(3, std_v1.size() + 3));
642  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v1_slice(vcl_v1_large, viennacl::slice(2, 3, std_v1.size()));
643 
644  viennacl::vector<NumericT> vcl_v2_native(std_v2.size());
645  viennacl::vector<NumericT> vcl_v2_large(4 * std_v2.size());
646  viennacl::vector_range< viennacl::vector<NumericT> > vcl_v2_range(vcl_v2_large, viennacl::range(8, std_v2.size() + 8));
647  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v2_slice(vcl_v2_large, viennacl::slice(6, 2, std_v2.size()));
648 
649  viennacl::matrix<NumericT, F> vcl_m1_native(std_m1.size(), std_m1[0].size());
650  viennacl::matrix<NumericT, F> vcl_m1_large(4 * std_m1.size(), 4 * std_m1[0].size());
651  viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m1_range(vcl_m1_large,
652  viennacl::range(8, std_m1.size() + 8),
653  viennacl::range(std_m1[0].size(), 2 * std_m1[0].size()) );
654  viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m1_slice(vcl_m1_large,
655  viennacl::slice(6, 2, std_m1.size()),
656  viennacl::slice(std_m1[0].size(), 2, std_m1[0].size()) );
657 
658  viennacl::matrix<NumericT, F> vcl_m2_native(std_m2.size(), std_m2[0].size());
659  viennacl::matrix<NumericT, F> vcl_m2_large(4 * std_m2.size(), 4 * std_m2[0].size());
660  viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m2_range(vcl_m2_large,
661  viennacl::range(8, std_m2.size() + 8),
662  viennacl::range(std_m2[0].size(), 2 * std_m2[0].size()) );
663  viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m2_slice(vcl_m2_large,
664  viennacl::slice(6, 2, std_m2.size()),
665  viennacl::slice(std_m2[0].size(), 2, std_m2[0].size()) );
666 
667 
668 /* std::cout << "Matrix resizing (to larger)" << std::endl;
669  matrix.resize(2*num_rows, 2*num_cols, true);
670  for (unsigned int i = 0; i < matrix.size1(); ++i)
671  {
672  for (unsigned int j = (i<result.size() ? rhs.size() : 0); j < matrix.size2(); ++j)
673  matrix(i,j) = 0;
674  }
675  vcl_matrix.resize(2*num_rows, 2*num_cols, true);
676  viennacl::copy(vcl_matrix, matrix);
677  if ( std::fabs(diff(matrix, vcl_matrix)) > epsilon )
678  {
679  std::cout << "# Error at operation: matrix resize (to larger)" << std::endl;
680  std::cout << " diff: " << std::fabs(diff(matrix, vcl_matrix)) << std::endl;
681  return EXIT_FAILURE;
682  }
683 
684  matrix(12, 14) = NumericT(1.9);
685  matrix(19, 16) = NumericT(1.0);
686  matrix (13, 15) = NumericT(-9);
687  vcl_matrix(12, 14) = NumericT(1.9);
688  vcl_matrix(19, 16) = NumericT(1.0);
689  vcl_matrix (13, 15) = NumericT(-9);
690 
691  std::cout << "Matrix resizing (to smaller)" << std::endl;
692  matrix.resize(result.size(), rhs.size(), true);
693  vcl_matrix.resize(result.size(), rhs.size(), true);
694  if ( std::fabs(diff(matrix, vcl_matrix)) > epsilon )
695  {
696  std::cout << "# Error at operation: matrix resize (to smaller)" << std::endl;
697  std::cout << " diff: " << std::fabs(diff(matrix, vcl_matrix)) << std::endl;
698  return EXIT_FAILURE;
699  }
700  */
701 
702  //
703  // Run a bunch of tests for rank-1-updates, matrix-vector products
704  //
705  std::cout << "------------ Testing rank-1-updates and matrix-vector products ------------------" << std::endl;
706 
707  std::cout << "* m = full, v1 = full, v2 = full" << std::endl;
708  retval = test_prod_rank1<NumericT>(epsilon,
709  std_m1, std_v1, std_v2, std_m2,
710  vcl_m1_native, vcl_v1_native, vcl_v2_native, vcl_m2_native);
711  if (retval == EXIT_FAILURE)
712  {
713  std::cout << " --- FAILED! ---" << std::endl;
714  return retval;
715  }
716  else
717  std::cout << " --- PASSED ---" << std::endl;
718 
719 
720  std::cout << "* m = full, v1 = full, v2 = range" << std::endl;
721  retval = test_prod_rank1<NumericT>(epsilon,
722  std_m1, std_v1, std_v2, std_m2,
723  vcl_m1_native, vcl_v1_native, vcl_v2_range, vcl_m2_native);
724  if (retval == EXIT_FAILURE)
725  {
726  std::cout << " --- FAILED! ---" << std::endl;
727  return retval;
728  }
729  else
730  std::cout << " --- PASSED ---" << std::endl;
731 
732 
733  std::cout << "* m = full, v1 = full, v2 = slice" << std::endl;
734  retval = test_prod_rank1<NumericT>(epsilon,
735  std_m1, std_v1, std_v2, std_m2,
736  vcl_m1_native, vcl_v1_native, vcl_v2_slice, vcl_m2_native);
737  if (retval == EXIT_FAILURE)
738  {
739  std::cout << " --- FAILED! ---" << std::endl;
740  return retval;
741  }
742  else
743  std::cout << " --- PASSED ---" << std::endl;
744 
745 
746  // v1 = range
747 
748 
749  std::cout << "* m = full, v1 = range, v2 = full" << std::endl;
750  retval = test_prod_rank1<NumericT>(epsilon,
751  std_m1, std_v1, std_v2, std_m2,
752  vcl_m1_native, vcl_v1_range, vcl_v2_native, vcl_m2_native);
753  if (retval == EXIT_FAILURE)
754  {
755  std::cout << " --- FAILED! ---" << std::endl;
756  return retval;
757  }
758  else
759  std::cout << " --- PASSED ---" << std::endl;
760 
761 
762  std::cout << "* m = full, v1 = range, v2 = range" << std::endl;
763  retval = test_prod_rank1<NumericT>(epsilon,
764  std_m1, std_v1, std_v2, std_m2,
765  vcl_m1_native, vcl_v1_range, vcl_v2_range, vcl_m2_native);
766  if (retval == EXIT_FAILURE)
767  {
768  std::cout << " --- FAILED! ---" << std::endl;
769  return retval;
770  }
771  else
772  std::cout << " --- PASSED ---" << std::endl;
773 
774 
775  std::cout << "* m = full, v1 = range, v2 = slice" << std::endl;
776  retval = test_prod_rank1<NumericT>(epsilon,
777  std_m1, std_v1, std_v2, std_m2,
778  vcl_m1_native, vcl_v1_range, vcl_v2_slice, vcl_m2_native);
779  if (retval == EXIT_FAILURE)
780  {
781  std::cout << " --- FAILED! ---" << std::endl;
782  return retval;
783  }
784  else
785  std::cout << " --- PASSED ---" << std::endl;
786 
787 
788 
789  // v1 = slice
790 
791  std::cout << "* m = full, v1 = slice, v2 = full" << std::endl;
792  retval = test_prod_rank1<NumericT>(epsilon,
793  std_m1, std_v1, std_v2, std_m2,
794  vcl_m1_native, vcl_v1_slice, vcl_v2_native, vcl_m2_native);
795  if (retval == EXIT_FAILURE)
796  {
797  std::cout << " --- FAILED! ---" << std::endl;
798  return retval;
799  }
800  else
801  std::cout << " --- PASSED ---" << std::endl;
802 
803 
804  std::cout << "* m = full, v1 = slice, v2 = range" << std::endl;
805  retval = test_prod_rank1<NumericT>(epsilon,
806  std_m1, std_v1, std_v2, std_m2,
807  vcl_m1_native, vcl_v1_slice, vcl_v2_range, vcl_m2_native);
808  if (retval == EXIT_FAILURE)
809  {
810  std::cout << " --- FAILED! ---" << std::endl;
811  return retval;
812  }
813  else
814  std::cout << " --- PASSED ---" << std::endl;
815 
816 
817  std::cout << "* m = full, v1 = slice, v2 = slice" << std::endl;
818  retval = test_prod_rank1<NumericT>(epsilon,
819  std_m1, std_v1, std_v2, std_m2,
820  vcl_m1_native, vcl_v1_slice, vcl_v2_slice, vcl_m2_native);
821  if (retval == EXIT_FAILURE)
822  {
823  std::cout << " --- FAILED! ---" << std::endl;
824  return retval;
825  }
826  else
827  std::cout << " --- PASSED ---" << std::endl;
828 
829 
831 
832  std::cout << "* m = range, v1 = full, v2 = full" << std::endl;
833  retval = test_prod_rank1<NumericT>(epsilon,
834  std_m1, std_v1, std_v2, std_m2,
835  vcl_m1_range, vcl_v1_native, vcl_v2_native, vcl_m2_range);
836  if (retval == EXIT_FAILURE)
837  {
838  std::cout << " --- FAILED! ---" << std::endl;
839  return retval;
840  }
841  else
842  std::cout << " --- PASSED ---" << std::endl;
843 
844 
845  std::cout << "* m = range, v1 = full, v2 = range" << std::endl;
846  retval = test_prod_rank1<NumericT>(epsilon,
847  std_m1, std_v1, std_v2, std_m2,
848  vcl_m1_range, vcl_v1_native, vcl_v2_range, vcl_m2_range);
849  if (retval == EXIT_FAILURE)
850  {
851  std::cout << " --- FAILED! ---" << std::endl;
852  return retval;
853  }
854  else
855  std::cout << " --- PASSED ---" << std::endl;
856 
857 
858  std::cout << "* m = range, v1 = full, v2 = slice" << std::endl;
859  retval = test_prod_rank1<NumericT>(epsilon,
860  std_m1, std_v1, std_v2, std_m2,
861  vcl_m1_range, vcl_v1_native, vcl_v2_slice, vcl_m2_range);
862  if (retval == EXIT_FAILURE)
863  {
864  std::cout << " --- FAILED! ---" << std::endl;
865  return retval;
866  }
867  else
868  std::cout << " --- PASSED ---" << std::endl;
869 
870 
871  // v1 = range
872 
873 
874  std::cout << "* m = range, v1 = range, v2 = full" << std::endl;
875  retval = test_prod_rank1<NumericT>(epsilon,
876  std_m1, std_v1, std_v2, std_m2,
877  vcl_m1_range, vcl_v1_range, vcl_v2_native, vcl_m2_range);
878  if (retval == EXIT_FAILURE)
879  {
880  std::cout << " --- FAILED! ---" << std::endl;
881  return retval;
882  }
883  else
884  std::cout << " --- PASSED ---" << std::endl;
885 
886 
887  std::cout << "* m = range, v1 = range, v2 = range" << std::endl;
888  retval = test_prod_rank1<NumericT>(epsilon,
889  std_m1, std_v1, std_v2, std_m2,
890  vcl_m1_range, vcl_v1_range, vcl_v2_range, vcl_m2_range);
891  if (retval == EXIT_FAILURE)
892  {
893  std::cout << " --- FAILED! ---" << std::endl;
894  return retval;
895  }
896  else
897  std::cout << " --- PASSED ---" << std::endl;
898 
899 
900  std::cout << "* m = range, v1 = range, v2 = slice" << std::endl;
901  retval = test_prod_rank1<NumericT>(epsilon,
902  std_m1, std_v1, std_v2, std_m2,
903  vcl_m1_range, vcl_v1_range, vcl_v2_slice, vcl_m2_range);
904  if (retval == EXIT_FAILURE)
905  {
906  std::cout << " --- FAILED! ---" << std::endl;
907  return retval;
908  }
909  else
910  std::cout << " --- PASSED ---" << std::endl;
911 
912 
913 
914  // v1 = slice
915 
916  std::cout << "* m = range, v1 = slice, v2 = full" << std::endl;
917  retval = test_prod_rank1<NumericT>(epsilon,
918  std_m1, std_v1, std_v2, std_m2,
919  vcl_m1_range, vcl_v1_slice, vcl_v2_native, vcl_m2_range);
920  if (retval == EXIT_FAILURE)
921  {
922  std::cout << " --- FAILED! ---" << std::endl;
923  return retval;
924  }
925  else
926  std::cout << " --- PASSED ---" << std::endl;
927 
928 
929  std::cout << "* m = range, v1 = slice, v2 = range" << std::endl;
930  retval = test_prod_rank1<NumericT>(epsilon,
931  std_m1, std_v1, std_v2, std_m2,
932  vcl_m1_range, vcl_v1_slice, vcl_v2_range, vcl_m2_range);
933  if (retval == EXIT_FAILURE)
934  {
935  std::cout << " --- FAILED! ---" << std::endl;
936  return retval;
937  }
938  else
939  std::cout << " --- PASSED ---" << std::endl;
940 
941 
942  std::cout << "* m = range, v1 = slice, v2 = slice" << std::endl;
943  retval = test_prod_rank1<NumericT>(epsilon,
944  std_m1, std_v1, std_v2, std_m2,
945  vcl_m1_range, vcl_v1_slice, vcl_v2_slice, vcl_m2_range);
946  if (retval == EXIT_FAILURE)
947  {
948  std::cout << " --- FAILED! ---" << std::endl;
949  return retval;
950  }
951  else
952  std::cout << " --- PASSED ---" << std::endl;
953 
954 
956 
957  std::cout << "* m = slice, v1 = full, v2 = full" << std::endl;
958  retval = test_prod_rank1<NumericT>(epsilon,
959  std_m1, std_v1, std_v2, std_m2,
960  vcl_m1_slice, vcl_v1_native, vcl_v2_native, vcl_m2_slice);
961  if (retval == EXIT_FAILURE)
962  {
963  std::cout << " --- FAILED! ---" << std::endl;
964  return retval;
965  }
966  else
967  std::cout << " --- PASSED ---" << std::endl;
968 
969 
970  std::cout << "* m = slice, v1 = full, v2 = range" << std::endl;
971  retval = test_prod_rank1<NumericT>(epsilon,
972  std_m1, std_v1, std_v2, std_m2,
973  vcl_m1_slice, vcl_v1_native, vcl_v2_range, vcl_m2_slice);
974  if (retval == EXIT_FAILURE)
975  {
976  std::cout << " --- FAILED! ---" << std::endl;
977  return retval;
978  }
979  else
980  std::cout << " --- PASSED ---" << std::endl;
981 
982 
983  std::cout << "* m = slice, v1 = full, v2 = slice" << std::endl;
984  retval = test_prod_rank1<NumericT>(epsilon,
985  std_m1, std_v1, std_v2, std_m2,
986  vcl_m1_slice, vcl_v1_native, vcl_v2_slice, vcl_m2_slice);
987  if (retval == EXIT_FAILURE)
988  {
989  std::cout << " --- FAILED! ---" << std::endl;
990  return retval;
991  }
992  else
993  std::cout << " --- PASSED ---" << std::endl;
994 
995 
996  // v1 = range
997 
998 
999  std::cout << "* m = slice, v1 = range, v2 = full" << std::endl;
1000  retval = test_prod_rank1<NumericT>(epsilon,
1001  std_m1, std_v1, std_v2, std_m2,
1002  vcl_m1_slice, vcl_v1_range, vcl_v2_native, vcl_m2_slice);
1003  if (retval == EXIT_FAILURE)
1004  {
1005  std::cout << " --- FAILED! ---" << std::endl;
1006  return retval;
1007  }
1008  else
1009  std::cout << " --- PASSED ---" << std::endl;
1010 
1011 
1012  std::cout << "* m = slice, v1 = range, v2 = range" << std::endl;
1013  retval = test_prod_rank1<NumericT>(epsilon,
1014  std_m1, std_v1, std_v2, std_m2,
1015  vcl_m1_slice, vcl_v1_range, vcl_v2_range, vcl_m2_slice);
1016  if (retval == EXIT_FAILURE)
1017  {
1018  std::cout << " --- FAILED! ---" << std::endl;
1019  return retval;
1020  }
1021  else
1022  std::cout << " --- PASSED ---" << std::endl;
1023 
1024 
1025  std::cout << "* m = slice, v1 = range, v2 = slice" << std::endl;
1026  retval = test_prod_rank1<NumericT>(epsilon,
1027  std_m1, std_v1, std_v2, std_m2,
1028  vcl_m1_slice, vcl_v1_range, vcl_v2_slice, vcl_m2_slice);
1029  if (retval == EXIT_FAILURE)
1030  {
1031  std::cout << " --- FAILED! ---" << std::endl;
1032  return retval;
1033  }
1034  else
1035  std::cout << " --- PASSED ---" << std::endl;
1036 
1037 
1038 
1039  // v1 = slice
1040 
1041  std::cout << "* m = slice, v1 = slice, v2 = full" << std::endl;
1042  retval = test_prod_rank1<NumericT>(epsilon,
1043  std_m1, std_v1, std_v2, std_m2,
1044  vcl_m1_slice, vcl_v1_slice, vcl_v2_native, vcl_m2_slice);
1045  if (retval == EXIT_FAILURE)
1046  {
1047  std::cout << " --- FAILED! ---" << std::endl;
1048  return retval;
1049  }
1050  else
1051  std::cout << " --- PASSED ---" << std::endl;
1052 
1053 
1054  std::cout << "* m = slice, v1 = slice, v2 = range" << std::endl;
1055  retval = test_prod_rank1<NumericT>(epsilon,
1056  std_m1, std_v1, std_v2, std_m2,
1057  vcl_m1_slice, vcl_v1_slice, vcl_v2_range, vcl_m2_slice);
1058  if (retval == EXIT_FAILURE)
1059  {
1060  std::cout << " --- FAILED! ---" << std::endl;
1061  return retval;
1062  }
1063  else
1064  std::cout << " --- PASSED ---" << std::endl;
1065 
1066 
1067  std::cout << "* m = slice, v1 = slice, v2 = slice" << std::endl;
1068  retval = test_prod_rank1<NumericT>(epsilon,
1069  std_m1, std_v1, std_v2, std_m2,
1070  vcl_m1_slice, vcl_v1_slice, vcl_v2_slice, vcl_m2_slice);
1071  if (retval == EXIT_FAILURE)
1072  {
1073  std::cout << " --- FAILED! ---" << std::endl;
1074  return retval;
1075  }
1076  else
1077  std::cout << " --- PASSED ---" << std::endl;
1078 
1079 
1080 
1081  //
1082  // Testing triangular solve() routines
1083  //
1084 
1085  std::cout << "------------ Testing triangular solves ------------------" << std::endl;
1086 
1087  std::cout << "* m = full, v1 = full" << std::endl;
1088  retval = test_solve<NumericT>(epsilon,
1089  std_m2, std_v1,
1090  vcl_m2_native, vcl_v1_native);
1091  if (retval == EXIT_FAILURE)
1092  {
1093  std::cout << " --- FAILED! ---" << std::endl;
1094  return retval;
1095  }
1096  else
1097  std::cout << " --- PASSED ---" << std::endl;
1098 
1099  std::cout << "* m = full, v1 = range" << std::endl;
1100  retval = test_solve<NumericT>(epsilon,
1101  std_m2, std_v1,
1102  vcl_m2_native, vcl_v1_range);
1103  if (retval == EXIT_FAILURE)
1104  {
1105  std::cout << " --- FAILED! ---" << std::endl;
1106  return retval;
1107  }
1108  else
1109  std::cout << " --- PASSED ---" << std::endl;
1110 
1111  std::cout << "* m = full, v1 = slice" << std::endl;
1112  retval = test_solve<NumericT>(epsilon,
1113  std_m2, std_v1,
1114  vcl_m2_native, vcl_v1_slice);
1115  if (retval == EXIT_FAILURE)
1116  {
1117  std::cout << " --- FAILED! ---" << std::endl;
1118  return retval;
1119  }
1120  else
1121  std::cout << " --- PASSED ---" << std::endl;
1122 
1124 
1125 
1126  std::cout << "* m = range, v1 = full" << std::endl;
1127  retval = test_solve<NumericT>(epsilon,
1128  std_m2, std_v1,
1129  vcl_m2_range, vcl_v1_native);
1130  if (retval == EXIT_FAILURE)
1131  {
1132  std::cout << " --- FAILED! ---" << std::endl;
1133  return retval;
1134  }
1135  else
1136  std::cout << " --- PASSED ---" << std::endl;
1137 
1138  std::cout << "* m = range, v1 = range" << std::endl;
1139  retval = test_solve<NumericT>(epsilon,
1140  std_m2, std_v1,
1141  vcl_m2_range, vcl_v1_range);
1142  if (retval == EXIT_FAILURE)
1143  {
1144  std::cout << " --- FAILED! ---" << std::endl;
1145  return retval;
1146  }
1147  else
1148  std::cout << " --- PASSED ---" << std::endl;
1149 
1150  std::cout << "* m = range, v1 = slice" << std::endl;
1151  retval = test_solve<NumericT>(epsilon,
1152  std_m2, std_v1,
1153  vcl_m2_range, vcl_v1_slice);
1154  if (retval == EXIT_FAILURE)
1155  {
1156  std::cout << " --- FAILED! ---" << std::endl;
1157  return retval;
1158  }
1159  else
1160  std::cout << " --- PASSED ---" << std::endl;
1161 
1163 
1164  std::cout << "* m = slice, v1 = full" << std::endl;
1165  retval = test_solve<NumericT>(epsilon,
1166  std_m2, std_v1,
1167  vcl_m2_slice, vcl_v1_native);
1168  if (retval == EXIT_FAILURE)
1169  {
1170  std::cout << " --- FAILED! ---" << std::endl;
1171  return retval;
1172  }
1173  else
1174  std::cout << " --- PASSED ---" << std::endl;
1175 
1176  std::cout << "* m = slice, v1 = range" << std::endl;
1177  retval = test_solve<NumericT>(epsilon,
1178  std_m2, std_v1,
1179  vcl_m2_slice, vcl_v1_range);
1180  if (retval == EXIT_FAILURE)
1181  {
1182  std::cout << " --- FAILED! ---" << std::endl;
1183  return retval;
1184  }
1185  else
1186  std::cout << " --- PASSED ---" << std::endl;
1187 
1188  std::cout << "* m = slice, v1 = slice" << std::endl;
1189  retval = test_solve<NumericT>(epsilon,
1190  std_m2, std_v1,
1191  vcl_m2_slice, vcl_v1_slice);
1192  if (retval == EXIT_FAILURE)
1193  {
1194  std::cout << " --- FAILED! ---" << std::endl;
1195  return retval;
1196  }
1197  else
1198  std::cout << " --- PASSED ---" << std::endl;
1199 
1200 
1201 
1202 
1203 
1204 
1205 
1207 
1208  //full solver:
1209  std::cout << "Full solver" << std::endl;
1210  unsigned int lu_dim = 100;
1211  std::vector<std::vector<NumericT> > square_matrix(lu_dim, std::vector<NumericT>(lu_dim));
1212  std::vector<NumericT> lu_rhs(lu_dim);
1213  std::vector<NumericT> lu_result(lu_dim);
1214  viennacl::matrix<NumericT, F> vcl_square_matrix(lu_dim, lu_dim);
1215  viennacl::vector<NumericT> vcl_lu_rhs(lu_dim);
1216 
1217  for (std::size_t i=0; i<lu_dim; ++i)
1218  for (std::size_t j=0; j<lu_dim; ++j)
1219  square_matrix[i][j] = -static_cast<NumericT>(0.5) * randomNumber();
1220 
1221  //put some more weight on diagonal elements:
1222  for (std::size_t j=0; j<lu_dim; ++j)
1223  {
1224  square_matrix[j][j] = static_cast<NumericT>(20.0) + randomNumber();
1225  lu_result[j] = NumericT(0.1) + randomNumber();
1226  }
1227 
1228  for (std::size_t i=0; i<lu_dim; ++i)
1229  for (std::size_t j=0; j<lu_dim; ++j)
1230  lu_rhs[i] += square_matrix[i][j] * lu_result[j];
1231 
1232  viennacl::copy(square_matrix, vcl_square_matrix);
1233  viennacl::copy(lu_rhs, vcl_lu_rhs);
1234 
1235  // ViennaCL:
1236  viennacl::linalg::lu_factorize(vcl_square_matrix);
1237  viennacl::linalg::lu_substitute(vcl_square_matrix, vcl_lu_rhs);
1238 
1239  if ( std::fabs(diff(lu_result, vcl_lu_rhs)) > epsilon )
1240  {
1241  std::cout << "# Error at operation: dense solver" << std::endl;
1242  std::cout << " diff: " << std::fabs(diff(lu_rhs, vcl_lu_rhs)) << std::endl;
1243  retval = EXIT_FAILURE;
1244  }
1245 
1246 
1247 
1248  return retval;
1249 }
1250 //
1251 // -------------------------------------------------------------
1252 //
1253 int main()
1254 {
1255  std::cout << std::endl;
1256  std::cout << "----------------------------------------------" << std::endl;
1257  std::cout << "----------------------------------------------" << std::endl;
1258  std::cout << "## Test :: Matrix" << std::endl;
1259  std::cout << "----------------------------------------------" << std::endl;
1260  std::cout << "----------------------------------------------" << std::endl;
1261  std::cout << std::endl;
1262 
1263  int retval = EXIT_SUCCESS;
1264 
1265 // std::cout << std::endl;
1266 // std::cout << "----------------------------------------------" << std::endl;
1267 // std::cout << std::endl;
1268 // {
1269 // typedef float NumericT;
1270 // NumericT epsilon = NumericT(1.0E-3);
1271 // std::cout << "# Testing setup:" << std::endl;
1272 // std::cout << " eps: " << epsilon << std::endl;
1273 // std::cout << " numeric: float" << std::endl;
1274 // std::cout << " layout: row-major" << std::endl;
1275 // retval = test<NumericT, viennacl::row_major>(epsilon);
1276 // if ( retval == EXIT_SUCCESS )
1277 // std::cout << "# Test passed" << std::endl;
1278 // else
1279 // return retval;
1280 // }
1281  std::cout << std::endl;
1282  std::cout << "----------------------------------------------" << std::endl;
1283  std::cout << std::endl;
1284  {
1285  typedef float NumericT;
1286  NumericT epsilon = NumericT(1.0E-3);
1287  std::cout << "# Testing setup:" << std::endl;
1288  std::cout << " eps: " << epsilon << std::endl;
1289  std::cout << " numeric: float" << std::endl;
1290  std::cout << " layout: column-major" << std::endl;
1291  retval = test<NumericT, viennacl::column_major>(epsilon);
1292  if ( retval == EXIT_SUCCESS )
1293  std::cout << "# Test passed" << std::endl;
1294  else
1295  return retval;
1296  }
1297  std::cout << std::endl;
1298  std::cout << "----------------------------------------------" << std::endl;
1299  std::cout << std::endl;
1300 
1301 
1302 #ifdef VIENNACL_WITH_OPENCL
1303  if ( viennacl::ocl::current_device().double_support() )
1304 #endif
1305  {
1306  {
1307  typedef double NumericT;
1308  NumericT epsilon = 1.0E-11;
1309  std::cout << "# Testing setup:" << std::endl;
1310  std::cout << " eps: " << epsilon << std::endl;
1311  std::cout << " numeric: double" << std::endl;
1312  std::cout << " layout: row-major" << std::endl;
1313  retval = test<NumericT, viennacl::row_major>(epsilon);
1314  if ( retval == EXIT_SUCCESS )
1315  std::cout << "# Test passed" << std::endl;
1316  else
1317  return retval;
1318  }
1319  std::cout << std::endl;
1320  std::cout << "----------------------------------------------" << std::endl;
1321  std::cout << std::endl;
1322  {
1323  typedef double NumericT;
1324  NumericT epsilon = 1.0E-11;
1325  std::cout << "# Testing setup:" << std::endl;
1326  std::cout << " eps: " << epsilon << std::endl;
1327  std::cout << " numeric: double" << std::endl;
1328  std::cout << " layout: column-major" << std::endl;
1329  retval = test<NumericT, viennacl::column_major>(epsilon);
1330  if ( retval == EXIT_SUCCESS )
1331  std::cout << "# Test passed" << std::endl;
1332  else
1333  return retval;
1334  }
1335  std::cout << std::endl;
1336  std::cout << "----------------------------------------------" << std::endl;
1337  std::cout << std::endl;
1338  }
1339 
1340  std::cout << std::endl;
1341  std::cout << "------- Test completed --------" << std::endl;
1342  std::cout << std::endl;
1343 
1344 
1345  return retval;
1346 }
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::matrix_slice
Class for representing strided submatrices of a bigger matrix A.
Definition: forwards.h:443
trans
std::vector< std::vector< NumericT > > trans(std::vector< std::vector< NumericT > > const &A)
Definition: blas3_solve.cpp:195
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::linalg::row_sum
viennacl::vector_expression< const viennacl::matrix_base< NumericT >, const viennacl::matrix_base< NumericT >, viennacl::op_row_sum > row_sum(viennacl::matrix_base< NumericT > const &A)
User interface function for computing the sum of all elements of each row of a matrix.
Definition: sum.hpp:77
viennacl::linalg::lu_substitute
void lu_substitute(matrix< NumericT, F1, AlignmentV1 > const &A, matrix< NumericT, F2, AlignmentV2 > &B)
LU substitution for the system LU = rhs.
Definition: lu.hpp:201
viennacl::linalg::lower_tag
A tag class representing a lower triangular matrix.
Definition: forwards.h:849
inplace_solve_upper
void inplace_solve_upper(std::vector< std::vector< NumericT > > const &A, std::vector< NumericT > &b, bool unit_diagonal)
Definition: matrix_vector.cpp:422
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
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
viennacl::linalg::solve
VectorT solve(MatrixT const &matrix, VectorT const &rhs, bicgstab_tag const &tag, PreconditionerT const &precond)
Definition: bicgstab.hpp:496
main
int main()
Definition: matrix_vector.cpp:1253
NumericT
float NumericT
Definition: bisect.cpp:40
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
inplace_solve
void inplace_solve(std::vector< std::vector< NumericT > > const &A, std::vector< NumericT > &b, viennacl::linalg::upper_tag)
Definition: matrix_vector.cpp:434
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::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
test_solve
int test_solve(Epsilon const &epsilon, STLMatrixType &std_m1, STLVectorType &std_v1, VCLMatrixType &vcl_m1, VCLVectorType1 &vcl_v1)
Definition: matrix_vector.cpp:481
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
viennacl::diag
vector_expression< const matrix_base< NumericT >, const int, op_matrix_diag > diag(const matrix_base< NumericT > &A, int k=0)
Definition: matrix.hpp:895
lu.hpp
Implementations of LU factorization for row-major and column-major dense matrices.
direct_solve.hpp
Implementations of dense direct solvers are found here.
viennacl::row
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
viennacl::linalg::outer_prod
viennacl::matrix_expression< const vector_base< NumericT >, const vector_base< NumericT >, op_prod > outer_prod(const vector_base< NumericT > &vec1, const vector_base< NumericT > &vec2)
Returns a proxy class for the operation mat += vec1 * vec2^T, i.e. a rank 1 update.
Definition: matrix_operations.hpp:840
sum.hpp
Stub routines for the summation of elements in a vector, or all elements in either a row or column of...
v2
viennacl::vector< int > v2
Definition: global_variables.cpp:61
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
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::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::unit_lower_tag
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:859
test_prod_rank1
int test_prod_rank1(Epsilon const &epsilon, STLMatrixType &std_m1, STLVectorType &std_v1, STLVectorType &std_v2, STLMatrixType &std_m2, VCLMatrixType &vcl_m1, VCLVectorType1 &vcl_v1, VCLVectorType2 &vcl_v2, VCLMatrixType &vcl_m2)
Definition: matrix_vector.cpp:104
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: matrix_vector.cpp:47
inplace_solve_lower
void inplace_solve_lower(std::vector< std::vector< NumericT > > const &A, std::vector< NumericT > &b, bool unit_diagonal)
Definition: matrix_vector.cpp:447
viennacl::matrix_range
Class for representing non-strided submatrices of a bigger matrix A.
Definition: forwards.h:440
viennacl::column
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_column > column(const matrix_base< NumericT, F > &A, unsigned int j)
Definition: matrix.hpp:918
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::linalg::column_sum
viennacl::vector_expression< const viennacl::matrix_base< NumericT >, const viennacl::matrix_base< NumericT >, viennacl::op_col_sum > column_sum(viennacl::matrix_base< NumericT > const &A)
User interface function for computing the sum of all elements of each column of a matrix...
Definition: sum.hpp:109
scalar.hpp
Implementation of the ViennaCL scalar class.
test
int test(Epsilon const &epsilon)
Definition: matrix_vector.cpp:606
viennacl::linalg::lu_factorize
void lu_factorize(matrix< NumericT, viennacl::row_major > &A)
LU factorization of a row-major dense matrix.
Definition: lu.hpp:42
viennacl::linalg::unit_upper_tag
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:864
solve
std::vector< NumericT > solve(std::vector< std::vector< NumericT > > const &A, std::vector< NumericT > const &b, TagT)
Definition: matrix_vector.cpp:471