ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
blas2.cpp
Go to the documentation of this file.
1 /* =========================================================================
2  Copyright (c) 2010-2014, 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 // include necessary system headers
19 #include <iostream>
20 
21 #include "viennacl.hpp"
22 #include "viennacl_private.hpp"
23 
24 #include "init_vector.hpp"
25 #include "init_matrix.hpp"
26 
27 //include basic scalar and vector types of ViennaCL
28 #include "viennacl/scalar.hpp"
29 #include "viennacl/vector.hpp"
30 #include "viennacl/matrix.hpp"
31 #include "viennacl/linalg/direct_solve.hpp"
32 #include "viennacl/linalg/prod.hpp"
33 
34 // GEMV
35 
36 VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLgemv(ViennaCLHostScalar alpha, ViennaCLMatrix A, ViennaCLVector x, ViennaCLHostScalar beta, ViennaCLVector y)
37 {
38  viennacl::backend::mem_handle v1_handle;
39  viennacl::backend::mem_handle v2_handle;
40  viennacl::backend::mem_handle A_handle;
41 
42  if (init_vector(v1_handle, x) != ViennaCLSuccess)
43  return ViennaCLGenericFailure;
44 
45  if (init_vector(v2_handle, y) != ViennaCLSuccess)
46  return ViennaCLGenericFailure;
47 
48  if (init_matrix(A_handle, A) != ViennaCLSuccess)
49  return ViennaCLGenericFailure;
50 
51  switch (x->precision)
52  {
53  case ViennaCLFloat:
54  {
55  typedef viennacl::vector_base<float>::size_type size_type;
56  typedef viennacl::vector_base<float>::size_type difference_type;
57 
58  viennacl::vector_base<float> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
59  viennacl::vector_base<float> v2(v2_handle, size_type(y->size), size_type(y->offset), difference_type(y->inc));
60 
61  viennacl::matrix_base<float> mat(A_handle,
62  size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
63  size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
64  v2 *= beta->value_float;
65  if (A->trans == ViennaCLTrans)
66  v2 += alpha->value_float * viennacl::linalg::prod(viennacl::trans(mat), v1);
67  else
68  v2 += alpha->value_float * viennacl::linalg::prod(mat, v1);
69 
70  return ViennaCLSuccess;
71  }
72 
73  case ViennaCLDouble:
74  {
75  typedef viennacl::vector_base<double>::size_type size_type;
76  typedef viennacl::vector_base<double>::size_type difference_type;
77 
78  viennacl::vector_base<double> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
79  viennacl::vector_base<double> v2(v2_handle, size_type(y->size), size_type(y->offset), difference_type(y->inc));
80 
81  viennacl::matrix_base<double> mat(A_handle,
82  size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
83  size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
84  v2 *= beta->value_double;
85  if (A->trans == ViennaCLTrans)
86  v2 += alpha->value_double * viennacl::linalg::prod(viennacl::trans(mat), v1);
87  else
88  v2 += alpha->value_double * viennacl::linalg::prod(mat, v1);
89 
90  return ViennaCLSuccess;
91  }
92 
93  default:
94  return ViennaCLGenericFailure;
95  }
96 }
97 
98 
99 // xTRSV
100 
101 VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLtrsv(ViennaCLMatrix A, ViennaCLVector x, ViennaCLUplo uplo)
102 {
103  viennacl::backend::mem_handle v1_handle;
104  viennacl::backend::mem_handle A_handle;
105 
106  if (init_vector(v1_handle, x) != ViennaCLSuccess)
107  return ViennaCLGenericFailure;
108 
109  if (init_matrix(A_handle, A) != ViennaCLSuccess)
110  return ViennaCLGenericFailure;
111 
112  switch (x->precision)
113  {
114  case ViennaCLFloat:
115  {
116  typedef viennacl::vector_base<float>::size_type size_type;
117  typedef viennacl::vector_base<float>::size_type difference_type;
118 
119  viennacl::vector_base<float> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
120 
121  viennacl::matrix_base<float> mat(A_handle,
122  size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
123  size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
124  if (A->trans == ViennaCLTrans)
125  {
126  if (uplo == ViennaCLUpper)
127  viennacl::linalg::inplace_solve(viennacl::trans(mat), v1, viennacl::linalg::upper_tag());
128  else
129  viennacl::linalg::inplace_solve(viennacl::trans(mat), v1, viennacl::linalg::lower_tag());
130  }
131  else
132  {
133  if (uplo == ViennaCLUpper)
134  viennacl::linalg::inplace_solve(mat, v1, viennacl::linalg::upper_tag());
135  else
136  viennacl::linalg::inplace_solve(mat, v1, viennacl::linalg::lower_tag());
137  }
138 
139  return ViennaCLSuccess;
140  }
141  case ViennaCLDouble:
142  {
143  typedef viennacl::vector_base<double>::size_type size_type;
144  typedef viennacl::vector_base<double>::size_type difference_type;
145 
146  viennacl::vector_base<double> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
147 
148  viennacl::matrix_base<double> mat(A_handle,
149  size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
150  size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
151  if (A->trans == ViennaCLTrans)
152  {
153  if (uplo == ViennaCLUpper)
154  viennacl::linalg::inplace_solve(viennacl::trans(mat), v1, viennacl::linalg::upper_tag());
155  else
156  viennacl::linalg::inplace_solve(viennacl::trans(mat), v1, viennacl::linalg::lower_tag());
157  }
158  else
159  {
160  if (uplo == ViennaCLUpper)
161  viennacl::linalg::inplace_solve(mat, v1, viennacl::linalg::upper_tag());
162  else
163  viennacl::linalg::inplace_solve(mat, v1, viennacl::linalg::lower_tag());
164  }
165 
166  return ViennaCLSuccess;
167  }
168 
169  default:
170  return ViennaCLGenericFailure;
171  }
172 }
173 
174 
175 // xGER
176 
177 VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLger(ViennaCLHostScalar alpha, ViennaCLVector x, ViennaCLVector y, ViennaCLMatrix A)
178 {
179  viennacl::backend::mem_handle v1_handle;
180  viennacl::backend::mem_handle v2_handle;
181  viennacl::backend::mem_handle A_handle;
182 
183  if (init_vector(v1_handle, x) != ViennaCLSuccess)
184  return ViennaCLGenericFailure;
185 
186  if (init_vector(v2_handle, y) != ViennaCLSuccess)
187  return ViennaCLGenericFailure;
188 
189  if (init_matrix(A_handle, A) != ViennaCLSuccess)
190  return ViennaCLGenericFailure;
191 
192  switch (x->precision)
193  {
194  case ViennaCLFloat:
195  {
196  typedef viennacl::vector_base<float>::size_type size_type;
197  typedef viennacl::vector_base<float>::size_type difference_type;
198 
199  viennacl::vector_base<float> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
200  viennacl::vector_base<float> v2(v2_handle, size_type(y->size), size_type(y->offset), difference_type(y->inc));
201 
202  viennacl::matrix_base<float> mat(A_handle,
203  size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
204  size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
205 
206  mat += alpha->value_float * viennacl::linalg::outer_prod(v1, v2);
207 
208  return ViennaCLSuccess;
209  }
210  case ViennaCLDouble:
211  {
212  typedef viennacl::vector_base<double>::size_type size_type;
213  typedef viennacl::vector_base<double>::size_type difference_type;
214 
215  viennacl::vector_base<double> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
216  viennacl::vector_base<double> v2(v2_handle, size_type(y->size), size_type(y->offset), difference_type(y->inc));
217 
218  viennacl::matrix_base<double> mat(A_handle,
219  size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
220  size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
221 
222  mat += alpha->value_double * viennacl::linalg::outer_prod(v1, v2);
223 
224  return ViennaCLSuccess;
225  }
226  default:
227  return ViennaCLGenericFailure;
228  }
229 }
230 
231 
ViennaCLtrsv
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLtrsv(ViennaCLMatrix A, ViennaCLVector x, ViennaCLUplo uplo)
Definition: blas2.cpp:101
viennacl::linalg::inplace_solve
void inplace_solve(const matrix_base< NumericT > &A, matrix_base< NumericT > &B, SolverTagT)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:217
viennacl::trans
viennacl::enable_if< viennacl::is_any_sparse_matrix< M1 >::value, matrix_expression< const M1, const M1, op_trans > >::type trans(const M1 &mat)
Returns an expression template class representing a transposed matrix.
Definition: sparse_matrix_operations.hpp:374
ViennaCLVector_impl::precision
ViennaCLPrecision precision
Definition: viennacl_private.hpp:97
ViennaCLHostScalar_impl::value_float
float value_float
Definition: viennacl_private.hpp:72
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
ViennaCLMatrix_impl::start1
ViennaCLInt start1
Definition: viennacl_private.hpp:130
matrix.hpp
Implementation of the dense matrix class.
ViennaCLMatrix_impl::stride2
ViennaCLInt stride2
Definition: viennacl_private.hpp:136
viennacl::linalg::lower_tag
A tag class representing a lower triangular matrix.
Definition: forwards.h:849
ViennaCLgemv
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLgemv(ViennaCLHostScalar alpha, ViennaCLMatrix A, ViennaCLVector x, ViennaCLHostScalar beta, ViennaCLVector y)
Definition: blas2.cpp:36
viennacl::matrix_base< float >
ViennaCLStatus
ViennaCLStatus
Definition: viennacl.hpp:97
ViennaCLMatrix_impl::internal_size2
ViennaCLInt internal_size2
Definition: viennacl_private.hpp:137
ViennaCLVector_impl::size
ViennaCLInt size
Definition: viennacl_private.hpp:110
ViennaCLMatrix_impl::order
ViennaCLOrder order
Definition: viennacl_private.hpp:117
VIENNACL_EXPORTED_FUNCTION
#define VIENNACL_EXPORTED_FUNCTION
Definition: viennacl.hpp:40
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
ViennaCLMatrix_impl
Definition: viennacl_private.hpp:113
ViennaCLVector_impl::inc
ViennaCLInt inc
Definition: viennacl_private.hpp:109
ViennaCLVector_impl
Definition: viennacl_private.hpp:94
viennacl_private.hpp
init_matrix.hpp
ViennaCLMatrix_impl::stride1
ViennaCLInt stride1
Definition: viennacl_private.hpp:131
ViennaCLUplo
ViennaCLUplo
Definition: viennacl.hpp:75
viennacl::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
ViennaCLMatrix_impl::start2
ViennaCLInt start2
Definition: viennacl_private.hpp:135
ViennaCLMatrix_impl::size1
ViennaCLInt size1
Definition: viennacl_private.hpp:129
ViennaCLMatrix_impl::internal_size1
ViennaCLInt internal_size1
Definition: viennacl_private.hpp:132
viennacl::vector_base< float >
init_vector.hpp
direct_solve.hpp
Implementations of dense direct solvers are found here.
ViennaCLVector_impl::offset
ViennaCLInt offset
Definition: viennacl_private.hpp:108
viennacl::linalg::outer_prod
viennacl::matrix_expression< const vector_base< NumericT >, const vector_base< NumericT >, op_prod > outer_prod(const vector_base< NumericT > &vec1, const vector_base< NumericT > &vec2)
Returns a proxy class for the operation mat += vec1 * vec2^T, i.e. a rank 1 update.
Definition: matrix_operations.hpp:840
ViennaCLHostScalar_impl
Definition: viennacl_private.hpp:67
v2
viennacl::vector< int > v2
Definition: global_variables.cpp:61
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
ViennaCLMatrix_impl::size2
ViennaCLInt size2
Definition: viennacl_private.hpp:134
viennacl::backend::mem_handle
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
ViennaCLMatrix_impl::trans
ViennaCLTranspose trans
Definition: viennacl_private.hpp:118
ViennaCLger
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLger(ViennaCLHostScalar alpha, ViennaCLVector x, ViennaCLVector y, ViennaCLMatrix A)
Definition: blas2.cpp:177
scalar.hpp
Implementation of the ViennaCL scalar class.
ViennaCLHostScalar_impl::value_double
double value_double
Definition: viennacl_private.hpp:73
viennacl::vector_base::size_type
SizeT size_type
Definition: vector_def.hpp:112