TriangularMatrixVector_BLAS.h
Go to the documentation of this file.
1 /*
2  Copyright (c) 2011, Intel Corporation. All rights reserved.
3 
4  Redistribution and use in source and binary forms, with or without modification,
5  are permitted provided that the following conditions are met:
6 
7  * Redistributions of source code must retain the above copyright notice, this
8  list of conditions and the following disclaimer.
9  * Redistributions in binary form must reproduce the above copyright notice,
10  this list of conditions and the following disclaimer in the documentation
11  and/or other materials provided with the distribution.
12  * Neither the name of Intel Corporation nor the names of its contributors may
13  be used to endorse or promote products derived from this software without
14  specific prior written permission.
15 
16  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
20  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
23  ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 
27  ********************************************************************************
28  * Content : Eigen bindings to BLAS F77
29  * Triangular matrix-vector product functionality based on ?TRMV.
30  ********************************************************************************
31 */
32 
33 #ifndef EIGEN_TRIANGULAR_MATRIX_VECTOR_BLAS_H
34 #define EIGEN_TRIANGULAR_MATRIX_VECTOR_BLAS_H
35 
36 namespace Eigen {
37 
38 namespace internal {
39 
40 /**********************************************************************
41 * This file implements triangular matrix-vector multiplication using BLAS
42 **********************************************************************/
43 
44 // trmv/hemv specialization
45 
46 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder>
48  triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,StorageOrder,BuiltIn> {};
49 
50 #define EIGEN_BLAS_TRMV_SPECIALIZE(Scalar) \
51 template<typename Index, int Mode, bool ConjLhs, bool ConjRhs> \
52 struct triangular_matrix_vector_product<Index,Mode,Scalar,ConjLhs,Scalar,ConjRhs,ColMajor,Specialized> { \
53  static void run(Index _rows, Index _cols, const Scalar* _lhs, Index lhsStride, \
54  const Scalar* _rhs, Index rhsIncr, Scalar* _res, Index resIncr, Scalar alpha) { \
55  triangular_matrix_vector_product_trmv<Index,Mode,Scalar,ConjLhs,Scalar,ConjRhs,ColMajor>::run( \
56  _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \
57  } \
58 }; \
59 template<typename Index, int Mode, bool ConjLhs, bool ConjRhs> \
60 struct triangular_matrix_vector_product<Index,Mode,Scalar,ConjLhs,Scalar,ConjRhs,RowMajor,Specialized> { \
61  static void run(Index _rows, Index _cols, const Scalar* _lhs, Index lhsStride, \
62  const Scalar* _rhs, Index rhsIncr, Scalar* _res, Index resIncr, Scalar alpha) { \
63  triangular_matrix_vector_product_trmv<Index,Mode,Scalar,ConjLhs,Scalar,ConjRhs,RowMajor>::run( \
64  _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \
65  } \
66 };
67 
72 
73 // implements col-major: res += alpha * op(triangular) * vector
74 #define EIGEN_BLAS_TRMV_CM(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX, BLASPOSTFIX) \
75 template<typename Index, int Mode, bool ConjLhs, bool ConjRhs> \
76 struct triangular_matrix_vector_product_trmv<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,ColMajor> { \
77  enum { \
78  IsLower = (Mode&Lower) == Lower, \
79  SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \
80  IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \
81  IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \
82  LowUp = IsLower ? Lower : Upper \
83  }; \
84  static void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \
85  const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \
86  { \
87  if (ConjLhs || IsZeroDiag) { \
88  triangular_matrix_vector_product<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,ColMajor,BuiltIn>::run( \
89  _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \
90  return; \
91  }\
92  Index size = (std::min)(_rows,_cols); \
93  Index rows = IsLower ? _rows : size; \
94  Index cols = IsLower ? size : _cols; \
95 \
96  typedef VectorX##EIGPREFIX VectorRhs; \
97  EIGTYPE *x, *y;\
98 \
99 /* Set x*/ \
100  Map<const VectorRhs, 0, InnerStride<> > rhs(_rhs,cols,InnerStride<>(rhsIncr)); \
101  VectorRhs x_tmp; \
102  if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \
103  x = x_tmp.data(); \
104 \
105 /* Square part handling */\
106 \
107  char trans, uplo, diag; \
108  BlasIndex m, n, lda, incx, incy; \
109  EIGTYPE const *a; \
110  EIGTYPE beta(1); \
111 \
112 /* Set m, n */ \
113  n = convert_index<BlasIndex>(size); \
114  lda = convert_index<BlasIndex>(lhsStride); \
115  incx = 1; \
116  incy = convert_index<BlasIndex>(resIncr); \
117 \
118 /* Set uplo, trans and diag*/ \
119  trans = 'N'; \
120  uplo = IsLower ? 'L' : 'U'; \
121  diag = IsUnitDiag ? 'U' : 'N'; \
122 \
123 /* call ?TRMV*/ \
124  BLASPREFIX##trmv##BLASPOSTFIX(&uplo, &trans, &diag, &n, (const BLASTYPE*)_lhs, &lda, (BLASTYPE*)x, &incx); \
125 \
126 /* Add op(a_tr)rhs into res*/ \
127  BLASPREFIX##axpy##BLASPOSTFIX(&n, (const BLASTYPE*)&numext::real_ref(alpha),(const BLASTYPE*)x, &incx, (BLASTYPE*)_res, &incy); \
128 /* Non-square case - doesn't fit to BLAS ?TRMV. Fall to default triangular product*/ \
129  if (size<(std::max)(rows,cols)) { \
130  if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \
131  x = x_tmp.data(); \
132  if (size<rows) { \
133  y = _res + size*resIncr; \
134  a = _lhs + size; \
135  m = convert_index<BlasIndex>(rows-size); \
136  n = convert_index<BlasIndex>(size); \
137  } \
138  else { \
139  x += size; \
140  y = _res; \
141  a = _lhs + size*lda; \
142  m = convert_index<BlasIndex>(size); \
143  n = convert_index<BlasIndex>(cols-size); \
144  } \
145  BLASPREFIX##gemv##BLASPOSTFIX(&trans, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)x, &incx, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)y, &incy); \
146  } \
147  } \
148 };
149 
150 #ifdef EIGEN_USE_MKL
151 EIGEN_BLAS_TRMV_CM(double, double, d, d,)
152 EIGEN_BLAS_TRMV_CM(dcomplex, MKL_Complex16, cd, z,)
153 EIGEN_BLAS_TRMV_CM(float, float, f, s,)
154 EIGEN_BLAS_TRMV_CM(scomplex, MKL_Complex8, cf, c,)
155 #else
156 EIGEN_BLAS_TRMV_CM(double, double, d, d, _)
157 EIGEN_BLAS_TRMV_CM(dcomplex, double, cd, z, _)
158 EIGEN_BLAS_TRMV_CM(float, float, f, s, _)
159 EIGEN_BLAS_TRMV_CM(scomplex, float, cf, c, _)
160 #endif
161 
162 // implements row-major: res += alpha * op(triangular) * vector
163 #define EIGEN_BLAS_TRMV_RM(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX, BLASPOSTFIX) \
164 template<typename Index, int Mode, bool ConjLhs, bool ConjRhs> \
165 struct triangular_matrix_vector_product_trmv<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,RowMajor> { \
166  enum { \
167  IsLower = (Mode&Lower) == Lower, \
168  SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \
169  IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \
170  IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \
171  LowUp = IsLower ? Lower : Upper \
172  }; \
173  static void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \
174  const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \
175  { \
176  if (IsZeroDiag) { \
177  triangular_matrix_vector_product<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,RowMajor,BuiltIn>::run( \
178  _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \
179  return; \
180  }\
181  Index size = (std::min)(_rows,_cols); \
182  Index rows = IsLower ? _rows : size; \
183  Index cols = IsLower ? size : _cols; \
184 \
185  typedef VectorX##EIGPREFIX VectorRhs; \
186  EIGTYPE *x, *y;\
187 \
188 /* Set x*/ \
189  Map<const VectorRhs, 0, InnerStride<> > rhs(_rhs,cols,InnerStride<>(rhsIncr)); \
190  VectorRhs x_tmp; \
191  if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \
192  x = x_tmp.data(); \
193 \
194 /* Square part handling */\
195 \
196  char trans, uplo, diag; \
197  BlasIndex m, n, lda, incx, incy; \
198  EIGTYPE const *a; \
199  EIGTYPE beta(1); \
200 \
201 /* Set m, n */ \
202  n = convert_index<BlasIndex>(size); \
203  lda = convert_index<BlasIndex>(lhsStride); \
204  incx = 1; \
205  incy = convert_index<BlasIndex>(resIncr); \
206 \
207 /* Set uplo, trans and diag*/ \
208  trans = ConjLhs ? 'C' : 'T'; \
209  uplo = IsLower ? 'U' : 'L'; \
210  diag = IsUnitDiag ? 'U' : 'N'; \
211 \
212 /* call ?TRMV*/ \
213  BLASPREFIX##trmv##BLASPOSTFIX(&uplo, &trans, &diag, &n, (const BLASTYPE*)_lhs, &lda, (BLASTYPE*)x, &incx); \
214 \
215 /* Add op(a_tr)rhs into res*/ \
216  BLASPREFIX##axpy##BLASPOSTFIX(&n, (const BLASTYPE*)&numext::real_ref(alpha),(const BLASTYPE*)x, &incx, (BLASTYPE*)_res, &incy); \
217 /* Non-square case - doesn't fit to BLAS ?TRMV. Fall to default triangular product*/ \
218  if (size<(std::max)(rows,cols)) { \
219  if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \
220  x = x_tmp.data(); \
221  if (size<rows) { \
222  y = _res + size*resIncr; \
223  a = _lhs + size*lda; \
224  m = convert_index<BlasIndex>(rows-size); \
225  n = convert_index<BlasIndex>(size); \
226  } \
227  else { \
228  x += size; \
229  y = _res; \
230  a = _lhs + size; \
231  m = convert_index<BlasIndex>(size); \
232  n = convert_index<BlasIndex>(cols-size); \
233  } \
234  BLASPREFIX##gemv##BLASPOSTFIX(&trans, &n, &m, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)x, &incx, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)y, &incy); \
235  } \
236  } \
237 };
238 
239 #ifdef EIGEN_USE_MKL
240 EIGEN_BLAS_TRMV_RM(double, double, d, d,)
241 EIGEN_BLAS_TRMV_RM(dcomplex, MKL_Complex16, cd, z,)
242 EIGEN_BLAS_TRMV_RM(float, float, f, s,)
243 EIGEN_BLAS_TRMV_RM(scomplex, MKL_Complex8, cf, c,)
244 #else
245 EIGEN_BLAS_TRMV_RM(double, double, d, d,_)
246 EIGEN_BLAS_TRMV_RM(dcomplex, double, cd, z,_)
247 EIGEN_BLAS_TRMV_RM(float, float, f, s,_)
248 EIGEN_BLAS_TRMV_RM(scomplex, float, cf, c,_)
249 #endif
250 
251 } // end namespase internal
252 
253 } // end namespace Eigen
254 
255 #endif // EIGEN_TRIANGULAR_MATRIX_VECTOR_BLAS_H
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
std::complex< float > scomplex
Definition: MKL_support.h:126
std::complex< double > dcomplex
Definition: MKL_support.h:125
#define EIGEN_BLAS_TRMV_SPECIALIZE(Scalar)
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
#define EIGEN_BLAS_TRMV_CM(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX, BLASPOSTFIX)
RealScalar s
constexpr descr< N - 1 > _(char const (&text)[N])
Definition: descr.h:109
#define EIGEN_BLAS_TRMV_RM(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX, BLASPOSTFIX)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:40:33