SelfadjointProduct.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SELFADJOINT_PRODUCT_H
11 #define EIGEN_SELFADJOINT_PRODUCT_H
12 
13 /**********************************************************************
14 * This file implements a self adjoint product: C += A A^T updating only
15 * half of the selfadjoint matrix C.
16 * It corresponds to the level 3 SYRK and level 2 SYR Blas routines.
17 **********************************************************************/
18 
19 namespace Eigen {
20 
21 
22 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
23 struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs>
24 {
25  static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha)
26  {
28  typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
30  for (Index i=0; i<size; ++i)
31  {
32  Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1)))
33  += (alpha * cj(vecY[i])) * ConjLhsType(OtherMap(vecX+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1)));
34  }
35  }
36 };
37 
38 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
39 struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs>
40 {
41  static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha)
42  {
44  }
45 };
46 
47 template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime>
49 
50 template<typename MatrixType, typename OtherType, int UpLo>
51 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
52 {
53  static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha)
54  {
55  typedef typename MatrixType::Scalar Scalar;
56  typedef internal::blas_traits<OtherType> OtherBlasTraits;
57  typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
58  typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
59  typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
60 
61  Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
62 
63  enum {
65  UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
66  };
68 
69  ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
70  (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
71 
72  if(!UseOtherDirectly)
73  Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
74 
75  selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
76  OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
77  (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
78  ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualOtherPtr, actualAlpha);
79  }
80 };
81 
82 template<typename MatrixType, typename OtherType, int UpLo>
83 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
84 {
85  static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha)
86  {
87  typedef typename MatrixType::Scalar Scalar;
88  typedef internal::blas_traits<OtherType> OtherBlasTraits;
89  typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
90  typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
91  typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
92 
93  Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
94 
95  enum {
96  IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0,
97  OtherIsRowMajor = _ActualOtherType::Flags&RowMajorBit ? 1 : 0
98  };
99 
100  Index size = mat.cols();
101  Index depth = actualOther.cols();
102 
104  MatrixType::MaxColsAtCompileTime, MatrixType::MaxColsAtCompileTime, _ActualOtherType::MaxColsAtCompileTime> BlockingType;
105 
106  BlockingType blocking(size, size, depth, 1, false);
107 
108 
110  Scalar, OtherIsRowMajor ? RowMajor : ColMajor, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
111  Scalar, OtherIsRowMajor ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
112  IsRowMajor ? RowMajor : ColMajor, MatrixType::InnerStrideAtCompileTime, UpLo>
113  ::run(size, depth,
114  actualOther.data(), actualOther.outerStride(), actualOther.data(), actualOther.outerStride(),
115  mat.data(), mat.innerStride(), mat.outerStride(), actualAlpha, blocking);
116  }
117 };
118 
119 // high level API
120 
121 template<typename MatrixType, unsigned int UpLo>
122 template<typename DerivedU>
125 {
126  selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
127 
128  return *this;
129 }
130 
131 } // end namespace Eigen
132 
133 #endif // EIGEN_SELFADJOINT_PRODUCT_H
SCALAR Scalar
Definition: bench_gemm.cpp:46
static void run(MatrixType &mat, const OtherType &other, const typename MatrixType::Scalar &alpha)
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
static double depth
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
const unsigned int RowMajorBit
Definition: Constants.h:66
internal::traits< SelfAdjointView >::Scalar Scalar
The type of coefficients in this matrix.
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
RealScalar alpha
static void run(MatrixType &mat, const OtherType &other, const typename MatrixType::Scalar &alpha)
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
static void run(Index size, Scalar *mat, Index stride, const Scalar *vecX, const Scalar *vecY, const Scalar &alpha)
EIGEN_DEVICE_FUNC SelfAdjointView & rankUpdate(const MatrixBase< DerivedU > &u, const MatrixBase< DerivedV > &v, const Scalar &alpha=Scalar(1))
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
Definition: Memory.h:768
#define EIGEN_DEVICE_FUNC
Definition: Macros.h:976
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
static void run(Index size, Scalar *mat, Index stride, const Scalar *vecX, const Scalar *vecY, const Scalar &alpha)
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:37