ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
blas3_solve.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 //#define NDEBUG
24 //#define VIENNACL_DEBUG_BUILD
25 
26 //
27 // *** System
28 //
29 #include <iostream>
30 
31 // We don't need debug mode in UBLAS:
32 #ifndef NDEBUG
33  #define BOOST_UBLAS_NDEBUG
34 #endif
35 
36 //
37 // *** Boost
38 //
39 #include <boost/numeric/ublas/io.hpp>
40 #include <boost/numeric/ublas/triangular.hpp>
41 #include <boost/numeric/ublas/matrix_sparse.hpp>
42 #include <boost/numeric/ublas/matrix.hpp>
43 #include <boost/numeric/ublas/matrix_proxy.hpp>
44 #include <boost/numeric/ublas/lu.hpp>
45 #include <boost/numeric/ublas/io.hpp>
46 
47 //
48 // *** ViennaCL
49 //
50 //#define VIENNACL_DEBUG_ALL
51 //#define VIENNACL_DEBUG_BUILD
52 #define VIENNACL_WITH_UBLAS 1
53 #include "viennacl/scalar.hpp"
54 #include "viennacl/matrix.hpp"
55 #include "viennacl/matrix_proxy.hpp"
56 #include "viennacl/vector.hpp"
57 #include "viennacl/linalg/prod.hpp"
58 #include "viennacl/linalg/norm_2.hpp"
59 #include "viennacl/linalg/direct_solve.hpp"
60 #include "viennacl/tools/random.hpp"
61 //
62 // -------------------------------------------------------------
63 //
64 using namespace boost::numeric;
65 //
66 // -------------------------------------------------------------
67 //
68 template<typename ScalarType>
69 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
70 {
71  viennacl::backend::finish();
72  if (s1 != s2)
73  return (s1 - s2) / std::max(fabs(s1), fabs(s2));
74  return 0;
75 }
76 
77 template<typename ScalarType>
78 ScalarType diff(ublas::vector<ScalarType> & v1, viennacl::vector<ScalarType> & v2)
79 {
80  ublas::vector<ScalarType> v2_cpu(v2.size());
81  viennacl::backend::finish();
82  viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
83  viennacl::backend::finish();
84 
85  for (std::size_t i=0;i<v1.size(); ++i)
86  {
87  if ( std::max( fabs(v2_cpu[i]), fabs(v1[i]) ) > 0 )
88  v2_cpu[i] = fabs(v2_cpu[i] - v1[i]) / std::max( fabs(v2_cpu[i]), fabs(v1[i]) );
89  else
90  v2_cpu[i] = 0.0;
91  }
92 
93  return norm_inf(v2_cpu);
94 }
95 
96 
97 template<typename ScalarType, typename VCLMatrixType>
98 ScalarType diff(ublas::matrix<ScalarType> & mat1, VCLMatrixType & mat2)
99 {
100  ublas::matrix<ScalarType> mat2_cpu(mat2.size1(), mat2.size2());
101  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
102  viennacl::copy(mat2, mat2_cpu);
103  ScalarType ret = 0;
104  ScalarType act = 0;
105 
106  for (unsigned int i = 0; i < mat2_cpu.size1(); ++i)
107  {
108  for (unsigned int j = 0; j < mat2_cpu.size2(); ++j)
109  {
110  act = std::fabs(mat2_cpu(i,j) - mat1(i,j)) / std::max( std::fabs(mat2_cpu(i, j)), std::fabs(mat1(i,j)) );
111  if (act > ret)
112  ret = act;
113  }
114  }
115  //std::cout << ret << std::endl;
116  return ret;
117 }
118 
119 
120 
121 //
122 // Triangular solvers
123 //
124 
125 
126 
127 template<typename RHSTypeRef, typename RHSTypeCheck, typename Epsilon >
128 void run_solver_check(RHSTypeRef & B_ref, RHSTypeCheck & B_check, int & retval, Epsilon const & epsilon)
129 {
130  double act_diff = fabs(diff(B_ref, B_check));
131  if ( act_diff > epsilon )
132  {
133  std::cout << " FAILED!" << std::endl;
134  std::cout << "# Error at operation: matrix-matrix solve" << std::endl;
135  std::cout << " diff: " << act_diff << std::endl;
136  retval = EXIT_FAILURE;
137  }
138  else
139  std::cout << " passed! " << act_diff << std::endl;
140 
141 }
142 
143 
144 template< typename NumericT, typename Epsilon,
145  typename ReferenceMatrixTypeA, typename ReferenceMatrixTypeB, typename ReferenceMatrixTypeC,
146  typename MatrixTypeA, typename MatrixTypeB, typename MatrixTypeC, typename MatrixTypeResult>
147 int test_solve(Epsilon const& epsilon,
148 
149  ReferenceMatrixTypeA const & A,
150  ReferenceMatrixTypeB const & B_start,
151  ReferenceMatrixTypeC const & C_start,
152 
153  MatrixTypeA const & vcl_A,
154  MatrixTypeB & vcl_B,
155  MatrixTypeC & vcl_C,
156  MatrixTypeResult const &
157  )
158 {
159  int retval = EXIT_SUCCESS;
160 
161  // --------------------------------------------------------------------------
162 
163  ReferenceMatrixTypeA result;
164  ReferenceMatrixTypeC C_trans;
165 
166  ReferenceMatrixTypeB B = B_start;
167  ReferenceMatrixTypeC C = C_start;
168 
169  MatrixTypeResult vcl_result;
170 
171  // Test: A \ B with various tags --------------------------------------------------------------------------
172  std::cout << "Testing A \\ B: " << std::endl;
173  std::cout << " * upper_tag: ";
174  result = ublas::solve(A, B, ublas::upper_tag());
175  vcl_result = viennacl::linalg::solve(vcl_A, vcl_B, viennacl::linalg::upper_tag());
176  run_solver_check(result, vcl_result, retval, epsilon);
177 
178  std::cout << " * unit_upper_tag: ";
179  result = ublas::solve(A, B, ublas::unit_upper_tag());
180  vcl_result = viennacl::linalg::solve(vcl_A, vcl_B, viennacl::linalg::unit_upper_tag());
181  run_solver_check(result, vcl_result, retval, epsilon);
182 
183  std::cout << " * lower_tag: ";
184  result = ublas::solve(A, B, ublas::lower_tag());
185  vcl_result = viennacl::linalg::solve(vcl_A, vcl_B, viennacl::linalg::lower_tag());
186  run_solver_check(result, vcl_result, retval, epsilon);
187 
188  std::cout << " * unit_lower_tag: ";
189  result = ublas::solve(A, B, ublas::unit_lower_tag());
190  vcl_result = viennacl::linalg::solve(vcl_A, vcl_B, viennacl::linalg::unit_lower_tag());
191  run_solver_check(result, vcl_result, retval, epsilon);
192 
193  if (retval == EXIT_SUCCESS)
194  std::cout << "Test A \\ B passed!" << std::endl;
195 
196  B = B_start;
197  C = C_start;
198 
199  // Test: A \ B^T --------------------------------------------------------------------------
200  std::cout << "Testing A \\ B^T: " << std::endl;
201  std::cout << " * upper_tag: ";
202  viennacl::copy(C, vcl_C); C_trans = trans(C);
203  //check solve():
204  result = ublas::solve(A, C_trans, ublas::upper_tag());
205  vcl_result = viennacl::linalg::solve(vcl_A, trans(vcl_C), viennacl::linalg::upper_tag());
206  run_solver_check(result, vcl_result, retval, epsilon);
207  //check compute kernels:
208  std::cout << " * upper_tag: ";
209  ublas::inplace_solve(A, C_trans, ublas::upper_tag());
210  viennacl::linalg::inplace_solve(vcl_A, trans(vcl_C), viennacl::linalg::upper_tag());
211  C = trans(C_trans); run_solver_check(C, vcl_C, retval, epsilon);
212 
213  std::cout << " * unit_upper_tag: ";
214  viennacl::copy(C, vcl_C); C_trans = trans(C);
215  ublas::inplace_solve(A, C_trans, ublas::unit_upper_tag());
216  viennacl::linalg::inplace_solve(vcl_A, trans(vcl_C), viennacl::linalg::unit_upper_tag());
217  C = trans(C_trans); run_solver_check(C, vcl_C, retval, epsilon);
218 
219  std::cout << " * lower_tag: ";
220  viennacl::copy(C, vcl_C); C_trans = trans(C);
221  ublas::inplace_solve(A, C_trans, ublas::lower_tag());
222  viennacl::linalg::inplace_solve(vcl_A, trans(vcl_C), viennacl::linalg::lower_tag());
223  C = trans(C_trans); run_solver_check(C, vcl_C, retval, epsilon);
224 
225  std::cout << " * unit_lower_tag: ";
226  viennacl::copy(C, vcl_C); C_trans = trans(C);
227  ublas::inplace_solve(A, C_trans, ublas::unit_lower_tag());
228  viennacl::linalg::inplace_solve(vcl_A, trans(vcl_C), viennacl::linalg::unit_lower_tag());
229  C = trans(C_trans); run_solver_check(C, vcl_C, retval, epsilon);
230 
231  if (retval == EXIT_SUCCESS)
232  std::cout << "Test A \\ B^T passed!" << std::endl;
233 
234  B = B_start;
235  C = C_start;
236 
237  // Test: A \ B with various tags --------------------------------------------------------------------------
238  std::cout << "Testing A^T \\ B: " << std::endl;
239  std::cout << " * upper_tag: ";
240  viennacl::copy(B, vcl_B);
241  result = ublas::solve(trans(A), B, ublas::upper_tag());
242  vcl_result = viennacl::linalg::solve(trans(vcl_A), vcl_B, viennacl::linalg::upper_tag());
243  run_solver_check(result, vcl_result, retval, epsilon);
244 
245  std::cout << " * unit_upper_tag: ";
246  viennacl::copy(B, vcl_B);
247  result = ublas::solve(trans(A), B, ublas::unit_upper_tag());
248  vcl_result = viennacl::linalg::solve(trans(vcl_A), vcl_B, viennacl::linalg::unit_upper_tag());
249  run_solver_check(result, vcl_result, retval, epsilon);
250 
251  std::cout << " * lower_tag: ";
252  viennacl::copy(B, vcl_B);
253  result = ublas::solve(trans(A), B, ublas::lower_tag());
254  vcl_result = viennacl::linalg::solve(trans(vcl_A), vcl_B, viennacl::linalg::lower_tag());
255  run_solver_check(result, vcl_result, retval, epsilon);
256 
257  std::cout << " * unit_lower_tag: ";
258  viennacl::copy(B, vcl_B);
259  result = ublas::solve(trans(A), B, ublas::unit_lower_tag());
260  vcl_result = viennacl::linalg::solve(trans(vcl_A), vcl_B, viennacl::linalg::unit_lower_tag());
261  run_solver_check(result, vcl_result, retval, epsilon);
262 
263  if (retval == EXIT_SUCCESS)
264  std::cout << "Test A^T \\ B passed!" << std::endl;
265 
266  B = B_start;
267  C = C_start;
268 
269  // Test: A^T \ B^T --------------------------------------------------------------------------
270  std::cout << "Testing A^T \\ B^T: " << std::endl;
271  std::cout << " * upper_tag: ";
272  viennacl::copy(C, vcl_C); C_trans = trans(C);
273  //check solve():
274  result = ublas::solve(trans(A), C_trans, ublas::upper_tag());
275  vcl_result = viennacl::linalg::solve(trans(vcl_A), trans(vcl_C), viennacl::linalg::upper_tag());
276  run_solver_check(result, vcl_result, retval, epsilon);
277  //check kernels:
278  std::cout << " * upper_tag: ";
279  ublas::inplace_solve(trans(A), C_trans, ublas::upper_tag());
280  viennacl::linalg::inplace_solve(trans(vcl_A), trans(vcl_C), viennacl::linalg::upper_tag());
281  C = trans(C_trans); run_solver_check(C, vcl_C, retval, epsilon);
282 
283  std::cout << " * unit_upper_tag: ";
284  viennacl::copy(C, vcl_C); C_trans = trans(C);
285  ublas::inplace_solve(trans(A), C_trans, ublas::unit_upper_tag());
286  viennacl::linalg::inplace_solve(trans(vcl_A), trans(vcl_C), viennacl::linalg::unit_upper_tag());
287  C = trans(C_trans); run_solver_check(C, vcl_C, retval, epsilon);
288 
289  std::cout << " * lower_tag: ";
290  viennacl::copy(C, vcl_C); C_trans = trans(C);
291  ublas::inplace_solve(trans(A), C_trans, ublas::lower_tag());
292  viennacl::linalg::inplace_solve(trans(vcl_A), trans(vcl_C), viennacl::linalg::lower_tag());
293  C = trans(C_trans); run_solver_check(C, vcl_C, retval, epsilon);
294 
295  std::cout << " * unit_lower_tag: ";
296  viennacl::copy(C, vcl_C); C_trans = trans(C);
297  ublas::inplace_solve(trans(A), C_trans, ublas::unit_lower_tag());
298  viennacl::linalg::inplace_solve(trans(vcl_A), trans(vcl_C), viennacl::linalg::unit_lower_tag());
299  C = trans(C_trans); run_solver_check(C, vcl_C, retval, epsilon);
300 
301  if (retval == EXIT_SUCCESS)
302  std::cout << "Test A^T \\ B^T passed!" << std::endl;
303 
304  return retval;
305 }
306 
307 
308 template< typename NumericT, typename F_A, typename F_B, typename Epsilon >
309 int test_solve(Epsilon const& epsilon)
310 {
311  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
312 
313  int ret = EXIT_SUCCESS;
314  std::size_t matrix_size = 135; //some odd number, not too large
315  std::size_t rhs_num = 67;
316 
317  std::cout << "--- Part 2: Testing matrix-matrix solver ---" << std::endl;
318 
319 
320  ublas::matrix<NumericT> A(matrix_size, matrix_size);
321  ublas::matrix<NumericT> B_start(matrix_size, rhs_num);
322  ublas::matrix<NumericT> C_start(rhs_num, matrix_size);
323 
324  for (std::size_t i = 0; i < A.size1(); ++i)
325  {
326  for (std::size_t j = 0; j < A.size2(); ++j)
327  A(i,j) = static_cast<NumericT>(-0.5) * randomNumber();
328  A(i,i) = NumericT(1.0) + NumericT(2.0) * randomNumber(); //some extra weight on diagonal for stability
329  }
330 
331  for (std::size_t i = 0; i < B_start.size1(); ++i)
332  for (std::size_t j = 0; j < B_start.size2(); ++j)
333  B_start(i,j) = randomNumber();
334 
335  for (std::size_t i = 0; i < C_start.size1(); ++i)
336  for (std::size_t j = 0; j < C_start.size2(); ++j)
337  C_start(i,j) = randomNumber();
338 
339 
340  // A
341  viennacl::range range1_A(matrix_size, 2*matrix_size);
342  viennacl::range range2_A(2*matrix_size, 3*matrix_size);
343  viennacl::slice slice1_A(matrix_size, 2, matrix_size);
344  viennacl::slice slice2_A(0, 3, matrix_size);
345 
346  viennacl::matrix<NumericT, F_A> vcl_A(matrix_size, matrix_size);
347  viennacl::copy(A, vcl_A);
348 
349  viennacl::matrix<NumericT, F_A> vcl_big_range_A(4*matrix_size, 4*matrix_size);
350  viennacl::matrix_range<viennacl::matrix<NumericT, F_A> > vcl_range_A(vcl_big_range_A, range1_A, range2_A);
351  viennacl::copy(A, vcl_range_A);
352 
353  viennacl::matrix<NumericT, F_A> vcl_big_slice_A(4*matrix_size, 4*matrix_size);
354  viennacl::matrix_slice<viennacl::matrix<NumericT, F_A> > vcl_slice_A(vcl_big_slice_A, slice1_A, slice2_A);
355  viennacl::copy(A, vcl_slice_A);
356 
357 
358  // B
359  viennacl::range range1_B(matrix_size, 2*matrix_size);
360  viennacl::range range2_B(2*rhs_num, 3*rhs_num);
361  viennacl::slice slice1_B(matrix_size, 2, matrix_size);
362  viennacl::slice slice2_B(0, 3, rhs_num);
363 
364  viennacl::matrix<NumericT, F_B> vcl_B(matrix_size, rhs_num);
365  viennacl::copy(B_start, vcl_B);
366 
367  viennacl::matrix<NumericT, F_B> vcl_big_range_B(4*matrix_size, 4*rhs_num);
368  viennacl::matrix_range<viennacl::matrix<NumericT, F_B> > vcl_range_B(vcl_big_range_B, range1_B, range2_B);
369  viennacl::copy(B_start, vcl_range_B);
370 
371  viennacl::matrix<NumericT, F_B> vcl_big_slice_B(4*matrix_size, 4*rhs_num);
372  viennacl::matrix_slice<viennacl::matrix<NumericT, F_B> > vcl_slice_B(vcl_big_slice_B, slice1_B, slice2_B);
373  viennacl::copy(B_start, vcl_slice_B);
374 
375 
376  // C
377  viennacl::range range1_C(rhs_num, 2*rhs_num);
378  viennacl::range range2_C(2*matrix_size, 3*matrix_size);
379  viennacl::slice slice1_C(rhs_num, 2, rhs_num);
380  viennacl::slice slice2_C(0, 3, matrix_size);
381 
382  viennacl::matrix<NumericT, F_B> vcl_C(rhs_num, matrix_size);
383  viennacl::copy(C_start, vcl_C);
384 
385  viennacl::matrix<NumericT, F_B> vcl_big_range_C(4*rhs_num, 4*matrix_size);
386  viennacl::matrix_range<viennacl::matrix<NumericT, F_B> > vcl_range_C(vcl_big_range_C, range1_C, range2_C);
387  viennacl::copy(C_start, vcl_range_C);
388 
389  viennacl::matrix<NumericT, F_B> vcl_big_slice_C(4*rhs_num, 4*matrix_size);
390  viennacl::matrix_slice<viennacl::matrix<NumericT, F_B> > vcl_slice_C(vcl_big_slice_C, slice1_C, slice2_C);
391  viennacl::copy(C_start, vcl_slice_C);
392 
393 
394  std::cout << "Now using A=matrix, B=matrix" << std::endl;
395  ret = test_solve<NumericT>(epsilon,
396  A, B_start, C_start,
397  vcl_A, vcl_B, vcl_C, vcl_B
398  );
399  if (ret != EXIT_SUCCESS)
400  return ret;
401 
402  std::cout << "Now using A=matrix, B=range" << std::endl;
403  ret = test_solve<NumericT>(epsilon,
404  A, B_start, C_start,
405  vcl_A, vcl_range_B, vcl_range_C, vcl_B
406  );
407  if (ret != EXIT_SUCCESS)
408  return ret;
409 
410  std::cout << "Now using A=matrix, B=slice" << std::endl;
411  ret = test_solve<NumericT>(epsilon,
412  A, B_start, C_start,
413  vcl_A, vcl_slice_B, vcl_slice_C, vcl_B
414  );
415  if (ret != EXIT_SUCCESS)
416  return ret;
417 
418 
419 
420  std::cout << "Now using A=range, B=matrix" << std::endl;
421  ret = test_solve<NumericT>(epsilon,
422  A, B_start, C_start,
423  vcl_range_A, vcl_B, vcl_C, vcl_B
424  );
425  if (ret != EXIT_SUCCESS)
426  return ret;
427 
428  std::cout << "Now using A=range, B=range" << std::endl;
429  ret = test_solve<NumericT>(epsilon,
430  A, B_start, C_start,
431  vcl_range_A, vcl_range_B, vcl_range_C, vcl_B
432  );
433  if (ret != EXIT_SUCCESS)
434  return ret;
435 
436  std::cout << "Now using A=range, B=slice" << std::endl;
437  ret = test_solve<NumericT>(epsilon,
438  A, B_start, C_start,
439  vcl_range_A, vcl_slice_B, vcl_slice_C, vcl_B
440  );
441  if (ret != EXIT_SUCCESS)
442  return ret;
443 
444 
445 
446 
447  std::cout << "Now using A=slice, B=matrix" << std::endl;
448  ret = test_solve<NumericT>(epsilon,
449  A, B_start, C_start,
450  vcl_slice_A, vcl_B, vcl_C, vcl_B
451  );
452  if (ret != EXIT_SUCCESS)
453  return ret;
454 
455  std::cout << "Now using A=slice, B=range" << std::endl;
456  ret = test_solve<NumericT>(epsilon,
457  A, B_start, C_start,
458  vcl_slice_A, vcl_range_B, vcl_range_C, vcl_B
459  );
460  if (ret != EXIT_SUCCESS)
461  return ret;
462 
463  std::cout << "Now using A=slice, B=slice" << std::endl;
464  ret = test_solve<NumericT>(epsilon,
465  A, B_start, C_start,
466  vcl_slice_A, vcl_slice_B, vcl_slice_C, vcl_B
467  );
468  if (ret != EXIT_SUCCESS)
469  return ret;
470 
471 
472 
473 
474  return ret;
475 
476 }
477 
478 
479 
480 //
481 // Control functions
482 //
483 
484 
485 template< typename NumericT, typename Epsilon >
486 int test(Epsilon const& epsilon)
487 {
488  int ret;
489 
490  std::cout << "////////////////////////////////" << std::endl;
491  std::cout << "/// Now testing A=row, B=row ///" << std::endl;
492  std::cout << "////////////////////////////////" << std::endl;
493  ret = test_solve<NumericT, viennacl::row_major, viennacl::row_major>(epsilon);
494  if (ret != EXIT_SUCCESS)
495  return ret;
496 
497 
498  std::cout << "////////////////////////////////" << std::endl;
499  std::cout << "/// Now testing A=row, B=col ///" << std::endl;
500  std::cout << "////////////////////////////////" << std::endl;
501  ret = test_solve<NumericT, viennacl::row_major, viennacl::column_major>(epsilon);
502  if (ret != EXIT_SUCCESS)
503  return ret;
504 
505  std::cout << "////////////////////////////////" << std::endl;
506  std::cout << "/// Now testing A=col, B=row ///" << std::endl;
507  std::cout << "////////////////////////////////" << std::endl;
508  ret = test_solve<NumericT, viennacl::column_major, viennacl::row_major>(epsilon);
509  if (ret != EXIT_SUCCESS)
510  return ret;
511 
512  std::cout << "////////////////////////////////" << std::endl;
513  std::cout << "/// Now testing A=col, B=col ///" << std::endl;
514  std::cout << "////////////////////////////////" << std::endl;
515  ret = test_solve<NumericT, viennacl::column_major, viennacl::column_major>(epsilon);
516  if (ret != EXIT_SUCCESS)
517  return ret;
518 
519 
520 
521  return ret;
522 }
523 
524 //
525 // -------------------------------------------------------------
526 //
527 int main()
528 {
529  std::cout << std::endl;
530  std::cout << "----------------------------------------------" << std::endl;
531  std::cout << "----------------------------------------------" << std::endl;
532  std::cout << "## Test :: BLAS 3 routines" << std::endl;
533  std::cout << "----------------------------------------------" << std::endl;
534  std::cout << "----------------------------------------------" << std::endl;
535  std::cout << std::endl;
536 
537  int retval = EXIT_SUCCESS;
538 
539  std::cout << std::endl;
540  std::cout << "----------------------------------------------" << std::endl;
541  std::cout << std::endl;
542  {
543  typedef float NumericT;
544  NumericT epsilon = NumericT(1.0E-3);
545  std::cout << "# Testing setup:" << std::endl;
546  std::cout << " eps: " << epsilon << std::endl;
547  std::cout << " numeric: float" << std::endl;
548  retval = test<NumericT>(epsilon);
549  if ( retval == EXIT_SUCCESS )
550  std::cout << "# Test passed" << std::endl;
551  else
552  return retval;
553  }
554  std::cout << std::endl;
555  std::cout << "----------------------------------------------" << std::endl;
556  std::cout << std::endl;
557 #ifdef VIENNACL_WITH_OPENCL
558  if ( viennacl::ocl::current_device().double_support() )
559 #endif
560  {
561  {
562  typedef double NumericT;
563  NumericT epsilon = 1.0E-11;
564  std::cout << "# Testing setup:" << std::endl;
565  std::cout << " eps: " << epsilon << std::endl;
566  std::cout << " numeric: double" << std::endl;
567  retval = test<NumericT>(epsilon);
568  if ( retval == EXIT_SUCCESS )
569  std::cout << "# Test passed" << std::endl;
570  else
571  return retval;
572  }
573  std::cout << std::endl;
574  std::cout << "----------------------------------------------" << std::endl;
575  std::cout << std::endl;
576  }
577 
578  std::cout << std::endl;
579  std::cout << "------- Test completed --------" << std::endl;
580  std::cout << std::endl;
581 
582 
583  return retval;
584 }
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::linalg::inplace_solve
void inplace_solve(const matrix_base< NumericT > &A, matrix_base< NumericT > &B, SolverTagT)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:217
main
int main()
Definition: blas3_solve.cpp:527
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::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
NumericT
float NumericT
Definition: bisect.cpp:40
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
viennacl::tools::uniform_random_numbers
Random number generator for returning uniformly distributed values in the closed interval [0...
Definition: random.hpp:44
viennacl::vector_base::begin
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
Definition: vector.hpp:841
viennacl::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
viennacl::ocl::device::double_support
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Definition: device.hpp:956
direct_solve.hpp
Implementations of dense direct solvers are found here.
matrix_proxy.hpp
Proxy classes for matrices.
test_solve
int test_solve(Epsilon const &epsilon, ReferenceMatrixTypeA const &A, ReferenceMatrixTypeB const &B_start, ReferenceMatrixTypeC const &C_start, MatrixTypeA const &vcl_A, MatrixTypeB &vcl_B, MatrixTypeC &vcl_C, MatrixTypeResult const &)
Definition: blas3_solve.cpp:147
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::vector_base::size
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
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
int test(Epsilon const &epsilon)
Definition: blas3_solve.cpp:486
viennacl::matrix_range
Class for representing non-strided submatrices of a bigger matrix A.
Definition: forwards.h:440
run_solver_check
void run_solver_check(RHSTypeRef &B_ref, RHSTypeCheck &B_check, int &retval, Epsilon const &epsilon)
Definition: blas3_solve.cpp:128
viennacl::vector_base::end
iterator end()
Returns an iterator pointing to the end of the vector (STL like)
Definition: vector.hpp:848
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
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: blas3_solve.cpp:69
scalar.hpp
Implementation of the ViennaCL scalar class.
viennacl::linalg::unit_upper_tag
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:864