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 
20 /** \file tests/src/sparse.cpp  Tests sparse matrix operations.
21 *   \test  Tests sparse matrix operations.
22 **/
23 
24 //
25 // *** System
26 //
27 #include <iostream>
28 #include <vector>
29 #include <map>
30 #include <cmath>
31 
32 //
33 // *** ViennaCL
34 //
35 #include "viennacl/scalar.hpp"
36 #include "viennacl/compressed_matrix.hpp"
37 #include "viennacl/compressed_compressed_matrix.hpp"
38 #include "viennacl/coordinate_matrix.hpp"
39 #include "viennacl/ell_matrix.hpp"
40 #include "viennacl/sliced_ell_matrix.hpp"
41 #include "viennacl/hyb_matrix.hpp"
42 #include "viennacl/vector.hpp"
43 #include "viennacl/vector_proxy.hpp"
44 #include "viennacl/linalg/prod.hpp"
45 #include "viennacl/linalg/norm_2.hpp"
46 #include "viennacl/linalg/ilu.hpp"
47 #include "viennacl/linalg/detail/ilu/common.hpp"
48 #include "viennacl/io/matrix_market.hpp"
49 #include "viennacl/tools/random.hpp"
50 
51 
52 
53 //
54 // -------------------------------------------------------------
55 //
56 template<typename ScalarType>
diff(ScalarType & s1,viennacl::scalar<ScalarType> & s2)57 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
58 {
59    if (s1 != s2)
60       return (s1 - s2) / std::max(fabs(s1), std::fabs(s2));
61    return 0;
62 }
63 
64 template<typename ScalarType>
diff(std::vector<ScalarType> & v1,viennacl::vector<ScalarType> & v2)65 ScalarType diff(std::vector<ScalarType> & v1, viennacl::vector<ScalarType> & v2)
66 {
67    std::vector<ScalarType> v2_cpu(v2.size());
68    viennacl::backend::finish();
69    viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
70 
71    for (unsigned int i=0;i<v1.size(); ++i)
72    {
73       if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
74       {
75         //if (std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) < 1e-10 )  //absolute tolerance (avoid round-off issues)
76         //  v2_cpu[i] = 0;
77         //else
78           v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
79       }
80       else
81          v2_cpu[i] = 0.0;
82 
83       if (v2_cpu[i] > 0.0001)
84       {
85         //std::cout << "Neighbor: "      << i-1 << ": " << v1[i-1] << " vs. " << v2_cpu[i-1] << std::endl;
86         std::cout << "Error at entry " << i   << ": Should: " << v1[i]   << " vs. Is: " << v2[i]   << std::endl;
87         //std::cout << "Neighbor: "      << i+1 << ": " << v1[i+1] << " vs. " << v2_cpu[i+1] << std::endl;
88         exit(EXIT_FAILURE);
89       }
90    }
91 
92    ScalarType norm_inf = 0;
93    for (std::size_t i=0; i<v2_cpu.size(); ++i)
94      norm_inf = std::max<ScalarType>(norm_inf, std::fabs(v2_cpu[i]));
95 
96    return norm_inf;
97 }
98 
99 
100 template<typename IndexT, typename NumericT, typename SparseMatrixT>
101 NumericT diff(std::vector<std::map<IndexT, NumericT> > & cpu_A, SparseMatrixT & vcl_A)
102 {
103   typedef typename std::map<IndexT, NumericT>::const_iterator  RowIterator;
104 
105   std::vector<std::map<IndexT, NumericT> > from_gpu(vcl_A.size1());
106 
107   viennacl::backend::finish();
108   viennacl::copy(vcl_A, from_gpu);
109 
110   NumericT error = 0;
111 
112   //step 1: compare all entries from cpu_A with vcl_A:
113   for (std::size_t i=0; i<cpu_A.size(); ++i)
114   {
115     //std::cout << "Row " << row_it.index1() << ": " << std::endl;
116     for (RowIterator it = cpu_A[i].begin(); it != cpu_A[i].end(); ++it)
117     {
118       //std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
119       NumericT current_error = 0;
120       NumericT val_cpu_A = it->second;
121       NumericT val_gpu_A = from_gpu[i][it->first];
122 
123       NumericT max_val = std::max(std::fabs(val_cpu_A), std::fabs(val_gpu_A));
124       if (max_val > 0)
125         current_error = std::fabs(val_cpu_A - val_gpu_A) / max_val;
126       if (current_error > error)
127         error = current_error;
128     }
129   }
130 
131   //step 2: compare all entries from gpu_matrix with cpu_matrix (sparsity pattern might differ):
132   //std::cout << "ViennaCL matrix: " << std::endl;
133   for (std::size_t i=0; i<from_gpu.size(); ++i)
134   {
135     //std::cout << "Row " << row_it.index1() << ": " << std::endl;
136     for (RowIterator it = from_gpu[i].begin(); it != from_gpu[i].end(); ++it)
137     {
138       //std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
139       NumericT current_error = 0;
140       NumericT val_gpu_A = it->second;
141       NumericT val_cpu_A = cpu_A[i][it->first];
142 
143       NumericT max_val = std::max(std::fabs(val_cpu_A), std::fabs(val_gpu_A));
144       if (max_val > 0)
145         current_error = std::fabs(val_cpu_A - val_gpu_A) / max_val;
146       if (current_error > error)
147         error = current_error;
148     }
149   }
150 
151   return error;
152 }
153 
154 
155 template<typename NumericT, typename VCL_MatrixT, typename Epsilon, typename STLVectorT, typename VCLVectorT>
strided_matrix_vector_product_test(Epsilon epsilon,STLVectorT & result,STLVectorT const & rhs,VCLVectorT & vcl_result,VCLVectorT & vcl_rhs)156 int strided_matrix_vector_product_test(Epsilon epsilon,
157                                        STLVectorT & result,     STLVectorT const & rhs,
158                                        VCLVectorT & vcl_result, VCLVectorT & vcl_rhs)
159 {
160   typedef typename std::map<unsigned int, NumericT>::const_iterator    RowIterator;
161     int retval = EXIT_SUCCESS;
162 
163     std::vector<std::map<unsigned int, NumericT> > std_A(5);
164     std_A[0][0] = NumericT(2.0); std_A[0][2] = NumericT(-1.0);
165     std_A[1][0] = NumericT(3.0); std_A[1][2] = NumericT(-5.0);
166     std_A[2][1] = NumericT(5.0); std_A[2][2] = NumericT(-2.0);
167     std_A[3][2] = NumericT(1.0); std_A[3][3] = NumericT(-6.0);
168     std_A[4][1] = NumericT(7.0); std_A[4][2] = NumericT(-5.0);
169     //the following computes project(result, slice(1, 3, 5)) = prod(std_A, project(rhs, slice(3, 2, 4)));
170     for (std::size_t i=0; i<5; ++i)
171     {
172       NumericT val = 0;
173       for (RowIterator it = std_A[i].begin(); it != std_A[i].end(); ++it)
174         val += it->second * rhs[3 + 2*it->first];
175       result[1 + 3*i] = val;
176     }
177 
178     VCL_MatrixT vcl_sparse_matrix2;
179     viennacl::copy(std_A, vcl_sparse_matrix2);
180     viennacl::vector<NumericT> vec(4);
181     vec(0) = rhs[3];
182     vec(1) = rhs[5];
183     vec(2) = rhs[7];
184     vec(3) = rhs[9];
185     viennacl::project(vcl_result, viennacl::slice(1, 3, 5)) = viennacl::linalg::prod(vcl_sparse_matrix2, viennacl::project(vcl_rhs, viennacl::slice(3, 2, 4)));
186 
187     if ( std::fabs(diff(result, vcl_result)) > epsilon )
188     {
189       std::cout << "# Error at operation: matrix-vector product with strided vectors, part 1" << std::endl;
190       std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
191       retval = EXIT_FAILURE;
192     }
193     vcl_result(1)  = NumericT(1.0);
194     vcl_result(4)  = NumericT(1.0);
195     vcl_result(7)  = NumericT(1.0);
196     vcl_result(10) = NumericT(1.0);
197     vcl_result(13) = NumericT(1.0);
198 
199     viennacl::project(vcl_result, viennacl::slice(1, 3, 5)) = viennacl::linalg::prod(vcl_sparse_matrix2, vec);
200 
201     if ( std::fabs(diff(result, vcl_result)) > epsilon )
202     {
203       std::cout << "# Error at operation: matrix-vector product with strided vectors, part 2" << std::endl;
204       std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
205       retval = EXIT_FAILURE;
206     }
207 
208     return retval;
209 }
210 
211 
212 template< typename NumericT, typename VCL_MATRIX, typename Epsilon >
resize_test(Epsilon const & epsilon)213 int resize_test(Epsilon const& epsilon)
214 {
215    int retval = EXIT_SUCCESS;
216 
217    std::vector<std::map<unsigned int, NumericT> > std_A(5);
218    VCL_MATRIX vcl_matrix;
219 
220    std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
221    std_A[1][0] = NumericT(1.0);  std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
222    std_A[2][0] = NumericT(2.0);  std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
223    std_A[3][0] = NumericT(3.0);  std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
224    std_A[4][0] = NumericT(4.0);  std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
225 
226    viennacl::copy(std_A, vcl_matrix);
227    std::vector<std::map<unsigned int, NumericT> > std_B(std_A.size());
228    viennacl::copy(vcl_matrix, std_B);
229 
230    std::cout << "Checking for equality after copy..." << std::endl;
231     if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
232     {
233         std::cout << "# Error at operation: equality after copy with sparse matrix" << std::endl;
234         std::cout << "  diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
235         return EXIT_FAILURE;
236     }
237 
238    std::cout << "Testing resize to larger..." << std::endl;
239    std_A.resize(10);
240    std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
241    std_A[1][0] = NumericT( 1.0); std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
242    std_A[2][0] = NumericT( 2.0); std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
243    std_A[3][0] = NumericT( 3.0); std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
244    std_A[4][0] = NumericT( 4.0); std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
245 
246    vcl_matrix.resize(10, 10, true);
247 
248     if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
249     {
250         std::cout << "# Error at operation: resize (to larger) with sparse matrix" << std::endl;
251         std::cout << "  diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
252         return EXIT_FAILURE;
253     }
254 
255    std_A[5][5] = NumericT(5.5); std_A[5][6] = NumericT(5.6); std_A[5][7] = NumericT(5.7); std_A[5][8] = NumericT(5.8); std_A[5][9] = NumericT(5.9);
256    std_A[6][5] = NumericT(6.5); std_A[6][6] = NumericT(6.6); std_A[6][7] = NumericT(6.7); std_A[6][8] = NumericT(6.8); std_A[6][9] = NumericT(6.9);
257    std_A[7][5] = NumericT(7.5); std_A[7][6] = NumericT(7.6); std_A[7][7] = NumericT(7.7); std_A[7][8] = NumericT(7.8); std_A[7][9] = NumericT(7.9);
258    std_A[8][5] = NumericT(8.5); std_A[8][6] = NumericT(8.6); std_A[8][7] = NumericT(8.7); std_A[8][8] = NumericT(8.8); std_A[8][9] = NumericT(8.9);
259    std_A[9][5] = NumericT(9.5); std_A[9][6] = NumericT(9.6); std_A[9][7] = NumericT(9.7); std_A[9][8] = NumericT(9.8); std_A[9][9] = NumericT(9.9);
260    viennacl::copy(std_A, vcl_matrix);
261 
262    std::cout << "Testing resize to smaller..." << std::endl;
263    std_A.clear();
264    std_A.resize(7);
265    std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
266    std_A[1][0] = NumericT( 1.0); std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
267    std_A[2][0] = NumericT( 2.0); std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
268    std_A[3][0] = NumericT( 3.0); std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
269    std_A[4][0] = NumericT( 4.0); std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
270    std_A[5][5] = NumericT( 5.5); std_A[5][6] = NumericT(5.6); //std_A[5][7] = NumericT(5.7); std_A[5][8] = NumericT(5.8); std_A[5][9] = NumericT(5.9);
271    std_A[6][5] = NumericT( 6.5); std_A[6][6] = NumericT(6.6); //std_A[6][7] = NumericT(6.7); std_A[6][8] = NumericT(6.8); std_A[6][9] = NumericT(6.9);
272 
273    vcl_matrix.resize(7, 7);
274 
275    //std::cout << std_A << std::endl;
276     if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
277     {
278         std::cout << "# Error at operation: resize (to smaller) with sparse matrix" << std::endl;
279         std::cout << "  diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
280         retval = EXIT_FAILURE;
281     }
282 
283    std::vector<NumericT> std_vec(std_A.size(), NumericT(3.1415));
284    viennacl::vector<NumericT> vcl_vec(std_A.size());
285 
286 
287   std::cout << "Testing unit lower triangular solve: compressed_matrix" << std::endl;
288   viennacl::copy(std_vec, vcl_vec);
289 
290   std::cout << "STL..." << std::endl;
291   //boost::numeric::ublas::inplace_solve((ublas_matrix), ublas_vec, boost::numeric::ublas::unit_lower_tag());
292   for (std::size_t i=1; i<std_A.size(); ++i)
293     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_A[i].begin(); it != std_A[i].end(); ++it)
294     {
295       if (it->first < static_cast<unsigned int>(i))
296         std_vec[i] -= it->second * std_vec[it->first];
297       else
298         continue;
299     }
300 
301   std::cout << "ViennaCL..." << std::endl;
302   viennacl::linalg::inplace_solve((vcl_matrix), vcl_vec, viennacl::linalg::unit_lower_tag());
303 
304   if ( std::fabs(diff(std_vec, vcl_vec)) > epsilon )
305   {
306       std::cout << "# Error at operation: unit lower triangular solve" << std::endl;
307       std::cout << "  diff: " << std::fabs(diff(std_vec, vcl_vec)) << std::endl;
308       retval = EXIT_FAILURE;
309   }
310   return retval;
311 }
312 
313 
314 //
315 // -------------------------------------------------------------
316 //
317 template< typename NumericT, typename Epsilon >
test(Epsilon const & epsilon)318 int test(Epsilon const& epsilon)
319 {
320   viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
321 
322   std::cout << "Testing resizing of compressed_matrix..." << std::endl;
323   int retval = resize_test<NumericT, viennacl::compressed_matrix<NumericT> >(epsilon);
324   if (retval != EXIT_SUCCESS)
325     return retval;
326 
327   // --------------------------------------------------------------------------
328   std::vector<NumericT> rhs;
329   std::vector<NumericT> result;
330   std::vector<std::map<unsigned int, NumericT> > std_matrix;
331 
332   if (viennacl::io::read_matrix_market_file(std_matrix, "../examples/testdata/mat65k.mtx") == EXIT_FAILURE)
333   {
334     std::cout << "Error reading Matrix file" << std::endl;
335     return EXIT_FAILURE;
336   }
337 
338   //unsigned int cg_mat_size = cg_mat.size();
339   std::cout << "done reading matrix" << std::endl;
340 
341 
342   rhs.resize(std_matrix.size());
343   for (std::size_t i=0; i<rhs.size(); ++i)
344   {
345     std_matrix[i][static_cast<unsigned int>(i)] = NumericT(0.5);   // Get rid of round-off errors by making row-sums unequal to zero:
346     rhs[i] = NumericT(1) + randomNumber();
347   }
348 
349   // add some random numbers to the double-compressed matrix:
350   std::vector<std::map<unsigned int, NumericT> > std_cc_matrix(std_matrix.size());
351   std_cc_matrix[42][199] = NumericT(3.1415);
352   std_cc_matrix[31][69] = NumericT(2.71);
353   std_cc_matrix[23][32] = NumericT(6);
354   std_cc_matrix[177][57] = NumericT(4);
355   std_cc_matrix[21][97] = NumericT(-4);
356   std_cc_matrix[92][25] = NumericT(2);
357   std_cc_matrix[89][62] = NumericT(11);
358   std_cc_matrix[ 1][ 7] = NumericT(8);
359   std_cc_matrix[85][41] = NumericT(13);
360   std_cc_matrix[66][28] = NumericT(8);
361   std_cc_matrix[21][74] = NumericT(-2);
362   viennacl::tools::sparse_matrix_adapter<NumericT> adapted_std_cc_matrix(std_cc_matrix, std_matrix.size(), std_matrix.size());
363 
364 
365   result = rhs;
366 
367 
368   viennacl::vector<NumericT> vcl_rhs(rhs.size());
369   viennacl::vector<NumericT> vcl_result(result.size());
370   viennacl::vector<NumericT> vcl_result2(result.size());
371   viennacl::compressed_matrix<NumericT> vcl_compressed_matrix(rhs.size(), rhs.size());
372   viennacl::compressed_compressed_matrix<NumericT> vcl_compressed_compressed_matrix(rhs.size(), rhs.size());
373   viennacl::coordinate_matrix<NumericT> vcl_coordinate_matrix(rhs.size(), rhs.size());
374   viennacl::ell_matrix<NumericT> vcl_ell_matrix;
375   viennacl::sliced_ell_matrix<NumericT> vcl_sliced_ell_matrix;
376   viennacl::hyb_matrix<NumericT> vcl_hyb_matrix;
377 
378   viennacl::copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
379   viennacl::copy(std_matrix, vcl_compressed_matrix);
380   viennacl::copy(adapted_std_cc_matrix, vcl_compressed_compressed_matrix);
381   viennacl::copy(std_matrix, vcl_coordinate_matrix);
382 
383   // --------------------------------------------------------------------------
384   std::cout << "Testing products: STL" << std::endl;
385   result = viennacl::linalg::prod(std_matrix, rhs);
386 
387   std::cout << "Testing products: compressed_matrix" << std::endl;
388   vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
389 
390   if ( std::fabs(diff(result, vcl_result)) > epsilon )
391   {
392     std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
393     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
394     return EXIT_FAILURE;
395   }
396 
397   std::cout << "Testing products: compressed_matrix, strided vectors" << std::endl;
398   retval = strided_matrix_vector_product_test<NumericT, viennacl::compressed_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
399   if (retval != EXIT_SUCCESS)
400     return retval;
401 
402   result = rhs;
403   result = viennacl::linalg::prod(std_matrix, rhs);
404   for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
405   vcl_result = vcl_rhs;
406   vcl_result += viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
407 
408   if ( std::fabs(diff(result, vcl_result)) > epsilon )
409   {
410     std::cout << "# Error at operation: matrix-vector product with compressed_matrix (+=)" << std::endl;
411     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
412     return EXIT_FAILURE;
413   }
414 
415   result = rhs;
416   result = viennacl::linalg::prod(std_matrix, rhs);
417   for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
418   vcl_result = vcl_rhs;
419   vcl_result -= viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
420 
421   if ( std::fabs(diff(result, vcl_result)) > epsilon )
422   {
423     std::cout << "# Error at operation: matrix-vector product with compressed_matrix (-=)" << std::endl;
424     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
425     return EXIT_FAILURE;
426   }
427 
428   //
429   // Triangular solvers for A \ b:
430   //
431 
432   std::cout << "Testing unit upper triangular solve: compressed_matrix" << std::endl;
433   result = rhs;
434   viennacl::copy(result, vcl_result);
435   //boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::unit_upper_tag());
436   for (std::size_t i2=0; i2<std_matrix.size(); ++i2)
437   {
438     std::size_t row = std_matrix.size() - i2 - 1;
439     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[row].begin(); it != std_matrix[row].end(); ++it)
440     {
441       if (it->first > static_cast<unsigned int>(row))
442         result[row] -= it->second * result[it->first];
443       else
444         continue;
445     }
446   }
447 
448   viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::unit_upper_tag());
449 
450   if ( std::fabs(diff(result, vcl_result)) > epsilon )
451   {
452     std::cout << "# Error at operation: unit upper triangular solve with compressed_matrix" << std::endl;
453     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
454     return EXIT_FAILURE;
455   }
456 
457   ////////////////////////////
458 
459   std::cout << "Testing upper triangular solve: compressed_matrix" << std::endl;
460   result = rhs;
461   viennacl::copy(result, vcl_result);
462   //boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::upper_tag());
463   for (std::size_t i2=0; i2<std_matrix.size(); ++i2)
464   {
465     std::size_t row = std_matrix.size() - i2 - 1;
466     NumericT diag = 0;
467     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[row].begin(); it != std_matrix[row].end(); ++it)
468     {
469       if (it->first > static_cast<unsigned int>(row))
470         result[row] -= it->second * result[it->first];
471       else if (it->first == static_cast<unsigned int>(row))
472         diag = it->second;
473       else
474         continue;
475     }
476     result[row] /= diag;
477   }
478 
479   viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::upper_tag());
480 
481   if ( std::fabs(diff(result, vcl_result)) > epsilon )
482   {
483     std::cout << "# Error at operation: upper triangular solve with compressed_matrix" << std::endl;
484     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
485     return EXIT_FAILURE;
486   }
487 
488   ////////////////////////////
489 
490   std::cout << "Testing unit lower triangular solve: compressed_matrix" << std::endl;
491   result = rhs;
492   viennacl::copy(result, vcl_result);
493   //boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::unit_lower_tag());
494   for (std::size_t i=1; i<std_matrix.size(); ++i)
495     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
496     {
497       if (it->first < static_cast<unsigned int>(i))
498         result[i] -= it->second * result[it->first];
499       else
500         continue;
501     }
502   viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::unit_lower_tag());
503 
504 
505   if ( std::fabs(diff(result, vcl_result)) > epsilon )
506   {
507     std::cout << "# Error at operation: unit lower triangular solve with compressed_matrix" << std::endl;
508     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
509     return EXIT_FAILURE;
510   }
511 
512 
513   std::cout << "Testing lower triangular solve: compressed_matrix" << std::endl;
514   result = rhs;
515   viennacl::copy(result, vcl_result);
516   //boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::lower_tag());
517   for (std::size_t i=0; i<std_matrix.size(); ++i)
518   {
519     NumericT diag = 0;
520     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
521     {
522       if (it->first < static_cast<unsigned int>(i))
523         result[i] -= it->second * result[it->first];
524       else if (it->first == static_cast<unsigned int>(i))
525         diag = it->second;
526       else
527         continue;
528     }
529     result[i] /= diag;
530   }
531   viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::lower_tag());
532 
533 
534   if ( std::fabs(diff(result, vcl_result)) > epsilon )
535   {
536     std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
537     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
538     return EXIT_FAILURE;
539   }
540 
541 
542 
543   //
544   // Triangular solvers for A^T \ b
545   //
546   std::vector<std::map<unsigned int, NumericT> > std_matrix_trans(std_matrix.size());
547 
548   // compute transpose:
549   for (std::size_t i=0; i<std_matrix.size(); ++i)
550     for (typename std::map<unsigned int, NumericT>::const_iterator it  = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
551       std_matrix_trans[i][it->first] = it->second;
552 
553   std::cout << "Testing transposed unit upper triangular solve: compressed_matrix" << std::endl;
554   result = rhs;
555   viennacl::copy(result, vcl_result);
556   //boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::unit_upper_tag());
557   for (std::size_t i2=0; i2<std_matrix_trans.size(); ++i2)
558   {
559     std::size_t row = std_matrix_trans.size() - i2 - 1;
560     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[row].begin(); it != std_matrix_trans[row].end(); ++it)
561     {
562       if (it->first > static_cast<unsigned int>(row))
563         result[row] -= it->second * result[it->first];
564       else
565         continue;
566     }
567   }
568   viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::unit_upper_tag());
569 
570   if ( std::fabs(diff(result, vcl_result)) > epsilon )
571   {
572     std::cout << "# Error at operation: unit upper triangular solve with compressed_matrix" << std::endl;
573     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
574     return EXIT_FAILURE;
575   }
576 
577   /////////////////////////
578 
579   std::cout << "Testing transposed upper triangular solve: compressed_matrix" << std::endl;
580   result = rhs;
581   viennacl::copy(result, vcl_result);
582   //boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::upper_tag());
583   for (std::size_t i2=0; i2<std_matrix_trans.size(); ++i2)
584   {
585     std::size_t row = std_matrix_trans.size() - i2 - 1;
586     NumericT diag = 0;
587     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[row].begin(); it != std_matrix_trans[row].end(); ++it)
588     {
589       if (it->first > static_cast<unsigned int>(row))
590         result[row] -= it->second * result[it->first];
591       else if (it->first == static_cast<unsigned int>(row))
592         diag = it->second;
593       else
594         continue;
595     }
596     result[row] /= diag;
597   }
598   viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::upper_tag());
599 
600   if ( std::fabs(diff(result, vcl_result)) > epsilon )
601   {
602     std::cout << "# Error at operation: upper triangular solve with compressed_matrix" << std::endl;
603     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
604     return EXIT_FAILURE;
605   }
606 
607   /////////////////////////
608 
609   std::cout << "Testing transposed unit lower triangular solve: compressed_matrix" << std::endl;
610   result = rhs;
611   viennacl::copy(result, vcl_result);
612   //boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::unit_lower_tag());
613   for (std::size_t i=1; i<std_matrix_trans.size(); ++i)
614     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[i].begin(); it != std_matrix_trans[i].end(); ++it)
615     {
616       if (it->first < static_cast<unsigned int>(i))
617         result[i] -= it->second * result[it->first];
618       else
619         continue;
620     }
621   viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::unit_lower_tag());
622 
623   if ( std::fabs(diff(result, vcl_result)) > epsilon )
624   {
625     std::cout << "# Error at operation: unit lower triangular solve with compressed_matrix" << std::endl;
626     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
627     return EXIT_FAILURE;
628   }
629 
630   /////////////////////////
631 
632   std::cout << "Testing transposed lower triangular solve: compressed_matrix" << std::endl;
633   result = rhs;
634   viennacl::copy(result, vcl_result);
635   //boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::lower_tag());
636   for (std::size_t i=0; i<std_matrix_trans.size(); ++i)
637   {
638     NumericT diag = 0;
639     for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[i].begin(); it != std_matrix_trans[i].end(); ++it)
640     {
641       if (it->first < static_cast<unsigned int>(i))
642         result[i] -= it->second * result[it->first];
643       else if (it->first == static_cast<unsigned int>(i))
644         diag = it->second;
645       else
646         continue;
647     }
648     result[i] /= diag;
649   }
650   viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::lower_tag());
651 
652   if ( std::fabs(diff(result, vcl_result)) > epsilon )
653   {
654     std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
655     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
656     return EXIT_FAILURE;
657   }
658 
659 
660   //
661   /////////////////////////
662   //
663 
664 
665   std::cout << "Testing products: compressed_compressed_matrix" << std::endl;
666   result     = viennacl::linalg::prod(std_cc_matrix, rhs);
667   vcl_result = viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
668 
669   if ( std::fabs(diff(result, vcl_result)) > epsilon )
670   {
671     std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (=)" << std::endl;
672     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
673     return EXIT_FAILURE;
674   }
675 
676   {
677     std::vector<std::map<unsigned int, NumericT> > temp(vcl_compressed_compressed_matrix.size1());
678     viennacl::copy(vcl_compressed_compressed_matrix, temp);
679 
680     // check that entries are correct by computing the product again:
681     result     = viennacl::linalg::prod(temp, rhs);
682 
683     if ( std::fabs(diff(result, vcl_result)) > epsilon )
684     {
685       std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (after copy back)" << std::endl;
686       std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
687       return EXIT_FAILURE;
688     }
689 
690   }
691 
692   result = rhs;
693   result = viennacl::linalg::prod(std_cc_matrix, rhs);
694   for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
695   vcl_result = vcl_rhs;
696   vcl_result += viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
697 
698   if ( std::fabs(diff(result, vcl_result)) > epsilon )
699   {
700     std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (+=)" << std::endl;
701     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
702     return EXIT_FAILURE;
703   }
704 
705   result = rhs;
706   result = viennacl::linalg::prod(std_cc_matrix, rhs);
707   for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
708   vcl_result = vcl_rhs;
709   vcl_result -= viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
710 
711   if ( std::fabs(diff(result, vcl_result)) > epsilon )
712   {
713     std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (-=)" << std::endl;
714     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
715     return EXIT_FAILURE;
716   }
717 
718 
719   //
720   /////////////////////////
721   //
722 
723 
724   std::cout << "Testing products: coordinate_matrix" << std::endl;
725   result     = viennacl::linalg::prod(std_matrix, rhs);
726   vcl_result = viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
727 
728   if ( std::fabs(diff(result, vcl_result)) > epsilon )
729   {
730     std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
731     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
732     return EXIT_FAILURE;
733   }
734 
735   std::cout << "Testing products: coordinate_matrix, strided vectors" << std::endl;
736   //std::cout << " --> SKIPPING <--" << std::endl;
737   retval = strided_matrix_vector_product_test<NumericT, viennacl::coordinate_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
738   if (retval != EXIT_SUCCESS)
739     return retval;
740 
741   result = rhs;
742   result = viennacl::linalg::prod(std_matrix, rhs);
743   for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
744   vcl_result = vcl_rhs;
745   vcl_result += viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
746 
747   if ( std::fabs(diff(result, vcl_result)) > epsilon )
748   {
749     std::cout << "# Error at operation: matrix-vector product with coordinate_matrix (+=)" << std::endl;
750     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
751     return EXIT_FAILURE;
752   }
753 
754   result = rhs;
755   result = viennacl::linalg::prod(std_matrix, rhs);
756   for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
757   vcl_result = vcl_rhs;
758   vcl_result -= viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
759 
760   if ( std::fabs(diff(result, vcl_result)) > epsilon )
761   {
762     std::cout << "# Error at operation: matrix-vector product with coordinate_matrix (-=)" << std::endl;
763     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
764     return EXIT_FAILURE;
765   }
766 
767   //
768   /////////////////////////
769   //
770 
771 
772   //std::cout << "Copying ell_matrix" << std::endl;
773   viennacl::copy(std_matrix, vcl_ell_matrix);
774   std_matrix.clear();
775   viennacl::copy(vcl_ell_matrix, std_matrix);// just to check that it works
776 
777 
778   std::cout << "Testing products: ell_matrix" << std::endl;
779   result     = viennacl::linalg::prod(std_matrix, rhs);
780   vcl_result.clear();
781   vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
782   //viennacl::linalg::prod_impl(vcl_ell_matrix, vcl_rhs, vcl_result);
783   //std::cout << vcl_result << "\n";
784   //std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
785   //std::cout << "First entry of result vector: " << vcl_result[0] << std::endl;
786 
787   if ( std::fabs(diff(result, vcl_result)) > epsilon )
788   {
789     std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
790     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
791     return EXIT_FAILURE;
792   }
793 
794   std::cout << "Testing products: ell_matrix, strided vectors" << std::endl;
795   retval = strided_matrix_vector_product_test<NumericT, viennacl::ell_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
796   if (retval != EXIT_SUCCESS)
797     return retval;
798 
799   result = rhs;
800   result = viennacl::linalg::prod(std_matrix, rhs);
801   for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
802   vcl_result = vcl_rhs;
803   vcl_result += viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
804 
805   if ( std::fabs(diff(result, vcl_result)) > epsilon )
806   {
807     std::cout << "# Error at operation: matrix-vector product with ell_matrix (+=)" << std::endl;
808     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
809     return EXIT_FAILURE;
810   }
811 
812   result = rhs;
813   result = viennacl::linalg::prod(std_matrix, rhs);
814   for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
815   vcl_result = vcl_rhs;
816   vcl_result -= viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
817 
818   if ( std::fabs(diff(result, vcl_result)) > epsilon )
819   {
820     std::cout << "# Error at operation: matrix-vector product with ell_matrix (-=)" << std::endl;
821     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
822     return EXIT_FAILURE;
823   }
824 
825   //
826   /////////////////////////
827   //
828 
829 
830   //std::cout << "Copying sliced_ell_matrix" << std::endl;
831   viennacl::copy(std_matrix, vcl_sliced_ell_matrix);
832 
833   std::cout << "Testing products: sliced_ell_matrix" << std::endl;
834   result     = viennacl::linalg::prod(std_matrix, rhs);
835   vcl_result.clear();
836   vcl_result = viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
837 
838   if ( std::fabs(diff(result, vcl_result)) > epsilon )
839   {
840     std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix" << std::endl;
841     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
842     return EXIT_FAILURE;
843   }
844 
845   std::cout << "Testing products: sliced_ell_matrix, strided vectors" << std::endl;
846   retval = strided_matrix_vector_product_test<NumericT, viennacl::sliced_ell_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
847   if (retval != EXIT_SUCCESS)
848     return retval;
849 
850   result = rhs;
851   result = viennacl::linalg::prod(std_matrix, rhs);
852   for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
853   vcl_result = vcl_rhs;
854   vcl_result += viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
855 
856   if ( std::fabs(diff(result, vcl_result)) > epsilon )
857   {
858     std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix (+=)" << std::endl;
859     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
860     return EXIT_FAILURE;
861   }
862 
863   result = rhs;
864   result = viennacl::linalg::prod(std_matrix, rhs);
865   for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
866   vcl_result = vcl_rhs;
867   vcl_result -= viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
868 
869   if ( std::fabs(diff(result, vcl_result)) > epsilon )
870   {
871     std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix (-=)" << std::endl;
872     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
873     return EXIT_FAILURE;
874   }
875 
876   //
877   /////////////////////////
878   //
879 
880 
881   //std::cout << "Copying hyb_matrix" << std::endl;
882   viennacl::copy(std_matrix, vcl_hyb_matrix);
883   std_matrix.clear();
884   viennacl::copy(vcl_hyb_matrix, std_matrix);// just to check that it works
885   viennacl::copy(std_matrix, vcl_hyb_matrix);
886 
887   std::cout << "Testing products: hyb_matrix" << std::endl;
888   result     = viennacl::linalg::prod(std_matrix, rhs);
889   vcl_result.clear();
890   vcl_result = viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
891   //viennacl::linalg::prod_impl(vcl_hyb_matrix, vcl_rhs, vcl_result);
892   //std::cout << vcl_result << "\n";
893   //std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
894   //std::cout << "First entry of result vector: " << vcl_result[0] << std::endl;
895 
896   if ( std::fabs(diff(result, vcl_result)) > epsilon )
897   {
898     std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
899     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
900     return EXIT_FAILURE;
901   }
902 
903   std::cout << "Testing products: hyb_matrix, strided vectors" << std::endl;
904   retval = strided_matrix_vector_product_test<NumericT, viennacl::hyb_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
905   if (retval != EXIT_SUCCESS)
906     return retval;
907 
908   result = rhs;
909   result = viennacl::linalg::prod(std_matrix, rhs);
910   for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
911   vcl_result = vcl_rhs;
912   vcl_result += viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
913 
914   if ( std::fabs(diff(result, vcl_result)) > epsilon )
915   {
916     std::cout << "# Error at operation: matrix-vector product with hyb_matrix (+=)" << std::endl;
917     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
918     return EXIT_FAILURE;
919   }
920 
921   result = rhs;
922   result = viennacl::linalg::prod(std_matrix, rhs);
923   for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
924   vcl_result = vcl_rhs;
925   vcl_result -= viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
926 
927   if ( std::fabs(diff(result, vcl_result)) > epsilon )
928   {
929     std::cout << "# Error at operation: matrix-vector product with hyb_matrix (-=)" << std::endl;
930     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
931     return EXIT_FAILURE;
932   }
933 
934 
935   // --------------------------------------------------------------------------
936   // --------------------------------------------------------------------------
937   NumericT alpha = static_cast<NumericT>(2.786);
938   NumericT beta = static_cast<NumericT>(1.432);
939   copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
940   copy(result.begin(), result.end(), vcl_result.begin());
941   copy(result.begin(), result.end(), vcl_result2.begin());
942 
943   std::cout << "Testing scaled additions of products and vectors" << std::endl;
944   std::vector<NumericT> result2(result);
945   result2 = viennacl::linalg::prod(std_matrix, rhs);
946   for (std::size_t i=0; i<result.size(); ++i)
947     result[i] = alpha * result2[i] + beta * result[i];
948   vcl_result2 = alpha * viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs) + beta * vcl_result;
949 
950   if ( std::fabs(diff(result, vcl_result2)) > epsilon )
951   {
952     std::cout << "# Error at operation: matrix-vector product (compressed_matrix) with scaled additions" << std::endl;
953     std::cout << "  diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
954     return EXIT_FAILURE;
955   }
956 
957 
958   vcl_result2.clear();
959   vcl_result2 = alpha * viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs) + beta * vcl_result;
960 
961   if ( std::fabs(diff(result, vcl_result2)) > epsilon )
962   {
963     std::cout << "# Error at operation: matrix-vector product (coordinate_matrix) with scaled additions" << std::endl;
964     std::cout << "  diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
965     return EXIT_FAILURE;
966   }
967 
968   vcl_result2.clear();
969   vcl_result2 = alpha * viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs) + beta * vcl_result;
970 
971   if ( std::fabs(diff(result, vcl_result2)) > epsilon )
972   {
973     std::cout << "# Error at operation: matrix-vector product (ell_matrix) with scaled additions" << std::endl;
974     std::cout << "  diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
975     return EXIT_FAILURE;
976   }
977 
978   vcl_result2.clear();
979   vcl_result2 = alpha * viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs) + beta * vcl_result;
980 
981   if ( std::fabs(diff(result, vcl_result2)) > epsilon )
982   {
983     std::cout << "# Error at operation: matrix-vector product (hyb_matrix) with scaled additions" << std::endl;
984     std::cout << "  diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
985     return EXIT_FAILURE;
986   }
987 
988   ////////////// Test of .clear() ////////////////
989   std_matrix.clear();
990 
991   std::cout << "Testing products after clear(): compressed_matrix" << std::endl;
992   vcl_compressed_matrix.clear();
993   result     = std::vector<NumericT>(result.size(), NumericT(1));
994   result     = viennacl::linalg::prod(std_matrix, rhs);
995   vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
996 
997   if ( std::fabs(diff(result, vcl_result)) > epsilon )
998   {
999     std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
1000     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1001     return EXIT_FAILURE;
1002   }
1003 
1004   std::cout << "Testing products after clear(): compressed_compressed_matrix" << std::endl;
1005   vcl_compressed_compressed_matrix.clear();
1006   result     = std::vector<NumericT>(result.size(), NumericT(1));
1007   result     = viennacl::linalg::prod(std_matrix, rhs);
1008   vcl_result = viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
1009 
1010   if ( std::fabs(diff(result, vcl_result)) > epsilon )
1011   {
1012     std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix" << std::endl;
1013     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1014     return EXIT_FAILURE;
1015   }
1016 
1017   std::cout << "Testing products after clear(): coordinate_matrix" << std::endl;
1018   vcl_coordinate_matrix.clear();
1019   result     = std::vector<NumericT>(result.size(), NumericT(1));
1020   result     = viennacl::linalg::prod(std_matrix, rhs);
1021   vcl_result = viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
1022 
1023   if ( std::fabs(diff(result, vcl_result)) > epsilon )
1024   {
1025     std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
1026     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1027     return EXIT_FAILURE;
1028   }
1029 
1030   std::cout << "Testing products after clear(): ell_matrix" << std::endl;
1031   vcl_ell_matrix.clear();
1032   result     = std::vector<NumericT>(result.size(), NumericT(1));
1033   result     = viennacl::linalg::prod(std_matrix, rhs);
1034   vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
1035 
1036   if ( std::fabs(diff(result, vcl_result)) > epsilon )
1037   {
1038     std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
1039     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1040     return EXIT_FAILURE;
1041   }
1042 
1043   std::cout << "Testing products after clear(): hyb_matrix" << std::endl;
1044   vcl_hyb_matrix.clear();
1045   result     = std::vector<NumericT>(result.size(), NumericT(1));
1046   result     = viennacl::linalg::prod(std_matrix, rhs);
1047   vcl_result = viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
1048 
1049   if ( std::fabs(diff(result, vcl_result)) > epsilon )
1050   {
1051     std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
1052     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1053     return EXIT_FAILURE;
1054   }
1055 
1056   std::cout << "Testing products after clear(): sliced_ell_matrix" << std::endl;
1057   vcl_sliced_ell_matrix.clear();
1058   result     = std::vector<NumericT>(result.size(), NumericT(1));
1059   result     = viennacl::linalg::prod(std_matrix, rhs);
1060   vcl_result = viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
1061 
1062   if ( std::fabs(diff(result, vcl_result)) > epsilon )
1063   {
1064     std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix" << std::endl;
1065     std::cout << "  diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1066     return EXIT_FAILURE;
1067   }
1068 
1069 
1070   // --------------------------------------------------------------------------
1071   return retval;
1072 }
1073 //
1074 // -------------------------------------------------------------
1075 //
main()1076 int main()
1077 {
1078   std::cout << std::endl;
1079   std::cout << "----------------------------------------------" << std::endl;
1080   std::cout << "----------------------------------------------" << std::endl;
1081   std::cout << "## Test :: Sparse Matrices" << std::endl;
1082   std::cout << "----------------------------------------------" << std::endl;
1083   std::cout << "----------------------------------------------" << std::endl;
1084   std::cout << std::endl;
1085 
1086   int retval = EXIT_SUCCESS;
1087 
1088   std::cout << std::endl;
1089   std::cout << "----------------------------------------------" << std::endl;
1090   std::cout << std::endl;
1091   {
1092     typedef float NumericT;
1093     NumericT epsilon = static_cast<NumericT>(1E-4);
1094     std::cout << "# Testing setup:" << std::endl;
1095     std::cout << "  eps:     " << epsilon << std::endl;
1096     std::cout << "  numeric: float" << std::endl;
1097     retval = test<NumericT>(epsilon);
1098     if ( retval == EXIT_SUCCESS )
1099         std::cout << "# Test passed" << std::endl;
1100     else
1101         return retval;
1102   }
1103   std::cout << std::endl;
1104   std::cout << "----------------------------------------------" << std::endl;
1105   std::cout << std::endl;
1106 
1107 #ifdef VIENNACL_WITH_OPENCL
1108   if ( viennacl::ocl::current_device().double_support() )
1109 #endif
1110   {
1111     {
1112       typedef double NumericT;
1113       NumericT epsilon = 1.0E-12;
1114       std::cout << "# Testing setup:" << std::endl;
1115       std::cout << "  eps:     " << epsilon << std::endl;
1116       std::cout << "  numeric: double" << std::endl;
1117       retval = test<NumericT>(epsilon);
1118       if ( retval == EXIT_SUCCESS )
1119         std::cout << "# Test passed" << std::endl;
1120       else
1121         return retval;
1122     }
1123     std::cout << std::endl;
1124     std::cout << "----------------------------------------------" << std::endl;
1125     std::cout << std::endl;
1126   }
1127 #ifdef VIENNACL_WITH_OPENCL
1128   else
1129     std::cout << "No double precision support, skipping test..." << std::endl;
1130 #endif
1131 
1132 
1133   std::cout << std::endl;
1134   std::cout << "------- Test completed --------" << std::endl;
1135   std::cout << std::endl;
1136 
1137   return retval;
1138 }
1139