ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
matrix_vector.cpp
Go to the documentation of this file.
1 /* =========================================================================
2  Copyright (c) 2010-2015, Institute for Microelectronics,
3  Institute for Analysis and Scientific Computing,
4  TU Wien.
5  Portions of this software are copyright by UChicago Argonne, LLC.
6 
7  -----------------
8  ViennaCL - The Vienna Computing Library
9  -----------------
10 
11  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12 
13  (A list of authors and contributors can be found in the PDF manual)
14 
15  License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17 
18 
23 //
24 // *** System
25 //
26 #include <iostream>
27 
28 //
29 // *** Boost
30 //
31 #include <boost/numeric/ublas/io.hpp>
32 #include <boost/numeric/ublas/triangular.hpp>
33 #include <boost/numeric/ublas/matrix_sparse.hpp>
34 #include <boost/numeric/ublas/matrix.hpp>
35 #include <boost/numeric/ublas/matrix_proxy.hpp>
36 #include <boost/numeric/ublas/lu.hpp>
37 #include <boost/numeric/ublas/io.hpp>
38 #include <boost/numeric/ublas/vector_proxy.hpp>
39 
40 //
41 // *** ViennaCL
42 //
43 //#define VIENNACL_DEBUG_ALL
44 #define VIENNACL_WITH_UBLAS 1
45 #include "viennacl/scalar.hpp"
46 #include "viennacl/matrix.hpp"
47 #include "viennacl/vector.hpp"
48 #include "viennacl/linalg/prod.hpp"
49 #include "viennacl/linalg/norm_2.hpp"
50 #include "viennacl/linalg/direct_solve.hpp"
51 #include "viennacl/linalg/lu.hpp"
52 #include "viennacl/linalg/sum.hpp"
53 #include "viennacl/tools/random.hpp"
54 
55 //
56 // -------------------------------------------------------------
57 //
58 using namespace boost::numeric;
59 //
60 // -------------------------------------------------------------
61 //
62 template<typename ScalarType>
63 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
64 {
65  viennacl::backend::finish();
66  if (s1 != s2)
67  return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
68  return 0;
69 }
70 
71 template<typename ScalarType, typename VCLVectorType>
72 ScalarType diff(ublas::vector<ScalarType> const & v1, VCLVectorType const & v2)
73 {
74  ublas::vector<ScalarType> v2_cpu(v2.size());
75  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
76  viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
77 
78  for (unsigned int i=0;i<v1.size(); ++i)
79  {
80  if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
81  v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
82  else
83  v2_cpu[i] = 0.0;
84  }
85 
86  return norm_inf(v2_cpu);
87 }
88 
89 template<typename ScalarType, typename VCLMatrixType>
90 ScalarType diff(ublas::matrix<ScalarType> const & mat1, VCLMatrixType const & mat2)
91 {
92  ublas::matrix<ScalarType> mat2_cpu(mat2.size1(), mat2.size2());
93  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
94  viennacl::copy(mat2, mat2_cpu);
95  ScalarType ret = 0;
96  ScalarType act = 0;
97 
98  for (unsigned int i = 0; i < mat2_cpu.size1(); ++i)
99  {
100  for (unsigned int j = 0; j < mat2_cpu.size2(); ++j)
101  {
102  act = std::fabs(mat2_cpu(i,j) - mat1(i,j)) / std::max( std::fabs(mat2_cpu(i, j)), std::fabs(mat1(i,j)) );
103  if (act > ret)
104  ret = act;
105  }
106  }
107  //std::cout << ret << std::endl;
108  return ret;
109 }
110 //
111 // -------------------------------------------------------------
112 //
113 
114 template<typename NumericT, typename Epsilon,
115  typename UblasMatrixType, typename UblasVectorType,
116  typename VCLMatrixType, typename VCLVectorType1, typename VCLVectorType2>
117 int test_prod_rank1(Epsilon const & epsilon,
118  UblasMatrixType & ublas_m1, UblasVectorType & ublas_v1, UblasVectorType & ublas_v2, UblasMatrixType & ublas_m2,
119  VCLMatrixType & vcl_m1, VCLVectorType1 & vcl_v1, VCLVectorType2 & vcl_v2, VCLMatrixType & vcl_m2)
120 {
121  int retval = EXIT_SUCCESS;
122 
123  // sync data:
124  ublas_v1 = ublas::scalar_vector<NumericT>(ublas_v1.size(), NumericT(0.1234));
125  ublas_v2 = ublas::scalar_vector<NumericT>(ublas_v2.size(), NumericT(0.4321));
126  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
127  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
128  viennacl::copy(ublas_m1, vcl_m1);
129 
130  // --------------------------------------------------------------------------
131  std::cout << "Rank 1 update" << std::endl;
132 
133  ublas_m1 += ublas::outer_prod(ublas_v1, ublas_v2);
134  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
135  if ( std::fabs(diff(ublas_m1, vcl_m1)) > epsilon )
136  {
137  std::cout << "# Error at operation: rank 1 update" << std::endl;
138  std::cout << " diff: " << std::fabs(diff(ublas_m1, vcl_m1)) << std::endl;
139  return EXIT_FAILURE;
140  }
141 
142 
143 
144  // --------------------------------------------------------------------------
145  std::cout << "Scaled rank 1 update - CPU Scalar" << std::endl;
146  ublas_m1 += NumericT(4.2) * ublas::outer_prod(ublas_v1, ublas_v2);
147  vcl_m1 += NumericT(2.1) * viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
148  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2) * NumericT(2.1); //check proper compilation
149  if ( std::fabs(diff(ublas_m1, vcl_m1)) > epsilon )
150  {
151  std::cout << "# Error at operation: scaled rank 1 update - CPU Scalar" << std::endl;
152  std::cout << " diff: " << std::fabs(diff(ublas_m1, vcl_m1)) << std::endl;
153  return EXIT_FAILURE;
154  }
155 
156  // --------------------------------------------------------------------------
157  std::cout << "Scaled rank 1 update - GPU Scalar" << std::endl;
158  ublas_m1 += NumericT(4.2) * ublas::outer_prod(ublas_v1, ublas_v2);
159  vcl_m1 += viennacl::scalar<NumericT>(NumericT(2.1)) * viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
160  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2) * viennacl::scalar<NumericT>(NumericT(2.1)); //check proper compilation
161  if ( std::fabs(diff(ublas_m1, vcl_m1)) > epsilon )
162  {
163  std::cout << "# Error at operation: scaled rank 1 update - GPU Scalar" << std::endl;
164  std::cout << " diff: " << std::fabs(diff(ublas_m1, vcl_m1)) << std::endl;
165  return EXIT_FAILURE;
166  }
167 
168  //reset vcl_matrix:
169  viennacl::copy(ublas_m1, vcl_m1);
170 
171  // --------------------------------------------------------------------------
172  std::cout << "Matrix-Vector product" << std::endl;
173  ublas_v1 = viennacl::linalg::prod(ublas_m1, ublas_v2);
174  vcl_v1 = viennacl::linalg::prod(vcl_m1, vcl_v2);
175 
176  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
177  {
178  std::cout << "# Error at operation: matrix-vector product" << std::endl;
179  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
180  retval = EXIT_FAILURE;
181  }
182  // --------------------------------------------------------------------------
183  std::cout << "Matrix-Vector product with scaled add" << std::endl;
184  NumericT alpha = static_cast<NumericT>(2.786);
185  NumericT beta = static_cast<NumericT>(1.432);
186  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
187  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
188 
189  ublas_v1 = alpha * viennacl::linalg::prod(ublas_m1, ublas_v2) + beta * ublas_v1;
190  vcl_v1 = alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) + beta * vcl_v1;
191 
192  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
193  {
194  std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
195  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
196  retval = EXIT_FAILURE;
197  }
198  // --------------------------------------------------------------------------
199  std::cout << "Matrix-Vector product with matrix expression" << std::endl;
200  ublas_v1 = ublas::prod(ublas_m1 + ublas_m1, ublas_v2);
201  vcl_v1 = viennacl::linalg::prod(vcl_m1 + vcl_m1, vcl_v2);
202 
203  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
204  {
205  std::cout << "# Error at operation: matrix-vector product" << std::endl;
206  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
207  retval = EXIT_FAILURE;
208  }
209  // --------------------------------------------------------------------------
210  std::cout << "Matrix-Vector product with vector expression" << std::endl;
211  ublas_v1 = ublas::prod(ublas_m1, NumericT(3) * ublas_v2);
212  vcl_v1 = viennacl::linalg::prod(vcl_m1, NumericT(3) * vcl_v2);
213 
214  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
215  {
216  std::cout << "# Error at operation: matrix-vector product" << std::endl;
217  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
218  retval = EXIT_FAILURE;
219  }
220  // --------------------------------------------------------------------------
221  std::cout << "Matrix-Vector product with matrix and vector expression" << std::endl;
222  ublas_v1 = ublas::prod(ublas_m1 + ublas_m1, ublas_v2 + ublas_v2);
223  vcl_v1 = viennacl::linalg::prod(vcl_m1 + vcl_m1, vcl_v2 + vcl_v2);
224 
225  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
226  {
227  std::cout << "# Error at operation: matrix-vector product" << std::endl;
228  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
229  retval = EXIT_FAILURE;
230  }
231  // --------------------------------------------------------------------------
232 
233  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
234  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
235 
236  std::cout << "Transposed Matrix-Vector product" << std::endl;
237  ublas_v2 = alpha * viennacl::linalg::prod(trans(ublas_m1), ublas_v1);
238  vcl_v2 = alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1);
239 
240  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
241  {
242  std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
243  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
244  retval = EXIT_FAILURE;
245  }
246 
247  std::cout << "Transposed Matrix-Vector product with scaled add" << std::endl;
248  ublas_v2 = alpha * viennacl::linalg::prod(trans(ublas_m1), ublas_v1) + beta * ublas_v2;
249  vcl_v2 = alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2;
250 
251  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
252  {
253  std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
254  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
255  retval = EXIT_FAILURE;
256  }
257  // --------------------------------------------------------------------------
258 
259  std::cout << "Row sum with matrix" << std::endl;
260  ublas_v1 = ublas::prod(ublas_m1, ublas::scalar_vector<NumericT>(ublas_m1.size2(), NumericT(1)));
261  vcl_v1 = viennacl::linalg::row_sum(vcl_m1);
262 
263  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
264  {
265  std::cout << "# Error at operation: row sum" << std::endl;
266  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
267  retval = EXIT_FAILURE;
268  }
269  // --------------------------------------------------------------------------
270 
271  std::cout << "Row sum with matrix expression" << std::endl;
272  ublas_v1 = ublas::prod(ublas_m1 + ublas_m1, ublas::scalar_vector<NumericT>(ublas_m1.size2(), NumericT(1)));
273  vcl_v1 = viennacl::linalg::row_sum(vcl_m1 + vcl_m1);
274 
275  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
276  {
277  std::cout << "# Error at operation: row sum (with expression)" << std::endl;
278  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
279  retval = EXIT_FAILURE;
280  }
281  // --------------------------------------------------------------------------
282 
283  std::cout << "Column sum with matrix" << std::endl;
284  ublas_v2 = ublas::prod(trans(ublas_m1), ublas::scalar_vector<NumericT>(ublas_m1.size1(), NumericT(1)));
285  vcl_v2 = viennacl::linalg::column_sum(vcl_m1);
286 
287  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
288  {
289  std::cout << "# Error at operation: column sum" << std::endl;
290  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
291  retval = EXIT_FAILURE;
292  }
293  // --------------------------------------------------------------------------
294 
295  std::cout << "Column sum with matrix expression" << std::endl;
296  ublas_v2 = ublas::prod(trans(ublas_m1 + ublas_m1), ublas::scalar_vector<NumericT>(ublas_m1.size1(), NumericT(1)));
297  vcl_v2 = viennacl::linalg::column_sum(vcl_m1 + vcl_m1);
298 
299  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
300  {
301  std::cout << "# Error at operation: column sum (with expression)" << std::endl;
302  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
303  retval = EXIT_FAILURE;
304  }
305  // --------------------------------------------------------------------------
306 
307 
308  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
309  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
310 
311  std::cout << "Row extraction from matrix" << std::endl;
312  ublas_v2 = row(ublas_m1, std::size_t(7));
313  vcl_v2 = row(vcl_m1, std::size_t(7));
314 
315  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
316  {
317  std::cout << "# Error at operation: diagonal extraction from matrix" << std::endl;
318  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
319  retval = EXIT_FAILURE;
320  }
321 
322  std::cout << "Column extraction from matrix" << std::endl;
323  ublas_v1 = column(ublas_m1, std::size_t(7));
324  vcl_v1 = column(vcl_m1, std::size_t(7));
325 
326  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
327  {
328  std::cout << "# Error at operation: diagonal extraction from matrix" << std::endl;
329  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
330  retval = EXIT_FAILURE;
331  }
332 
333  // --------------------------------------------------------------------------
334 
335  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
336  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
337  viennacl::copy(ublas_m2, vcl_m2);
338  UblasMatrixType A = ublas_m2;
339 
340  std::cout << "Diagonal extraction from matrix" << std::endl;
341  for (std::size_t i=0; i<ublas_m1.size2(); ++i)
342  ublas_v2[i] = ublas_m1(i + 3, i);
343  vcl_v2 = diag(vcl_m1, static_cast<int>(-3));
344 
345  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
346  {
347  std::cout << "# Error at operation: diagonal extraction from matrix" << std::endl;
348  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
349  retval = EXIT_FAILURE;
350  }
351 
352  std::cout << "Matrix diagonal assignment from vector" << std::endl;
353  A = ublas::scalar_matrix<NumericT>(A.size1(), A.size2(), NumericT(0));
354  for (std::size_t i=0; i<ublas_m1.size2(); ++i)
355  A(i + (A.size1() - ublas_m1.size2()), i) = ublas_v2[i];
356  vcl_m2 = diag(vcl_v2, static_cast<int>(ublas_m1.size2()) - static_cast<int>(A.size1()));
357 
358  if ( std::fabs(diff(A, vcl_m2)) > epsilon )
359  {
360  std::cout << "# Error at operation: Matrix assignment from diagonal" << std::endl;
361  std::cout << " diff: " << std::fabs(diff(A, vcl_m2)) << std::endl;
362  retval = EXIT_FAILURE;
363  }
364  // --------------------------------------------------------------------------
365 
366  return retval;
367 }
368 
369 
370 
371 template<typename NumericT, typename Epsilon,
372  typename UblasMatrixType, typename UblasVectorType,
373  typename VCLMatrixType, typename VCLVectorType1>
374 int test_solve(Epsilon const & epsilon,
375  UblasMatrixType & ublas_m1, UblasVectorType & ublas_v1,
376  VCLMatrixType & vcl_m1, VCLVectorType1 & vcl_v1)
377 {
378  int retval = EXIT_SUCCESS;
379 
380  // sync data:
381  //viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
382  viennacl::copy(ublas_v1, vcl_v1);
383  viennacl::copy(ublas_m1, vcl_m1);
384 
386 
387  //upper triangular:
388  std::cout << "Upper triangular solver" << std::endl;
389  ublas_v1 = ublas::solve(ublas_m1, ublas_v1, ublas::upper_tag());
390  vcl_v1 = viennacl::linalg::solve(vcl_m1, vcl_v1, viennacl::linalg::upper_tag());
391  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
392  {
393  std::cout << "# Error at operation: upper triangular solver" << std::endl;
394  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
395  retval = EXIT_FAILURE;
396  }
397 
398  //upper unit triangular:
399  std::cout << "Upper unit triangular solver" << std::endl;
400  viennacl::copy(ublas_v1, vcl_v1);
401  ublas_v1 = ublas::solve(ublas_m1, ublas_v1, ublas::unit_upper_tag());
402  vcl_v1 = viennacl::linalg::solve(vcl_m1, vcl_v1, viennacl::linalg::unit_upper_tag());
403  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
404  {
405  std::cout << "# Error at operation: unit upper triangular solver" << std::endl;
406  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
407  retval = EXIT_FAILURE;
408  }
409 
410  //lower triangular:
411  std::cout << "Lower triangular solver" << std::endl;
412  viennacl::copy(ublas_v1, vcl_v1);
413  ublas_v1 = ublas::solve(ublas_m1, ublas_v1, ublas::lower_tag());
414  vcl_v1 = viennacl::linalg::solve(vcl_m1, vcl_v1, viennacl::linalg::lower_tag());
415  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
416  {
417  std::cout << "# Error at operation: lower triangular solver" << std::endl;
418  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
419  retval = EXIT_FAILURE;
420  }
421 
422  //lower unit triangular:
423  std::cout << "Lower unit triangular solver" << std::endl;
424  viennacl::copy(ublas_v1, vcl_v1);
425  ublas_v1 = ublas::solve(ublas_m1, ublas_v1, ublas::unit_lower_tag());
426  vcl_v1 = viennacl::linalg::solve(vcl_m1, vcl_v1, viennacl::linalg::unit_lower_tag());
427  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
428  {
429  std::cout << "# Error at operation: unit lower triangular solver" << std::endl;
430  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
431  retval = EXIT_FAILURE;
432  }
433 
434 
435 
436 
437 
438  //transposed upper triangular:
439  std::cout << "Transposed upper triangular solver" << std::endl;
440  viennacl::copy(ublas_v1, vcl_v1);
441  ublas_v1 = ublas::solve(trans(ublas_m1), ublas_v1, ublas::upper_tag());
442  vcl_v1 = viennacl::linalg::solve(trans(vcl_m1), vcl_v1, viennacl::linalg::upper_tag());
443  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
444  {
445  std::cout << "# Error at operation: upper triangular solver" << std::endl;
446  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
447  retval = EXIT_FAILURE;
448  }
449 
450  //transposed upper unit triangular:
451  std::cout << "Transposed unit upper triangular solver" << std::endl;
452  viennacl::copy(ublas_v1, vcl_v1);
453  ublas_v1 = ublas::solve(trans(ublas_m1), ublas_v1, ublas::unit_upper_tag());
454  vcl_v1 = viennacl::linalg::solve(trans(vcl_m1), vcl_v1, viennacl::linalg::unit_upper_tag());
455  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
456  {
457  std::cout << "# Error at operation: unit upper triangular solver" << std::endl;
458  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
459  retval = EXIT_FAILURE;
460  }
461 
462  //transposed lower triangular:
463  std::cout << "Transposed lower triangular solver" << std::endl;
464  viennacl::copy(ublas_v1, vcl_v1);
465  ublas_v1 = ublas::solve(trans(ublas_m1), ublas_v1, ublas::lower_tag());
466  vcl_v1 = viennacl::linalg::solve(trans(vcl_m1), vcl_v1, viennacl::linalg::lower_tag());
467  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
468  {
469  std::cout << "# Error at operation: lower triangular solver" << std::endl;
470  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
471  retval = EXIT_FAILURE;
472  }
473 
474  //transposed lower unit triangular:
475  std::cout << "Transposed unit lower triangular solver" << std::endl;
476  viennacl::copy(ublas_v1, vcl_v1);
477  ublas_v1 = ublas::solve(trans(ublas_m1), ublas_v1, ublas::unit_lower_tag());
478  vcl_v1 = viennacl::linalg::solve(trans(vcl_m1), vcl_v1, viennacl::linalg::unit_lower_tag());
479  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
480  {
481  std::cout << "# Error at operation: unit lower triangular solver" << std::endl;
482  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
483  retval = EXIT_FAILURE;
484  }
485 
486  return retval;
487 }
488 
489 
490 //
491 // -------------------------------------------------------------
492 //
493 template< typename NumericT, typename F, typename Epsilon >
494 int test(Epsilon const& epsilon)
495 {
496  int retval = EXIT_SUCCESS;
497 
498  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
499 
500  std::size_t num_rows = 141; //note: use num_rows > num_cols + 3 for diag() tests to work
501  std::size_t num_cols = 103;
502 
503  // --------------------------------------------------------------------------
504  ublas::vector<NumericT> ublas_v1(num_rows);
505  for (std::size_t i = 0; i < ublas_v1.size(); ++i)
506  ublas_v1(i) = randomNumber();
507  ublas::vector<NumericT> ublas_v2 = ublas::scalar_vector<NumericT>(num_cols, NumericT(3.1415));
508 
509 
510  ublas::matrix<NumericT> ublas_m1(ublas_v1.size(), ublas_v2.size());
511 
512  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
513  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
514  ublas_m1(i,j) = static_cast<NumericT>(0.1) * randomNumber();
515 
516 
517  ublas::matrix<NumericT> ublas_m2(ublas_v1.size(), ublas_v1.size());
518 
519  for (std::size_t i = 0; i < ublas_m2.size1(); ++i)
520  {
521  for (std::size_t j = 0; j < ublas_m2.size2(); ++j)
522  ublas_m2(i,j) = static_cast<NumericT>(-0.1) * randomNumber();
523  ublas_m2(i, i) = static_cast<NumericT>(2) + randomNumber();
524  }
525 
526 
527  viennacl::vector<NumericT> vcl_v1_native(ublas_v1.size());
528  viennacl::vector<NumericT> vcl_v1_large(4 * ublas_v1.size());
529  viennacl::vector_range< viennacl::vector<NumericT> > vcl_v1_range(vcl_v1_large, viennacl::range(3, ublas_v1.size() + 3));
530  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v1_slice(vcl_v1_large, viennacl::slice(2, 3, ublas_v1.size()));
531 
532  viennacl::vector<NumericT> vcl_v2_native(ublas_v2.size());
533  viennacl::vector<NumericT> vcl_v2_large(4 * ublas_v2.size());
534  viennacl::vector_range< viennacl::vector<NumericT> > vcl_v2_range(vcl_v2_large, viennacl::range(8, ublas_v2.size() + 8));
535  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v2_slice(vcl_v2_large, viennacl::slice(6, 2, ublas_v2.size()));
536 
537  viennacl::matrix<NumericT, F> vcl_m1_native(ublas_m1.size1(), ublas_m1.size2());
538  viennacl::matrix<NumericT, F> vcl_m1_large(4 * ublas_m1.size1(), 4 * ublas_m1.size2());
539  viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m1_range(vcl_m1_large,
540  viennacl::range(8, ublas_m1.size1() + 8),
541  viennacl::range(ublas_m1.size2(), 2 * ublas_m1.size2()) );
542  viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m1_slice(vcl_m1_large,
543  viennacl::slice(6, 2, ublas_m1.size1()),
544  viennacl::slice(ublas_m1.size2(), 2, ublas_m1.size2()) );
545 
546  viennacl::matrix<NumericT, F> vcl_m2_native(ublas_m2.size1(), ublas_m2.size2());
547  viennacl::matrix<NumericT, F> vcl_m2_large(4 * ublas_m2.size1(), 4 * ublas_m2.size2());
548  viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m2_range(vcl_m2_large,
549  viennacl::range(8, ublas_m2.size1() + 8),
550  viennacl::range(ublas_m2.size2(), 2 * ublas_m2.size2()) );
551  viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m2_slice(vcl_m2_large,
552  viennacl::slice(6, 2, ublas_m2.size1()),
553  viennacl::slice(ublas_m2.size2(), 2, ublas_m2.size2()) );
554 
555 
556 /* std::cout << "Matrix resizing (to larger)" << std::endl;
557  matrix.resize(2*num_rows, 2*num_cols, true);
558  for (unsigned int i = 0; i < matrix.size1(); ++i)
559  {
560  for (unsigned int j = (i<result.size() ? rhs.size() : 0); j < matrix.size2(); ++j)
561  matrix(i,j) = 0;
562  }
563  vcl_matrix.resize(2*num_rows, 2*num_cols, true);
564  viennacl::copy(vcl_matrix, matrix);
565  if ( std::fabs(diff(matrix, vcl_matrix)) > epsilon )
566  {
567  std::cout << "# Error at operation: matrix resize (to larger)" << std::endl;
568  std::cout << " diff: " << std::fabs(diff(matrix, vcl_matrix)) << std::endl;
569  return EXIT_FAILURE;
570  }
571 
572  matrix(12, 14) = NumericT(1.9);
573  matrix(19, 16) = NumericT(1.0);
574  matrix (13, 15) = NumericT(-9);
575  vcl_matrix(12, 14) = NumericT(1.9);
576  vcl_matrix(19, 16) = NumericT(1.0);
577  vcl_matrix (13, 15) = NumericT(-9);
578 
579  std::cout << "Matrix resizing (to smaller)" << std::endl;
580  matrix.resize(result.size(), rhs.size(), true);
581  vcl_matrix.resize(result.size(), rhs.size(), true);
582  if ( std::fabs(diff(matrix, vcl_matrix)) > epsilon )
583  {
584  std::cout << "# Error at operation: matrix resize (to smaller)" << std::endl;
585  std::cout << " diff: " << std::fabs(diff(matrix, vcl_matrix)) << std::endl;
586  return EXIT_FAILURE;
587  }
588  */
589 
590  //
591  // Run a bunch of tests for rank-1-updates, matrix-vector products
592  //
593  std::cout << "------------ Testing rank-1-updates and matrix-vector products ------------------" << std::endl;
594 
595  std::cout << "* m = full, v1 = full, v2 = full" << std::endl;
596  retval = test_prod_rank1<NumericT>(epsilon,
597  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
598  vcl_m1_native, vcl_v1_native, vcl_v2_native, vcl_m2_native);
599  if (retval == EXIT_FAILURE)
600  {
601  std::cout << " --- FAILED! ---" << std::endl;
602  return retval;
603  }
604  else
605  std::cout << " --- PASSED ---" << std::endl;
606 
607 
608  std::cout << "* m = full, v1 = full, v2 = range" << std::endl;
609  retval = test_prod_rank1<NumericT>(epsilon,
610  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
611  vcl_m1_native, vcl_v1_native, vcl_v2_range, vcl_m2_native);
612  if (retval == EXIT_FAILURE)
613  {
614  std::cout << " --- FAILED! ---" << std::endl;
615  return retval;
616  }
617  else
618  std::cout << " --- PASSED ---" << std::endl;
619 
620 
621  std::cout << "* m = full, v1 = full, v2 = slice" << std::endl;
622  retval = test_prod_rank1<NumericT>(epsilon,
623  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
624  vcl_m1_native, vcl_v1_native, vcl_v2_slice, vcl_m2_native);
625  if (retval == EXIT_FAILURE)
626  {
627  std::cout << " --- FAILED! ---" << std::endl;
628  return retval;
629  }
630  else
631  std::cout << " --- PASSED ---" << std::endl;
632 
633 
634  // v1 = range
635 
636 
637  std::cout << "* m = full, v1 = range, v2 = full" << std::endl;
638  retval = test_prod_rank1<NumericT>(epsilon,
639  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
640  vcl_m1_native, vcl_v1_range, vcl_v2_native, vcl_m2_native);
641  if (retval == EXIT_FAILURE)
642  {
643  std::cout << " --- FAILED! ---" << std::endl;
644  return retval;
645  }
646  else
647  std::cout << " --- PASSED ---" << std::endl;
648 
649 
650  std::cout << "* m = full, v1 = range, v2 = range" << std::endl;
651  retval = test_prod_rank1<NumericT>(epsilon,
652  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
653  vcl_m1_native, vcl_v1_range, vcl_v2_range, vcl_m2_native);
654  if (retval == EXIT_FAILURE)
655  {
656  std::cout << " --- FAILED! ---" << std::endl;
657  return retval;
658  }
659  else
660  std::cout << " --- PASSED ---" << std::endl;
661 
662 
663  std::cout << "* m = full, v1 = range, v2 = slice" << std::endl;
664  retval = test_prod_rank1<NumericT>(epsilon,
665  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
666  vcl_m1_native, vcl_v1_range, vcl_v2_slice, vcl_m2_native);
667  if (retval == EXIT_FAILURE)
668  {
669  std::cout << " --- FAILED! ---" << std::endl;
670  return retval;
671  }
672  else
673  std::cout << " --- PASSED ---" << std::endl;
674 
675 
676 
677  // v1 = slice
678 
679  std::cout << "* m = full, v1 = slice, v2 = full" << std::endl;
680  retval = test_prod_rank1<NumericT>(epsilon,
681  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
682  vcl_m1_native, vcl_v1_slice, vcl_v2_native, vcl_m2_native);
683  if (retval == EXIT_FAILURE)
684  {
685  std::cout << " --- FAILED! ---" << std::endl;
686  return retval;
687  }
688  else
689  std::cout << " --- PASSED ---" << std::endl;
690 
691 
692  std::cout << "* m = full, v1 = slice, v2 = range" << std::endl;
693  retval = test_prod_rank1<NumericT>(epsilon,
694  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
695  vcl_m1_native, vcl_v1_slice, vcl_v2_range, vcl_m2_native);
696  if (retval == EXIT_FAILURE)
697  {
698  std::cout << " --- FAILED! ---" << std::endl;
699  return retval;
700  }
701  else
702  std::cout << " --- PASSED ---" << std::endl;
703 
704 
705  std::cout << "* m = full, v1 = slice, v2 = slice" << std::endl;
706  retval = test_prod_rank1<NumericT>(epsilon,
707  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
708  vcl_m1_native, vcl_v1_slice, vcl_v2_slice, vcl_m2_native);
709  if (retval == EXIT_FAILURE)
710  {
711  std::cout << " --- FAILED! ---" << std::endl;
712  return retval;
713  }
714  else
715  std::cout << " --- PASSED ---" << std::endl;
716 
717 
719 
720  std::cout << "* m = range, v1 = full, v2 = full" << std::endl;
721  retval = test_prod_rank1<NumericT>(epsilon,
722  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
723  vcl_m1_range, vcl_v1_native, vcl_v2_native, vcl_m2_range);
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 = range, v1 = full, v2 = range" << std::endl;
734  retval = test_prod_rank1<NumericT>(epsilon,
735  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
736  vcl_m1_range, vcl_v1_native, vcl_v2_range, vcl_m2_range);
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  std::cout << "* m = range, v1 = full, v2 = slice" << std::endl;
747  retval = test_prod_rank1<NumericT>(epsilon,
748  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
749  vcl_m1_range, vcl_v1_native, vcl_v2_slice, vcl_m2_range);
750  if (retval == EXIT_FAILURE)
751  {
752  std::cout << " --- FAILED! ---" << std::endl;
753  return retval;
754  }
755  else
756  std::cout << " --- PASSED ---" << std::endl;
757 
758 
759  // v1 = range
760 
761 
762  std::cout << "* m = range, v1 = range, v2 = full" << std::endl;
763  retval = test_prod_rank1<NumericT>(epsilon,
764  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
765  vcl_m1_range, vcl_v1_range, vcl_v2_native, vcl_m2_range);
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 = range, v1 = range, v2 = range" << std::endl;
776  retval = test_prod_rank1<NumericT>(epsilon,
777  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
778  vcl_m1_range, vcl_v1_range, vcl_v2_range, vcl_m2_range);
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  std::cout << "* m = range, v1 = range, v2 = slice" << std::endl;
789  retval = test_prod_rank1<NumericT>(epsilon,
790  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
791  vcl_m1_range, vcl_v1_range, vcl_v2_slice, vcl_m2_range);
792  if (retval == EXIT_FAILURE)
793  {
794  std::cout << " --- FAILED! ---" << std::endl;
795  return retval;
796  }
797  else
798  std::cout << " --- PASSED ---" << std::endl;
799 
800 
801 
802  // v1 = slice
803 
804  std::cout << "* m = range, v1 = slice, v2 = full" << std::endl;
805  retval = test_prod_rank1<NumericT>(epsilon,
806  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
807  vcl_m1_range, vcl_v1_slice, vcl_v2_native, vcl_m2_range);
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 = range, v1 = slice, v2 = range" << std::endl;
818  retval = test_prod_rank1<NumericT>(epsilon,
819  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
820  vcl_m1_range, vcl_v1_slice, vcl_v2_range, vcl_m2_range);
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 
830  std::cout << "* m = range, v1 = slice, v2 = slice" << std::endl;
831  retval = test_prod_rank1<NumericT>(epsilon,
832  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
833  vcl_m1_range, vcl_v1_slice, vcl_v2_slice, vcl_m2_range);
834  if (retval == EXIT_FAILURE)
835  {
836  std::cout << " --- FAILED! ---" << std::endl;
837  return retval;
838  }
839  else
840  std::cout << " --- PASSED ---" << std::endl;
841 
842 
844 
845  std::cout << "* m = slice, v1 = full, v2 = full" << std::endl;
846  retval = test_prod_rank1<NumericT>(epsilon,
847  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
848  vcl_m1_slice, vcl_v1_native, vcl_v2_native, vcl_m2_slice);
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 = slice, v1 = full, v2 = range" << std::endl;
859  retval = test_prod_rank1<NumericT>(epsilon,
860  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
861  vcl_m1_slice, vcl_v1_native, vcl_v2_range, vcl_m2_slice);
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  std::cout << "* m = slice, v1 = full, v2 = slice" << std::endl;
872  retval = test_prod_rank1<NumericT>(epsilon,
873  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
874  vcl_m1_slice, vcl_v1_native, vcl_v2_slice, vcl_m2_slice);
875  if (retval == EXIT_FAILURE)
876  {
877  std::cout << " --- FAILED! ---" << std::endl;
878  return retval;
879  }
880  else
881  std::cout << " --- PASSED ---" << std::endl;
882 
883 
884  // v1 = range
885 
886 
887  std::cout << "* m = slice, v1 = range, v2 = full" << std::endl;
888  retval = test_prod_rank1<NumericT>(epsilon,
889  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
890  vcl_m1_slice, vcl_v1_range, vcl_v2_native, vcl_m2_slice);
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 = slice, v1 = range, v2 = range" << std::endl;
901  retval = test_prod_rank1<NumericT>(epsilon,
902  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
903  vcl_m1_slice, vcl_v1_range, vcl_v2_range, vcl_m2_slice);
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  std::cout << "* m = slice, v1 = range, v2 = slice" << std::endl;
914  retval = test_prod_rank1<NumericT>(epsilon,
915  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
916  vcl_m1_slice, vcl_v1_range, vcl_v2_slice, vcl_m2_slice);
917  if (retval == EXIT_FAILURE)
918  {
919  std::cout << " --- FAILED! ---" << std::endl;
920  return retval;
921  }
922  else
923  std::cout << " --- PASSED ---" << std::endl;
924 
925 
926 
927  // v1 = slice
928 
929  std::cout << "* m = slice, v1 = slice, v2 = full" << std::endl;
930  retval = test_prod_rank1<NumericT>(epsilon,
931  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
932  vcl_m1_slice, vcl_v1_slice, vcl_v2_native, vcl_m2_slice);
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 = slice, v1 = slice, v2 = range" << std::endl;
943  retval = test_prod_rank1<NumericT>(epsilon,
944  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
945  vcl_m1_slice, vcl_v1_slice, vcl_v2_range, vcl_m2_slice);
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 
955  std::cout << "* m = slice, v1 = slice, v2 = slice" << std::endl;
956  retval = test_prod_rank1<NumericT>(epsilon,
957  ublas_m1, ublas_v1, ublas_v2, ublas_m2,
958  vcl_m1_slice, vcl_v1_slice, vcl_v2_slice, vcl_m2_slice);
959  if (retval == EXIT_FAILURE)
960  {
961  std::cout << " --- FAILED! ---" << std::endl;
962  return retval;
963  }
964  else
965  std::cout << " --- PASSED ---" << std::endl;
966 
967 
968 
969  //
970  // Testing triangular solve() routines
971  //
972 
973  std::cout << "------------ Testing triangular solves ------------------" << std::endl;
974 
975  std::cout << "* m = full, v1 = full" << std::endl;
976  retval = test_solve<NumericT>(epsilon,
977  ublas_m2, ublas_v1,
978  vcl_m2_native, vcl_v1_native);
979  if (retval == EXIT_FAILURE)
980  {
981  std::cout << " --- FAILED! ---" << std::endl;
982  return retval;
983  }
984  else
985  std::cout << " --- PASSED ---" << std::endl;
986 
987  std::cout << "* m = full, v1 = range" << std::endl;
988  retval = test_solve<NumericT>(epsilon,
989  ublas_m2, ublas_v1,
990  vcl_m2_native, vcl_v1_range);
991  if (retval == EXIT_FAILURE)
992  {
993  std::cout << " --- FAILED! ---" << std::endl;
994  return retval;
995  }
996  else
997  std::cout << " --- PASSED ---" << std::endl;
998 
999  std::cout << "* m = full, v1 = slice" << std::endl;
1000  retval = test_solve<NumericT>(epsilon,
1001  ublas_m2, ublas_v1,
1002  vcl_m2_native, vcl_v1_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 
1012 
1013 
1014  std::cout << "* m = range, v1 = full" << std::endl;
1015  retval = test_solve<NumericT>(epsilon,
1016  ublas_m2, ublas_v1,
1017  vcl_m2_range, vcl_v1_native);
1018  if (retval == EXIT_FAILURE)
1019  {
1020  std::cout << " --- FAILED! ---" << std::endl;
1021  return retval;
1022  }
1023  else
1024  std::cout << " --- PASSED ---" << std::endl;
1025 
1026  std::cout << "* m = range, v1 = range" << std::endl;
1027  retval = test_solve<NumericT>(epsilon,
1028  ublas_m2, ublas_v1,
1029  vcl_m2_range, vcl_v1_range);
1030  if (retval == EXIT_FAILURE)
1031  {
1032  std::cout << " --- FAILED! ---" << std::endl;
1033  return retval;
1034  }
1035  else
1036  std::cout << " --- PASSED ---" << std::endl;
1037 
1038  std::cout << "* m = range, v1 = slice" << std::endl;
1039  retval = test_solve<NumericT>(epsilon,
1040  ublas_m2, ublas_v1,
1041  vcl_m2_range, vcl_v1_slice);
1042  if (retval == EXIT_FAILURE)
1043  {
1044  std::cout << " --- FAILED! ---" << std::endl;
1045  return retval;
1046  }
1047  else
1048  std::cout << " --- PASSED ---" << std::endl;
1049 
1051 
1052  std::cout << "* m = slice, v1 = full" << std::endl;
1053  retval = test_solve<NumericT>(epsilon,
1054  ublas_m2, ublas_v1,
1055  vcl_m2_slice, vcl_v1_native);
1056  if (retval == EXIT_FAILURE)
1057  {
1058  std::cout << " --- FAILED! ---" << std::endl;
1059  return retval;
1060  }
1061  else
1062  std::cout << " --- PASSED ---" << std::endl;
1063 
1064  std::cout << "* m = slice, v1 = range" << std::endl;
1065  retval = test_solve<NumericT>(epsilon,
1066  ublas_m2, ublas_v1,
1067  vcl_m2_slice, vcl_v1_range);
1068  if (retval == EXIT_FAILURE)
1069  {
1070  std::cout << " --- FAILED! ---" << std::endl;
1071  return retval;
1072  }
1073  else
1074  std::cout << " --- PASSED ---" << std::endl;
1075 
1076  std::cout << "* m = slice, v1 = slice" << std::endl;
1077  retval = test_solve<NumericT>(epsilon,
1078  ublas_m2, ublas_v1,
1079  vcl_m2_slice, vcl_v1_slice);
1080  if (retval == EXIT_FAILURE)
1081  {
1082  std::cout << " --- FAILED! ---" << std::endl;
1083  return retval;
1084  }
1085  else
1086  std::cout << " --- PASSED ---" << std::endl;
1087 
1088 
1089 
1090 
1091 
1092 
1093 
1095 
1096  //full solver:
1097  std::cout << "Full solver" << std::endl;
1098  unsigned int lu_dim = 100;
1099  ublas::matrix<NumericT> square_matrix(lu_dim, lu_dim);
1100  ublas::vector<NumericT> lu_rhs(lu_dim);
1101  viennacl::matrix<NumericT, F> vcl_square_matrix(lu_dim, lu_dim);
1102  viennacl::vector<NumericT> vcl_lu_rhs(lu_dim);
1103 
1104  for (std::size_t i=0; i<lu_dim; ++i)
1105  for (std::size_t j=0; j<lu_dim; ++j)
1106  square_matrix(i,j) = -static_cast<NumericT>(0.5) * randomNumber();
1107 
1108  //put some more weight on diagonal elements:
1109  for (std::size_t j=0; j<lu_dim; ++j)
1110  {
1111  square_matrix(j,j) = static_cast<NumericT>(20.0) + randomNumber();
1112  lu_rhs(j) = randomNumber();
1113  }
1114 
1115  viennacl::copy(square_matrix, vcl_square_matrix);
1116  viennacl::copy(lu_rhs, vcl_lu_rhs);
1117 
1118  //ublas::
1119  ublas::lu_factorize(square_matrix);
1120  ublas::inplace_solve (square_matrix, lu_rhs, ublas::unit_lower_tag ());
1121  ublas::inplace_solve (square_matrix, lu_rhs, ublas::upper_tag ());
1122 
1123  // ViennaCL:
1124  viennacl::linalg::lu_factorize(vcl_square_matrix);
1125  //viennacl::copy(square_matrix, vcl_square_matrix);
1126  viennacl::linalg::lu_substitute(vcl_square_matrix, vcl_lu_rhs);
1127 
1128  if ( std::fabs(diff(lu_rhs, vcl_lu_rhs)) > epsilon )
1129  {
1130  std::cout << "# Error at operation: dense solver" << std::endl;
1131  std::cout << " diff: " << std::fabs(diff(lu_rhs, vcl_lu_rhs)) << std::endl;
1132  retval = EXIT_FAILURE;
1133  }
1134 
1135 
1136 
1137  return retval;
1138 }
1139 //
1140 // -------------------------------------------------------------
1141 //
1142 int main()
1143 {
1144  std::cout << std::endl;
1145  std::cout << "----------------------------------------------" << std::endl;
1146  std::cout << "----------------------------------------------" << std::endl;
1147  std::cout << "## Test :: Matrix" << std::endl;
1148  std::cout << "----------------------------------------------" << std::endl;
1149  std::cout << "----------------------------------------------" << std::endl;
1150  std::cout << std::endl;
1151 
1152  int retval = EXIT_SUCCESS;
1153 
1154 // std::cout << std::endl;
1155 // std::cout << "----------------------------------------------" << std::endl;
1156 // std::cout << std::endl;
1157 // {
1158 // typedef float NumericT;
1159 // NumericT epsilon = NumericT(1.0E-3);
1160 // std::cout << "# Testing setup:" << std::endl;
1161 // std::cout << " eps: " << epsilon << std::endl;
1162 // std::cout << " numeric: float" << std::endl;
1163 // std::cout << " layout: row-major" << std::endl;
1164 // retval = test<NumericT, viennacl::row_major>(epsilon);
1165 // if ( retval == EXIT_SUCCESS )
1166 // std::cout << "# Test passed" << std::endl;
1167 // else
1168 // return retval;
1169 // }
1170  std::cout << std::endl;
1171  std::cout << "----------------------------------------------" << std::endl;
1172  std::cout << std::endl;
1173  {
1174  typedef float NumericT;
1175  NumericT epsilon = NumericT(1.0E-3);
1176  std::cout << "# Testing setup:" << std::endl;
1177  std::cout << " eps: " << epsilon << std::endl;
1178  std::cout << " numeric: float" << std::endl;
1179  std::cout << " layout: column-major" << std::endl;
1180  retval = test<NumericT, viennacl::column_major>(epsilon);
1181  if ( retval == EXIT_SUCCESS )
1182  std::cout << "# Test passed" << std::endl;
1183  else
1184  return retval;
1185  }
1186  std::cout << std::endl;
1187  std::cout << "----------------------------------------------" << std::endl;
1188  std::cout << std::endl;
1189 
1190 
1191 #ifdef VIENNACL_WITH_OPENCL
1192  if ( viennacl::ocl::current_device().double_support() )
1193 #endif
1194  {
1195  {
1196  typedef double NumericT;
1197  NumericT epsilon = 1.0E-11;
1198  std::cout << "# Testing setup:" << std::endl;
1199  std::cout << " eps: " << epsilon << std::endl;
1200  std::cout << " numeric: double" << std::endl;
1201  std::cout << " layout: row-major" << std::endl;
1202  retval = test<NumericT, viennacl::row_major>(epsilon);
1203  if ( retval == EXIT_SUCCESS )
1204  std::cout << "# Test passed" << std::endl;
1205  else
1206  return retval;
1207  }
1208  std::cout << std::endl;
1209  std::cout << "----------------------------------------------" << std::endl;
1210  std::cout << std::endl;
1211  {
1212  typedef double NumericT;
1213  NumericT epsilon = 1.0E-11;
1214  std::cout << "# Testing setup:" << std::endl;
1215  std::cout << " eps: " << epsilon << std::endl;
1216  std::cout << " numeric: double" << std::endl;
1217  std::cout << " layout: column-major" << std::endl;
1218  retval = test<NumericT, viennacl::column_major>(epsilon);
1219  if ( retval == EXIT_SUCCESS )
1220  std::cout << "# Test passed" << std::endl;
1221  else
1222  return retval;
1223  }
1224  std::cout << std::endl;
1225  std::cout << "----------------------------------------------" << std::endl;
1226  std::cout << std::endl;
1227  }
1228 
1229  std::cout << std::endl;
1230  std::cout << "------- Test completed --------" << std::endl;
1231  std::cout << std::endl;
1232 
1233 
1234  return retval;
1235 }
viennacl::linalg::cuda::inplace_solve
void inplace_solve(matrix_base< NumericT > const &A, matrix_base< NumericT > &B, SolverTagT tag)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:253
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
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::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
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:1142
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
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
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:885
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:900
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
prod
void prod(std::vector< std::map< IndexT, NumericT > > const &stl_A, std::vector< std::map< IndexT, NumericT > > const &stl_B, std::vector< std::map< IndexT, NumericT > > &stl_C)
Definition: sparse_prod.cpp:114
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
test_solve
int test_solve(Epsilon const &epsilon, UblasMatrixType &ublas_m1, UblasVectorType &ublas_v1, VCLMatrixType &vcl_m1, VCLVectorType1 &vcl_v1)
Definition: matrix_vector.cpp:374
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.
test_prod_rank1
int test_prod_rank1(Epsilon const &epsilon, UblasMatrixType &ublas_m1, UblasVectorType &ublas_v1, UblasVectorType &ublas_v2, UblasMatrixType &ublas_m2, VCLMatrixType &vcl_m1, VCLVectorType1 &vcl_v1, VCLVectorType2 &vcl_v2, VCLMatrixType &vcl_m2)
Definition: matrix_vector.cpp:117
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
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: matrix_vector.cpp:63
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:908
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:494
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