24 #define VIENNACL_WITH_UBLAS
48 #include "boost/numeric/ublas/vector.hpp"
49 #include "boost/numeric/ublas/matrix.hpp"
50 #include "boost/numeric/ublas/matrix_proxy.hpp"
51 #include "boost/numeric/ublas/vector_proxy.hpp"
52 #include "boost/numeric/ublas/io.hpp"
56 using namespace boost::numeric;
58 template<
typename MatrixType,
typename VCLMatrixType>
61 typedef typename MatrixType::value_type value_type;
63 boost::numeric::ublas::matrix<value_type> vcl_A_cpu(vcl_A.size1(), vcl_A.size2());
67 for (std::size_t i=0; i<ublas_A.size1(); ++i)
69 for (std::size_t j=0; j<ublas_A.size2(); ++j)
71 if (std::fabs(ublas_A(i,j) - vcl_A_cpu(i,j)) > 0)
73 if ( (std::fabs(ublas_A(i,j) - vcl_A_cpu(i,j)) /
std::max(std::fabs(ublas_A(i,j)), std::fabs(vcl_A_cpu(i,j))) > epsilon) || std::fabs(vcl_A_cpu(i,j) - vcl_A_cpu(i,j)) > 0 )
75 std::cout <<
"Error at index (" << i <<
", " << j <<
"): " << ublas_A(i,j) <<
" vs " << vcl_A_cpu(i,j) << std::endl;
76 std::cout << std::endl <<
"TEST failed!" << std::endl;
83 std::cout <<
"PASSED!" << std::endl;
90 template<
typename UBLASMatrixType,
91 typename ViennaCLMatrixType1,
typename ViennaCLMatrixType2,
typename ViennaCLMatrixType3>
93 UBLASMatrixType & ublas_A, UBLASMatrixType & ublas_B, UBLASMatrixType & ublas_C,
94 ViennaCLMatrixType1 & vcl_A, ViennaCLMatrixType2 & vcl_B, ViennaCLMatrixType3 vcl_C)
99 cpu_value_type alpha = cpu_value_type(3.1415);
102 cpu_value_type beta = cpu_value_type(2.7182);
109 std::cout <<
"Checking for zero_matrix initializer..." << std::endl;
110 ublas_A = ublas::zero_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2());
115 std::cout <<
"Checking for scalar_matrix initializer..." << std::endl;
116 ublas_A = ublas::scalar_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2(), alpha);
121 ublas_A = ublas::scalar_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2(), gpu_beta);
133 std::cout << std::endl;
141 std::cout <<
"Testing matrix assignment... ";
157 std::cout <<
"Testing upper left copy to GPU... ";
164 std::cout <<
"Testing lower right copy to GPU... ";
174 std::cout <<
"Testing upper left copy to A... ";
178 std::cout <<
"Testing lower right copy to C... ";
189 std::cout <<
"Assignment: ";
200 std::cout <<
"Inplace add: ";
210 std::cout <<
"Inplace sub: ";
221 std::cout <<
"Add: ";
223 ublas_C = ublas_A + ublas_B;
231 std::cout <<
"Sub: ";
233 ublas_C = ublas_A - ublas_B;
241 std::cout <<
"Composite assignments: ";
243 ublas_C += alpha * ublas_A - beta * ublas_B + ublas_A / beta - ublas_B / alpha;
252 std::cout <<
"--- Testing elementwise operations (binary) ---" << std::endl;
253 std::cout <<
"x = element_prod(x, y)... ";
263 std::cout <<
"x = element_prod(x + y, y)... ";
273 std::cout <<
"x = element_prod(x, x + y)... ";
283 std::cout <<
"x = element_prod(x - y, y + x)... ";
285 ublas_C =
element_prod(ublas_A - ublas_B, ublas_B + ublas_A);
295 std::cout <<
"x = element_div(x, y)... ";
305 std::cout <<
"x = element_div(x + y, y)... ";
315 std::cout <<
"x = element_div(x, x + y)... ";
325 std::cout <<
"x = element_div(x - y, y + x)... ";
327 ublas_C =
element_div(ublas_A - ublas_B, ublas_B + ublas_A);
336 std::cout <<
"--- Testing elementwise operations (unary) ---" << std::endl;
337 #define GENERATE_UNARY_OP_TEST(OPNAME) \
338 ublas_A = ublas::scalar_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2(), cpu_value_type(0.21)); \
339 ublas_B = cpu_value_type(3.1415) * ublas_A; \
340 viennacl::copy(ublas_A, vcl_A); \
341 viennacl::copy(ublas_B, vcl_B); \
343 for (std::size_t i=0; i<ublas_C.size1(); ++i) \
344 for (std::size_t j=0; j<ublas_C.size2(); ++j) \
345 ublas_C(i,j) = static_cast<cpu_value_type>(OPNAME(ublas_A(i,j))); \
346 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_##OPNAME(vcl_A)); \
347 viennacl::scheduler::execute(my_statement); \
348 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS) \
349 return EXIT_FAILURE; \
352 for (std::size_t i=0; i<ublas_C.size1(); ++i) \
353 for (std::size_t j=0; j<ublas_C.size2(); ++j) \
354 ublas_C(i,j) = static_cast<cpu_value_type>(OPNAME(ublas_A(i,j) / cpu_value_type(2))); \
355 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_##OPNAME(vcl_A / cpu_value_type(2))); \
356 viennacl::scheduler::execute(my_statement); \
357 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS) \
358 return EXIT_FAILURE; \
376 #undef GENERATE_UNARY_OP_TEST
378 if (ublas_C.size1() == ublas_C.size2())
380 std::cout <<
"z = element_fabs(x - trans(y))... ";
382 for (std::size_t i=0; i<ublas_C.size1(); ++i)
383 for (std::size_t j=0; j<ublas_C.size2(); ++j)
384 ublas_C(i,j) = std::fabs(ublas_A(i,j) - ublas_B(j,i));
392 std::cout <<
"z = element_fabs(trans(x) - (y))... ";
394 for (std::size_t i=0; i<ublas_C.size1(); ++i)
395 for (std::size_t j=0; j<ublas_C.size2(); ++j)
396 ublas_C(i,j) = std::fabs(ublas_A(j,i) - ublas_B(i,j));
404 std::cout <<
"z = element_fabs(trans(x) - trans(y))... ";
406 for (std::size_t i=0; i<ublas_C.size1(); ++i)
407 for (std::size_t j=0; j<ublas_C.size2(); ++j)
408 ublas_C(i,j) = std::fabs(ublas_A(j,i) - ublas_B(j,i));
416 std::cout <<
"z += trans(x)... ";
418 for (std::size_t i=0; i<ublas_C.size1(); ++i)
419 for (std::size_t j=0; j<ublas_C.size2(); ++j)
420 ublas_C(i,j) += ublas_A(j,i);
428 std::cout <<
"z -= trans(x)... ";
430 for (std::size_t i=0; i<ublas_C.size1(); ++i)
431 for (std::size_t j=0; j<ublas_C.size2(); ++j)
432 ublas_C(i,j) -= ublas_A(j,i);
442 std::cout << std::endl;
443 std::cout <<
"----------------------------------------------" << std::endl;
444 std::cout << std::endl;
453 template<
typename T,
typename ScalarType>
457 typedef boost::numeric::ublas::matrix<ScalarType> MatrixType;
461 std::size_t dim_rows = 131;
462 std::size_t dim_cols = 33;
467 MatrixType ublas_A(dim_rows, dim_cols);
468 MatrixType ublas_B(dim_rows, dim_cols);
469 MatrixType ublas_C(dim_rows, dim_cols);
471 for (std::size_t i=0; i<ublas_A.size1(); ++i)
472 for (std::size_t j=0; j<ublas_A.size2(); ++j)
474 ublas_A(i,j) =
ScalarType((i+2) + (j+1)*(i+2));
475 ublas_B(i,j) =
ScalarType((j+2) + (j+1)*(j+2));
476 ublas_C(i,j) =
ScalarType((i+1) + (i+1)*(i+2));
479 MatrixType ublas_A_large(4 * dim_rows, 4 * dim_cols);
480 for (std::size_t i=0; i<ublas_A_large.size1(); ++i)
481 for (std::size_t j=0; j<ublas_A_large.size2(); ++j)
482 ublas_A_large(i,j) =
ScalarType(i * ublas_A_large.size2() + j);
485 VCLMatrixType vcl_A_full(4 * dim_rows, 4 * dim_cols);
486 VCLMatrixType vcl_B_full(4 * dim_rows, 4 * dim_cols);
487 VCLMatrixType vcl_C_full(4 * dim_rows, 4 * dim_cols);
496 VCLMatrixType vcl_A(dim_rows, dim_cols);
510 VCLMatrixType vcl_B(dim_rows, dim_cols);
524 VCLMatrixType vcl_C(dim_rows, dim_cols);
547 std::cout << std::endl;
548 std::cout <<
"//" << std::endl;
549 std::cout <<
"////////// Test: Copy CTOR //////////" << std::endl;
550 std::cout <<
"//" << std::endl;
553 std::cout <<
"Testing matrix created from range... ";
554 VCLMatrixType vcl_temp = vcl_range_A;
556 std::cout <<
"PASSED!" << std::endl;
559 std::cout <<
"ublas_A: " << ublas_A << std::endl;
560 std::cout <<
"vcl_temp: " << vcl_temp << std::endl;
561 std::cout <<
"vcl_range_A: " << vcl_range_A << std::endl;
562 std::cout <<
"vcl_A: " << vcl_A << std::endl;
563 std::cout << std::endl <<
"TEST failed!" << std::endl;
567 std::cout <<
"Testing matrix created from slice... ";
568 VCLMatrixType vcl_temp2 = vcl_range_B;
570 std::cout <<
"PASSED!" << std::endl;
573 std::cout << std::endl <<
"TEST failed!" << std::endl;
578 std::cout <<
"//" << std::endl;
579 std::cout <<
"////////// Test: Initializer for matrix type //////////" << std::endl;
580 std::cout <<
"//" << std::endl;
583 ublas::matrix<ScalarType> ublas_dummy1 = ublas::identity_matrix<ScalarType>(ublas_A.size1());
584 ublas::matrix<ScalarType> ublas_dummy2 = ublas::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
585 ublas::matrix<ScalarType> ublas_dummy3 = ublas::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
591 std::cout <<
"Testing initializer CTOR... ";
596 std::cout <<
"PASSED!" << std::endl;
599 std::cout << std::endl <<
"TEST failed!" << std::endl;
603 ublas_dummy1 = ublas::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
604 ublas_dummy2 = ublas::identity_matrix<ScalarType>(ublas_A.size1());
605 ublas_dummy3 = ublas::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
611 std::cout <<
"Testing initializer assignment... ";
616 std::cout <<
"PASSED!" << std::endl;
619 std::cout << std::endl <<
"TEST failed!" << std::endl;
630 std::cout <<
"Testing A=matrix, B=matrix, C=matrix ..." << std::endl;
635 ublas_A, ublas_B, ublas_C,
636 vcl_A, vcl_B, vcl_C) != EXIT_SUCCESS)
641 std::cout <<
"Testing A=matrix, B=matrix, C=range ..." << std::endl;
646 ublas_A, ublas_B, ublas_C,
647 vcl_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
652 std::cout <<
"Testing A=matrix, B=matrix, C=slice ..." << std::endl;
657 ublas_A, ublas_B, ublas_C,
658 vcl_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
663 std::cout <<
"Testing A=matrix, B=range, C=matrix ..." << std::endl;
668 ublas_A, ublas_B, ublas_C,
669 vcl_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
674 std::cout <<
"Testing A=matrix, B=range, C=range ..." << std::endl;
679 ublas_A, ublas_B, ublas_C,
680 vcl_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
685 std::cout <<
"Testing A=matrix, B=range, C=slice ..." << std::endl;
690 ublas_A, ublas_B, ublas_C,
691 vcl_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
697 std::cout <<
"Testing A=matrix, B=slice, C=matrix ..." << std::endl;
702 ublas_A, ublas_B, ublas_C,
703 vcl_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
708 std::cout <<
"Testing A=matrix, B=slice, C=range ..." << std::endl;
713 ublas_A, ublas_B, ublas_C,
714 vcl_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
719 std::cout <<
"Testing A=matrix, B=slice, C=slice ..." << std::endl;
724 ublas_A, ublas_B, ublas_C,
725 vcl_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
733 std::cout <<
"Testing A=range, B=matrix, C=matrix ..." << std::endl;
738 ublas_A, ublas_B, ublas_C,
739 vcl_range_A, vcl_B, vcl_C) != EXIT_SUCCESS)
744 std::cout <<
"Testing A=range, B=matrix, C=range ..." << std::endl;
749 ublas_A, ublas_B, ublas_C,
750 vcl_range_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
755 std::cout <<
"Testing A=range, B=matrix, C=slice ..." << std::endl;
760 ublas_A, ublas_B, ublas_C,
761 vcl_range_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
768 std::cout <<
"Testing A=range, B=range, C=matrix ..." << std::endl;
773 ublas_A, ublas_B, ublas_C,
774 vcl_range_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
779 std::cout <<
"Testing A=range, B=range, C=range ..." << std::endl;
784 ublas_A, ublas_B, ublas_C,
785 vcl_range_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
790 std::cout <<
"Testing A=range, B=range, C=slice ..." << std::endl;
795 ublas_A, ublas_B, ublas_C,
796 vcl_range_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
803 std::cout <<
"Testing A=range, B=slice, C=matrix ..." << std::endl;
808 ublas_A, ublas_B, ublas_C,
809 vcl_range_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
814 std::cout <<
"Testing A=range, B=slice, C=range ..." << std::endl;
819 ublas_A, ublas_B, ublas_C,
820 vcl_range_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
825 std::cout <<
"Testing A=range, B=slice, C=slice ..." << std::endl;
830 ublas_A, ublas_B, ublas_C,
831 vcl_range_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
838 std::cout <<
"Testing A=slice, B=matrix, C=matrix ..." << std::endl;
843 ublas_A, ublas_B, ublas_C,
844 vcl_slice_A, vcl_B, vcl_C) != EXIT_SUCCESS)
849 std::cout <<
"Testing A=slice, B=matrix, C=range ..." << std::endl;
854 ublas_A, ublas_B, ublas_C,
855 vcl_slice_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
860 std::cout <<
"Testing A=slice, B=matrix, C=slice ..." << std::endl;
865 ublas_A, ublas_B, ublas_C,
866 vcl_slice_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
873 std::cout <<
"Testing A=slice, B=range, C=matrix ..." << std::endl;
878 ublas_A, ublas_B, ublas_C,
879 vcl_slice_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
884 std::cout <<
"Testing A=slice, B=range, C=range ..." << std::endl;
889 ublas_A, ublas_B, ublas_C,
890 vcl_slice_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
895 std::cout <<
"Testing A=slice, B=range, C=slice ..." << std::endl;
900 ublas_A, ublas_B, ublas_C,
901 vcl_slice_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
908 std::cout <<
"Testing A=slice, B=slice, C=matrix ..." << std::endl;
913 ublas_A, ublas_B, ublas_C,
914 vcl_slice_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
919 std::cout <<
"Testing A=slice, B=slice, C=range ..." << std::endl;
924 ublas_A, ublas_B, ublas_C,
925 vcl_slice_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
930 std::cout <<
"Testing A=slice, B=slice, C=slice ..." << std::endl;
935 ublas_A, ublas_B, ublas_C,
936 vcl_slice_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
947 std::cout << std::endl;
948 std::cout <<
"----------------------------------------------" << std::endl;
949 std::cout <<
"----------------------------------------------" << std::endl;
950 std::cout <<
"## Test :: Matrix Range" << std::endl;
951 std::cout <<
"----------------------------------------------" << std::endl;
952 std::cout <<
"----------------------------------------------" << std::endl;
953 std::cout << std::endl;
955 double epsilon = 1e-4;
956 std::cout <<
"# Testing setup:" << std::endl;
957 std::cout <<
" eps: " << epsilon << std::endl;
958 std::cout <<
" numeric: float" << std::endl;
959 std::cout <<
" --- row-major ---" << std::endl;
960 if (run_test<viennacl::row_major, float>(epsilon) != EXIT_SUCCESS)
962 std::cout <<
" --- column-major ---" << std::endl;
963 if (run_test<viennacl::column_major, float>(epsilon) != EXIT_SUCCESS)
967 #ifdef VIENNACL_WITH_OPENCL
972 std::cout <<
"# Testing setup:" << std::endl;
973 std::cout <<
" eps: " << epsilon << std::endl;
974 std::cout <<
" numeric: double" << std::endl;
976 if (run_test<viennacl::row_major, double>(epsilon) != EXIT_SUCCESS)
978 if (run_test<viennacl::column_major, double>(epsilon) != EXIT_SUCCESS)
982 std::cout << std::endl;
983 std::cout <<
"------- Test completed --------" << std::endl;
984 std::cout << std::endl;
viennacl::vector_expression< const vector_base< T >, const vector_base< T >, op_element_binary< op_div > > element_div(vector_base< T > const &v1, vector_base< T > const &v2)
int main(int, const char **)
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Class for representing strided submatrices of a bigger matrix A.
std::vector< std::vector< NumericT > > trans(std::vector< std::vector< NumericT > > const &A)
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Implementation of the dense matrix class.
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)
Represents a vector consisting of 1 at a given index and zeros otherwise. To be used as an initialize...
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.
A tag class representing inplace addition.
bool check_for_equality(MatrixType const &ublas_A, VCLMatrixType const &vcl_A, double epsilon)
Represents a vector consisting of scalars 's' only, i.e. v[i] = s for all i. To be used as an initial...
#define GENERATE_UNARY_OP_TEST(OPNAME)
Represents a vector consisting of zeros only. To be used as an initializer for viennacl::vector, vector_range, or vector_slize only.
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
int run_test(double epsilon, UBLASMatrixType &ublas_A, UBLASMatrixType &ublas_B, UBLASMatrixType &ublas_C, ViennaCLMatrixType1 &vcl_A, ViennaCLMatrixType2 &vcl_B, ViennaCLMatrixType3 vcl_C)
Proxy classes for vectors.
Proxy classes for matrices.
A tag class representing inplace subtraction.
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 range class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
viennacl::vector_expression< const vector_base< T >, const vector_base< T >, op_element_binary< op_prod > > element_prod(vector_base< T > const &v1, vector_base< T > const &v2)
Provides the datastructures for dealing with a single statement such as 'x = y + z;'.
Class for representing non-strided submatrices of a bigger matrix A.
The main class for representing a statement such as x = inner_prod(y,z); at runtime.
A slice class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Implementation of the ViennaCL scalar class.