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 /** \example sparse.cpp
19 *
20 * This tutorial demonstrates the use of sparse matrices.
21 * The primary operation for sparse matrices in ViennaCL is the sparse matrix-vector product.
22 *
23 * We start with including the respective headers:
24 **/
25
26 // system headers
27 #include <iostream>
28
29 // ublas headers
30 #include <boost/numeric/ublas/matrix_sparse.hpp>
31 #include <boost/numeric/ublas/matrix.hpp>
32 #include <boost/numeric/ublas/operation.hpp>
33 #include <boost/numeric/ublas/operation_sparse.hpp>
34 #include <boost/numeric/ublas/io.hpp>
35
36 // Must be set if you want to use ViennaCL algorithms on ublas objects
37 #define VIENNACL_WITH_UBLAS 1
38
39 // ViennaCL includes
40 #include "viennacl/scalar.hpp"
41 #include "viennacl/vector.hpp"
42 #include "viennacl/compressed_matrix.hpp"
43 #include "viennacl/linalg/prod.hpp"
44
45
46 // Additional helper functions for this tutorial:
47 #include "vector-io.hpp"
48
49 // Shortcut for writing 'ublas::' instead of 'boost::numeric::ublas::'
50 using namespace boost::numeric;
51
52 /**
53 * We setup a sparse matrix in uBLAS and populate it with values.
54 * Then, the respective ViennaCL sparse matrix is created and initialized with data from the uBLAS matrix.
55 * After a direct manipulation of the ViennaCL matrix, matrix-vector products are computed with both matrices.
56 **/
main()57 int main()
58 {
59 typedef float ScalarType;
60
61 std::size_t size = 5;
62
63 /**
64 * Set up some ublas objects
65 **/
66 ublas::vector<ScalarType> rhs = ublas::scalar_vector<ScalarType>(size, ScalarType(size));
67 ublas::compressed_matrix<ScalarType> ublas_matrix(size, size);
68
69 ublas_matrix(0,0) = 2.0f; ublas_matrix(0,1) = -1.0f;
70 ublas_matrix(1,0) = -1.0f; ublas_matrix(1,1) = 2.0f; ublas_matrix(1,2) = -1.0f;
71 ublas_matrix(2,1) = -1.0f; ublas_matrix(2,2) = 2.0f; ublas_matrix(2,3) = -1.0f;
72 ublas_matrix(3,2) = -1.0f; ublas_matrix(3,3) = 2.0f; ublas_matrix(3,4) = -1.0f;
73 ublas_matrix(4,3) = -1.0f; ublas_matrix(4,4) = 2.0f;
74
75 std::cout << "ublas matrix: " << ublas_matrix << std::endl;
76
77 /**
78 * Set up some ViennaCL objects and initialize with data from uBLAS objects
79 **/
80 viennacl::vector<ScalarType> vcl_rhs(size);
81 viennacl::compressed_matrix<ScalarType> vcl_compressed_matrix(size, size);
82
83 viennacl::copy(rhs, vcl_rhs);
84 viennacl::copy(ublas_matrix, vcl_compressed_matrix);
85
86 // just get the data directly from the GPU and print it:
87 ublas::compressed_matrix<ScalarType> temp(size, size);
88 viennacl::copy(vcl_compressed_matrix, temp);
89 std::cout << "ViennaCL: " << temp << std::endl;
90
91 // now modify GPU data directly:
92 std::cout << "Modifying vcl_compressed_matrix a bit: " << std::endl;
93 vcl_compressed_matrix(0, 0) = 3.0f;
94 vcl_compressed_matrix(2, 3) = -3.0f;
95 vcl_compressed_matrix(4, 2) = -3.0f; //this is a new nonzero entry
96 vcl_compressed_matrix(4, 3) = -3.0f;
97
98 // and print it again:
99 viennacl::copy(vcl_compressed_matrix, temp);
100 std::cout << "ViennaCL matrix copied to uBLAS matrix: " << temp << std::endl;
101
102 /**
103 * Compute matrix-vector products and output the results (should match):
104 **/
105 std::cout << "ublas: " << ublas::prod(temp, rhs) << std::endl;
106 std::cout << "ViennaCL: " << viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs) << std::endl;
107
108 /**
109 * That's it. Print a success message and exit.
110 **/
111 std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
112
113 return EXIT_SUCCESS;
114 }
115
116