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 matrix-range.cpp
19 *
20 * This tutorial explains the use of matrix ranges with simple BLAS level 1 and 2 operations.
21 *
22 * We start with including the necessary headers:
23 **/
24
25
26 // activate ublas support in ViennaCL
27 #define VIENNACL_WITH_UBLAS
28
29 // System headers
30 #include <iostream>
31 #include <string>
32
33
34 // ViennaCL headers
35 #include "viennacl/scalar.hpp"
36 #include "viennacl/matrix.hpp"
37 #include "viennacl/linalg/prod.hpp"
38 #include "viennacl/matrix_proxy.hpp"
39
40
41 // Boost headers
42 #include "boost/numeric/ublas/vector.hpp"
43 #include "boost/numeric/ublas/matrix.hpp"
44 #include "boost/numeric/ublas/matrix_proxy.hpp"
45 #include "boost/numeric/ublas/io.hpp"
46
47 /**
48 * In the main() routine we set up Boost.uBLAS as well as ViennaCL objects.
49 * A few standard operations on submatrices are performed by using the matrix_range<> view available in both libraries.
50 **/
main(int,const char **)51 int main (int, const char **)
52 {
53 // feel free to change this to 'double' if supported by your hardware
54 typedef float ScalarType;
55 typedef boost::numeric::ublas::matrix<ScalarType> MatrixType;
56
57 typedef viennacl::matrix<ScalarType, viennacl::row_major> VCLMatrixType;
58
59 /**
60 * Setup ublas objects and fill with data:
61 **/
62 std::size_t dim_large = 5;
63 std::size_t dim_small = 3;
64
65 MatrixType ublas_A(dim_large, dim_large);
66 MatrixType ublas_B(dim_small, dim_small);
67 MatrixType ublas_C(dim_large, dim_small);
68 MatrixType ublas_D(dim_small, dim_large);
69
70
71 for (std::size_t i=0; i<ublas_A.size1(); ++i)
72 for (std::size_t j=0; j<ublas_A.size2(); ++j)
73 ublas_A(i,j) = static_cast<ScalarType>((i+1) + (j+1)*(i+1));
74
75 for (std::size_t i=0; i<ublas_B.size1(); ++i)
76 for (std::size_t j=0; j<ublas_B.size2(); ++j)
77 ublas_B(i,j) = static_cast<ScalarType>((i+1) + (j+1)*(i+1));
78
79 for (std::size_t i=0; i<ublas_C.size1(); ++i)
80 for (std::size_t j=0; j<ublas_C.size2(); ++j)
81 ublas_C(i,j) = static_cast<ScalarType>((j+2) + (j+1)*(i+1));
82
83 for (std::size_t i=0; i<ublas_D.size1(); ++i)
84 for (std::size_t j=0; j<ublas_D.size2(); ++j)
85 ublas_D(i,j) = static_cast<ScalarType>((j+2) + (j+1)*(i+1));
86
87 /**
88 * Extract submatrices using the ranges in ublas
89 **/
90 boost::numeric::ublas::range ublas_r1(0, dim_small); //the first 'dim_small' entries
91 boost::numeric::ublas::range ublas_r2(dim_large - dim_small, dim_large); //the last 'dim_small' entries
92 boost::numeric::ublas::matrix_range<MatrixType> ublas_A_sub1(ublas_A, ublas_r1, ublas_r1); //upper left part of A
93 boost::numeric::ublas::matrix_range<MatrixType> ublas_A_sub2(ublas_A, ublas_r2, ublas_r2); //lower right part of A
94
95 boost::numeric::ublas::matrix_range<MatrixType> ublas_C_sub(ublas_C, ublas_r1, ublas_r1); //upper left part of C
96 boost::numeric::ublas::matrix_range<MatrixType> ublas_D_sub(ublas_D, ublas_r1, ublas_r1); //upper left part of D
97
98 /**
99 * Setup ViennaCL objects and copy data from uBLAS objects
100 **/
101 VCLMatrixType vcl_A(dim_large, dim_large);
102 VCLMatrixType vcl_B(dim_small, dim_small);
103 VCLMatrixType vcl_C(dim_large, dim_small);
104 VCLMatrixType vcl_D(dim_small, dim_large);
105
106 viennacl::copy(ublas_A, vcl_A);
107 viennacl::copy(ublas_B, vcl_B);
108 viennacl::copy(ublas_C, vcl_C);
109 viennacl::copy(ublas_D, vcl_D);
110
111 /**
112 * Extract submatrices using the ranges in ViennaCL. Similar to the code above for uBLAS.
113 **/
114 viennacl::range vcl_r1(0, dim_small); //the first 'dim_small' entries
115 viennacl::range vcl_r2(dim_large - dim_small, dim_large); //the last 'dim_small' entries
116 viennacl::matrix_range<VCLMatrixType> vcl_A_sub1(vcl_A, vcl_r1, vcl_r1); //upper left part of A
117 viennacl::matrix_range<VCLMatrixType> vcl_A_sub2(vcl_A, vcl_r2, vcl_r2); //lower right part of A
118
119 viennacl::matrix_range<VCLMatrixType> vcl_C_sub(vcl_C, vcl_r1, vcl_r1); //upper left part of C
120 viennacl::matrix_range<VCLMatrixType> vcl_D_sub(vcl_D, vcl_r1, vcl_r1); //upper left part of D
121
122 /**
123 * First use case: Copy from ublas to submatrices and back:
124 **/
125
126 ublas_A_sub1 = ublas_B;
127 viennacl::copy(ublas_B, vcl_A_sub1);
128 viennacl::copy(vcl_A_sub1, ublas_B);
129
130 /**
131 * Second use case: Addition of matrices.
132 **/
133
134 // range to range:
135 ublas_A_sub2 += ublas_A_sub2;
136 vcl_A_sub2 += vcl_A_sub2;
137
138 // range to matrix:
139 ublas_B += ublas_A_sub2;
140 vcl_B += vcl_A_sub2;
141
142
143 /**
144 * Third use case: Matrix range with matrix-matrix product:
145 **/
146 ublas_A_sub1 += prod(ublas_C_sub, ublas_D_sub);
147 vcl_A_sub1 += viennacl::linalg::prod(vcl_C_sub, vcl_D_sub);
148
149 /**
150 * Print result matrices:
151 **/
152 std::cout << "Result ublas: " << ublas_A << std::endl;
153 std::cout << "Result ViennaCL: " << vcl_A << std::endl;
154
155 /**
156 * That's it. Print a success message:
157 **/
158 std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
159
160 return EXIT_SUCCESS;
161 }
162
163