gemv.hpp
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief -
6  *
7  * \author O. Krause
8  * \date 2010
9  *
10  *
11  * \par Copyright 1995-2015 Shark Development Team
12  *
13  * <BR><HR>
14  * This file is part of Shark.
15  * <http://image.diku.dk/shark/>
16  *
17  * Shark is free software: you can redistribute it and/or modify
18  * it under the terms of the GNU Lesser General Public License as published
19  * by the Free Software Foundation, either version 3 of the License, or
20  * (at your option) any later version.
21  *
22  * Shark is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU Lesser General Public License for more details.
26  *
27  * You should have received a copy of the GNU Lesser General Public License
28  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
29  *
30  */
31 //===========================================================================
32 #ifndef SHARK_LINALG_BLAS_KERNELS_CBLAS_GEMV_HPP
33 #define SHARK_LINALG_BLAS_KERNELS_CBLAS_GEMV_HPP
34 
35 #include "cblas_inc.hpp"
36 
37 namespace shark {namespace blas {namespace bindings {
38 
39 inline void gemv(CBLAS_ORDER const Order,
40  CBLAS_TRANSPOSE const TransA, int const M, int const N,
41  double alpha, float const *A, int const lda,
42  float const *X, int const incX,
43  double beta, float *Y, int const incY
44 ) {
45  cblas_sgemv(Order, TransA, M, N, alpha, A, lda,
46  X, incX,
47  beta, Y, incY);
48 }
49 
50 inline void gemv(CBLAS_ORDER const Order,
51  CBLAS_TRANSPOSE const TransA, int const M, int const N,
52  double alpha, double const *A, int const lda,
53  double const *X, int const incX,
54  double beta, double *Y, int const incY
55 ) {
56  cblas_dgemv(Order, TransA, M, N, alpha, A, lda,
57  X, incX,
58  beta, Y, incY);
59 }
60 
61 inline void gemv(CBLAS_ORDER const Order,
62  CBLAS_TRANSPOSE const TransA, int const M, int const N,
63  double alpha,
64  std::complex<float> const *A, int const lda,
65  std::complex<float> const *X, int const incX,
66  double beta,
67  std::complex<float> *Y, int const incY
68 ) {
69  std::complex<float> alphaArg(alpha,0);
70  std::complex<float> betaArg(beta,0);
71  cblas_cgemv(Order, TransA, M, N,
72  reinterpret_cast<cblas_float_complex_type const *>(&alphaArg),
73  reinterpret_cast<cblas_float_complex_type const *>(A), lda,
74  reinterpret_cast<cblas_float_complex_type const *>(X), incX,
75  reinterpret_cast<cblas_float_complex_type const *>(&betaArg),
76  reinterpret_cast<cblas_float_complex_type *>(Y), incY);
77 }
78 
79 inline void gemv(CBLAS_ORDER const Order,
80  CBLAS_TRANSPOSE const TransA, int const M, int const N,
81  double alpha,
82  std::complex<double> const *A, int const lda,
83  std::complex<double> const *X, int const incX,
84  double beta,
85  std::complex<double> *Y, int const incY
86 ) {
87  std::complex<double> alphaArg(alpha,0);
88  std::complex<double> betaArg(beta,0);
89  cblas_zgemv(Order, TransA, M, N,
90  reinterpret_cast<cblas_double_complex_type const *>(&alphaArg),
91  reinterpret_cast<cblas_double_complex_type const *>(A), lda,
92  reinterpret_cast<cblas_double_complex_type const *>(X), incX,
93  reinterpret_cast<cblas_double_complex_type const *>(&betaArg),
94  reinterpret_cast<cblas_double_complex_type *>(Y), incY);
95 }
96 
97 
98 // y <- alpha * op (A) * x + beta * y
99 // op (A) == A || A^T || A^H
100 template <typename MatrA, typename VectorX, typename VectorY>
101 void gemv(
102  matrix_expression<MatrA> const &A,
103  vector_expression<VectorX> const &x,
104  vector_expression<VectorY> &y,
105  typename VectorY::value_type alpha,
106  boost::mpl::true_
107 ){
108  std::size_t m = A().size1();
109  std::size_t n = A().size2();
110 
111  SIZE_CHECK(x().size() == A().size2());
112  SIZE_CHECK(y().size() == A().size1());
113 
114  CBLAS_ORDER const stor_ord= (CBLAS_ORDER)storage_order<typename MatrA::orientation>::value;
115 
116  gemv(stor_ord, CblasNoTrans, (int)m, (int)n, alpha,
117  traits::storage(A),
118  traits::leading_dimension(A),
119  traits::storage(x),
120  traits::stride(x),
121  typename VectorY::value_type(1),
122  traits::storage(y),
123  traits::stride(y)
124  );
125 }
126 
127 template<class Storage1, class Storage2, class Storage3, class T1, class T2, class T3>
128 struct optimized_gemv_detail{
129  typedef boost::mpl::false_ type;
130 };
131 template<>
132 struct optimized_gemv_detail<
133  dense_tag, dense_tag, dense_tag,
134  double, double, double
135 >{
136  typedef boost::mpl::true_ type;
137 };
138 template<>
139 struct optimized_gemv_detail<
140  dense_tag, dense_tag, dense_tag,
141  float, float, float
142 >{
143  typedef boost::mpl::true_ type;
144 };
145 
146 template<>
147 struct optimized_gemv_detail<
148  dense_tag, dense_tag, dense_tag,
149  std::complex<double>, std::complex<double>, std::complex<double>
150 >{
151  typedef boost::mpl::true_ type;
152 };
153 template<>
154 struct optimized_gemv_detail<
155  dense_tag, dense_tag, dense_tag,
156  std::complex<float>, std::complex<float>, std::complex<float>
157 >{
158  typedef boost::mpl::true_ type;
159 };
160 
161 template<class M1, class M2, class M3>
162 struct has_optimized_gemv
163 : public optimized_gemv_detail<
164  typename M1::storage_category,
165  typename M2::storage_category,
166  typename M3::storage_category,
167  typename M1::value_type,
168  typename M2::value_type,
169  typename M3::value_type
170 >{};
171 
172 }}}
173 #endif