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
19 /** \file tests/src/blas3_solve.cpp Tests the BLAS level 3 triangular solvers
20 * \test Tests the BLAS level 3 triangular solvers
21 **/
22
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>
diff(ScalarType & s1,viennacl::scalar<ScalarType> & s2)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>
diff(std::vector<ScalarType> & v1,viennacl::vector<ScalarType> & v2)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>
inplace_solve_lower(std::vector<std::vector<NumericT>> const & A,std::vector<std::vector<NumericT>> & B,bool unit_diagonal)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>
inplace_solve(std::vector<std::vector<NumericT>> const & A,std::vector<std::vector<NumericT>> & B,viennacl::linalg::lower_tag)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>
inplace_solve(std::vector<std::vector<NumericT>> const & A,std::vector<std::vector<NumericT>> & B,viennacl::linalg::unit_lower_tag)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>
inplace_solve_upper(std::vector<std::vector<NumericT>> const & A,std::vector<std::vector<NumericT>> & B,bool unit_diagonal)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>
inplace_solve(std::vector<std::vector<NumericT>> const & A,std::vector<std::vector<NumericT>> & B,viennacl::linalg::upper_tag)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>
inplace_solve(std::vector<std::vector<NumericT>> const & A,std::vector<std::vector<NumericT>> & B,viennacl::linalg::unit_upper_tag)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>
solve(std::vector<std::vector<NumericT>> const & A,std::vector<std::vector<NumericT>> const & B,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 >
run_solver_check(RHSTypeRef & B_ref,RHSTypeCheck & B_check,int & retval,Epsilon const & 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>
trans(std::vector<std::vector<NumericT>> const & A)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>
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 &)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 >
test_solve(Epsilon const & 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 >
test(Epsilon const & 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 //
main()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 }
648