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