35 #include <boost/numeric/ublas/io.hpp>
36 #include <boost/numeric/ublas/triangular.hpp>
37 #include <boost/numeric/ublas/matrix_sparse.hpp>
38 #include <boost/numeric/ublas/matrix.hpp>
39 #include <boost/numeric/ublas/matrix_proxy.hpp>
40 #include <boost/numeric/ublas/lu.hpp>
41 #include <boost/numeric/ublas/io.hpp>
42 #include <boost/numeric/ublas/operation_sparse.hpp>
48 #define VIENNACL_WITH_UBLAS 1
73 template<
typename ScalarType>
77 return (s1 - s2) /
std::max(fabs(s1), std::fabs(s2));
81 template<
typename ScalarType>
84 ublas::vector<ScalarType> v2_cpu(v2.
size());
88 for (
unsigned int i=0;i<v1.size(); ++i)
90 if (
std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
95 v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) /
std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
100 if (v2_cpu[i] > 0.0001)
103 std::cout <<
"Error at entry " << i <<
": " << v1[i] <<
" vs. " << v2_cpu[i] << std::endl;
113 template<
typename ScalarType,
typename VCL_MATRIX>
114 ScalarType diff(ublas::compressed_matrix<ScalarType> & cpu_matrix, VCL_MATRIX & gpu_matrix)
116 typedef ublas::compressed_matrix<ScalarType> CPU_MATRIX;
126 for (
typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
127 row_it != cpu_matrix.end1();
131 for (
typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
132 col_it != row_it.end();
138 if (
std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
139 std::fabs(from_gpu(col_it.index1(), col_it.index2())) ) > 0 )
140 current_error = std::fabs(cpu_matrix(col_it.index1(), col_it.index2()) - from_gpu(col_it.index1(), col_it.index2()))
141 /
std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
142 std::fabs(from_gpu(col_it.index1(), col_it.index2())) );
143 if (current_error > error)
144 error = current_error;
150 for (
typename CPU_MATRIX::const_iterator1 row_it = from_gpu.begin1();
151 row_it != from_gpu.end1();
155 for (
typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
156 col_it != row_it.end();
162 if (
std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
163 std::fabs(from_gpu(col_it.index1(), col_it.index2())) ) > 0 )
164 current_error = std::fabs(cpu_matrix(col_it.index1(), col_it.index2()) - from_gpu(col_it.index1(), col_it.index2()))
165 /
std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
166 std::fabs(from_gpu(col_it.index1(), col_it.index2())) );
167 if (current_error > error)
168 error = current_error;
178 template<
typename NumericT,
typename Epsilon >
179 int test(Epsilon
const& epsilon)
181 int retval = EXIT_SUCCESS;
189 ublas::vector<NumericT> rhs;
190 ublas::vector<NumericT> result;
191 ublas::compressed_matrix<NumericT> ublas_matrix;
195 std::cout <<
"Error reading Matrix file" << std::endl;
199 std::cout <<
"done reading matrix" << std::endl;
202 rhs.resize(ublas_matrix.size2());
203 for (std::size_t i=0; i<rhs.size(); ++i)
206 rhs[i] =
NumericT(1) + randomNumber();
225 std::cout <<
"Testing products: compressed_matrix" << std::endl;
233 if ( std::fabs(
diff(result, vcl_result)) > epsilon )
235 std::cout <<
"# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
236 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result)) << std::endl;
237 retval = EXIT_FAILURE;
240 std::cout <<
"Testing products: coordinate_matrix" << std::endl;
249 if ( std::fabs(
diff(result, vcl_result)) > epsilon )
251 std::cout <<
"# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
252 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result)) << std::endl;
253 retval = EXIT_FAILURE;
262 if ( std::fabs(
diff(result, vcl_result2)) > epsilon )
264 std::cout <<
"# Error at operation: matrix-vector product (coordinate_matrix) with scaled additions" << std::endl;
265 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result2)) << std::endl;
266 retval = EXIT_FAILURE;
271 ublas_matrix.clear();
275 std::cout <<
"Testing products: ell_matrix" << std::endl;
285 if ( std::fabs(
diff(result, vcl_result)) > epsilon )
287 std::cout <<
"# Error at operation: matrix-vector product with ell_matrix" << std::endl;
288 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result)) << std::endl;
289 retval = EXIT_FAILURE;
294 ublas_matrix.clear();
298 std::cout <<
"Testing products: hyb_matrix" << std::endl;
307 if ( std::fabs(
diff(result, vcl_result)) > epsilon )
309 std::cout <<
"# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
310 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result)) << std::endl;
311 retval = EXIT_FAILURE;
317 copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
318 copy(result.begin(), result.end(), vcl_result.begin());
319 copy(result.begin(), result.end(), vcl_result2.begin());
320 copy(ublas_matrix, vcl_compressed_matrix);
321 copy(ublas_matrix, vcl_coordinate_matrix);
322 copy(ublas_matrix, vcl_ell_matrix);
323 copy(ublas_matrix, vcl_hyb_matrix);
325 std::cout <<
"Testing scaled additions of products and vectors: compressed_matrix" << std::endl;
334 if ( std::fabs(
diff(result, vcl_result2)) > epsilon )
336 std::cout <<
"# Error at operation: matrix-vector product (compressed_matrix) with scaled additions" << std::endl;
337 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result2)) << std::endl;
338 retval = EXIT_FAILURE;
342 std::cout <<
"Testing scaled additions of products and vectors: coordinate_matrix" << std::endl;
343 copy(result.begin(), result.end(), vcl_result.begin());
352 if ( std::fabs(
diff(result, vcl_result2)) > epsilon )
354 std::cout <<
"# Error at operation: matrix-vector product (coordinate_matrix) with scaled additions" << std::endl;
355 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result2)) << std::endl;
356 retval = EXIT_FAILURE;
359 std::cout <<
"Testing scaled additions of products and vectors: ell_matrix" << std::endl;
360 copy(result.begin(), result.end(), vcl_result.begin());
369 if ( std::fabs(
diff(result, vcl_result2)) > epsilon )
371 std::cout <<
"# Error at operation: matrix-vector product (ell_matrix) with scaled additions" << std::endl;
372 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result2)) << std::endl;
373 retval = EXIT_FAILURE;
376 std::cout <<
"Testing scaled additions of products and vectors: hyb_matrix" << std::endl;
377 copy(result.begin(), result.end(), vcl_result.begin());
386 if ( std::fabs(
diff(result, vcl_result2)) > epsilon )
388 std::cout <<
"# Error at operation: matrix-vector product (hyb_matrix) with scaled additions" << std::endl;
389 std::cout <<
" diff: " << std::fabs(
diff(result, vcl_result2)) << std::endl;
390 retval = EXIT_FAILURE;
402 std::cout << std::endl;
403 std::cout <<
"----------------------------------------------" << std::endl;
404 std::cout <<
"----------------------------------------------" << std::endl;
405 std::cout <<
"## Test :: Sparse Matrices" << std::endl;
406 std::cout <<
"----------------------------------------------" << std::endl;
407 std::cout <<
"----------------------------------------------" << std::endl;
408 std::cout << std::endl;
410 int retval = EXIT_SUCCESS;
412 std::cout << std::endl;
413 std::cout <<
"----------------------------------------------" << std::endl;
414 std::cout << std::endl;
417 NumericT epsilon =
static_cast<NumericT
>(1E-4);
418 std::cout <<
"# Testing setup:" << std::endl;
419 std::cout <<
" eps: " << epsilon << std::endl;
420 std::cout <<
" numeric: float" << std::endl;
421 retval = test<NumericT>(epsilon);
422 if ( retval == EXIT_SUCCESS )
423 std::cout <<
"# Test passed" << std::endl;
427 std::cout << std::endl;
428 std::cout <<
"----------------------------------------------" << std::endl;
429 std::cout << std::endl;
431 #ifdef VIENNACL_WITH_OPENCL
437 NumericT epsilon = 1.0E-13;
438 std::cout <<
"# Testing setup:" << std::endl;
439 std::cout <<
" eps: " << epsilon << std::endl;
440 std::cout <<
" numeric: double" << std::endl;
441 retval = test<NumericT>(epsilon);
442 if ( retval == EXIT_SUCCESS )
443 std::cout <<
"# Test passed" << std::endl;
447 std::cout << std::endl;
448 std::cout <<
"----------------------------------------------" << std::endl;
449 std::cout << std::endl;
451 #ifdef VIENNACL_WITH_OPENCL
453 std::cout <<
"No double precision support, skipping test..." << std::endl;
457 std::cout << std::endl;
458 std::cout <<
"------- Test completed --------" << std::endl;
459 std::cout << std::endl;
Sparse matrix class using a hybrid format composed of the ELL and CSR format for storing the nonzeros...
A reader and writer for the matrix market format is implemented here.
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
int test(Epsilon const &epsilon)
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Some helper routines for reading/writing/printing scheduler expressions.
A tag class representing assignment.
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
void execute(statement const &s)
viennacl::scalar< int > s2
viennacl::scalar< float > s1
T max(const T &lhs, const T &rhs)
Maximum.
viennacl::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Implementation of the coordinate_matrix class.
viennacl::vector< float > v1
Implementation of the hyb_matrix class.
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Sparse matrix class using the ELLPACK format for storing the nonzeros.
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
Implementations of incomplete factorization preconditioners. Convenience header file.
Implementation of the compressed_matrix class.
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Implementation of the ell_matrix class.
Proxy classes for vectors.
viennacl::vector< int > v2
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
T norm_inf(std::vector< T, A > const &v1)
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) ...
A small collection of sequential random number generators.
size_type size() const
Returns the length of the vector (cf. std::vector)
Provides the datastructures for dealing with a single statement such as 'x = y + z;'.
The main class for representing a statement such as x = inner_prod(y,z); at runtime.
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
long read_matrix_market_file(MatrixT &mat, const char *file, long index_base=1)
Reads a sparse matrix from a file (MatrixMarket format)
iterator end()
Returns an iterator pointing to the end of the vector (STL like)
Common routines used within ILU-type preconditioners.
Implementation of the ViennaCL scalar class.
A sparse square matrix, where entries are stored as triplets (i,j, val), where i and j are the row an...