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 least-squares.cpp
19 *
20 * This tutorial shows how least Squares problems for matrices from ViennaCL or Boost.uBLAS can be solved solved.
21 *
22 * We start with including the respective header files:
23 **/
24
25 // activate ublas support in ViennaCL
26 #define VIENNACL_WITH_UBLAS
27
28 //
29 // include necessary system headers
30 //
31 #include <iostream>
32
33 // Boost headers
34 #include <boost/numeric/ublas/triangular.hpp>
35 #include <boost/numeric/ublas/vector.hpp>
36 #include <boost/numeric/ublas/vector_proxy.hpp>
37 #include <boost/numeric/ublas/matrix.hpp>
38 #include <boost/numeric/ublas/matrix_proxy.hpp>
39 #include <boost/numeric/ublas/lu.hpp>
40 #include <boost/numeric/ublas/io.hpp>
41
42
43 // ViennaCL headers
44 #include "viennacl/vector.hpp"
45 #include "viennacl/matrix.hpp"
46 #include "viennacl/matrix_proxy.hpp"
47 #include "viennacl/linalg/qr.hpp"
48 #include "viennacl/linalg/lu.hpp"
49 #include "viennacl/linalg/direct_solve.hpp"
50
51
52 /**
53 * The minimization problem of finding x such that \f$ \Vert Ax - b \Vert \f$ is solved as follows:
54 * - Compute the QR-factorization of A = QR.
55 * - Compute \f$ b' = Q^{\mathrm{T}} b \f$ for the equivalent minimization problem \f$ \Vert Rx - Q^{\mathrm{T}} b \f$.
56 * - Solve the triangular system \f$ \tilde{R} x = b' \f$, where \f$ \tilde{R} \f$ is the upper square matrix of R.
57 *
58 **/
main(int,const char **)59 int main (int, const char **)
60 {
61 typedef float ScalarType; //feel free to change this to 'double' if supported by your hardware
62
63 typedef boost::numeric::ublas::matrix<ScalarType> MatrixType;
64 typedef boost::numeric::ublas::vector<ScalarType> VectorType;
65 typedef viennacl::matrix<ScalarType, viennacl::column_major> VCLMatrixType;
66 typedef viennacl::vector<ScalarType> VCLVectorType;
67
68 /**
69 * Create vectors and matrices with data:
70 **/
71 VectorType ublas_b(4);
72 ublas_b(0) = -4;
73 ublas_b(1) = 2;
74 ublas_b(2) = 5;
75 ublas_b(3) = -1;
76
77 MatrixType ublas_A(4, 3);
78
79 ublas_A(0, 0) = 2; ublas_A(0, 1) = -1; ublas_A(0, 2) = 1;
80 ublas_A(1, 0) = 1; ublas_A(1, 1) = -5; ublas_A(1, 2) = 2;
81 ublas_A(2, 0) = -3; ublas_A(2, 1) = 1; ublas_A(2, 2) = -4;
82 ublas_A(3, 0) = 1; ublas_A(3, 1) = -1; ublas_A(3, 2) = 1;
83
84 /**
85 * Setup the matrix and vector with ViennaCL objects and copy the data from the uBLAS objects:
86 **/
87 VCLVectorType vcl_b(ublas_b.size());
88 VCLMatrixType vcl_A(ublas_A.size1(), ublas_A.size2());
89
90 viennacl::copy(ublas_b, vcl_b);
91 viennacl::copy(ublas_A, vcl_A);
92
93
94 /**
95 * <h2>Option 1: Using Boost.uBLAS</h2>
96 *
97 * The implementation in ViennaCL accepts both uBLAS and ViennaCL types.
98 * We start with a single-threaded implementation using Boost.uBLAS.
99 **/
100
101 std::cout << "--- Boost.uBLAS ---" << std::endl;
102 /**
103 * The first (and computationally most expensive) step is to compute the QR factorization of A.
104 * Since we do not need A later, we directly overwrite A with the householder reflectors and the upper triangular matrix R.
105 * The returned vector holds the scalar coefficients (betas) for the Householder reflections \f$ I - \beta v v^{\mathrm{T}} \f$
106 **/
107 std::vector<ScalarType> ublas_betas = viennacl::linalg::inplace_qr(ublas_A);
108
109 /**
110 * Compute the modified RHS of the minimization problem from the QR factorization, but do not form \f$ Q^{\mathrm{T}} \f$ explicitly:
111 * b' := Q^T b
112 **/
113 viennacl::linalg::inplace_qr_apply_trans_Q(ublas_A, ublas_betas, ublas_b);
114
115 /**
116 * Final step: triangular solve: Rx = b'', where b'' are the first three entries in b'
117 * We only need the upper left square part of A, which defines the upper triangular matrix R
118 **/
119 boost::numeric::ublas::range ublas_range(0, 3);
120 boost::numeric::ublas::matrix_range<MatrixType> ublas_R(ublas_A, ublas_range, ublas_range);
121 boost::numeric::ublas::vector_range<VectorType> ublas_b2(ublas_b, ublas_range);
122 boost::numeric::ublas::inplace_solve(ublas_R, ublas_b2, boost::numeric::ublas::upper_tag());
123
124 std::cout << "Result: " << ublas_b2 << std::endl;
125
126 /**
127 * <h2>Option 2: Use ViennaCL types</h2>
128 *
129 * ViennaCL is used for the computationally intensive BLAS 3 computations.
130 * Boost.uBLAS is used for the panel factorization on the host (CPU).
131 */
132
133 std::cout << "--- ViennaCL (hybrid implementation) ---" << std::endl;
134 std::vector<ScalarType> hybrid_betas = viennacl::linalg::inplace_qr(vcl_A);
135
136 /**
137 * compute modified RHS of the minimization problem: \f$ b' := Q^T b \f$
138 **/
139 viennacl::linalg::inplace_qr_apply_trans_Q(vcl_A, hybrid_betas, vcl_b);
140
141 /**
142 * Final step: triangular solve: Rx = b'.
143 * We only need the upper part of A such that R is a square matrix
144 **/
145 viennacl::range vcl_range(0, 3);
146 viennacl::matrix_range<VCLMatrixType> vcl_R(vcl_A, vcl_range, vcl_range);
147 viennacl::vector_range<VCLVectorType> vcl_b2(vcl_b, vcl_range);
148 viennacl::linalg::inplace_solve(vcl_R, vcl_b2, viennacl::linalg::upper_tag());
149
150 std::cout << "Result: " << vcl_b2 << std::endl;
151
152 /**
153 * That's it.
154 **/
155 std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
156
157 return EXIT_SUCCESS;
158 }
159
160