SelfadjointRank2Update.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_SELFADJOINTRANK2UPTADE_H
11 #define EIGEN_SELFADJOINTRANK2UPTADE_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 /* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu'
18  * It corresponds to the Level2 syr2 BLAS routine
19  */
20 
21 template<typename Scalar, typename Index, typename UType, typename VType, int UpLo>
23 
24 template<typename Scalar, typename Index, typename UType, typename VType>
26 {
27  static EIGEN_DEVICE_FUNC
28  void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha)
29  {
30  const Index size = u.size();
31  for (Index i=0; i<size; ++i)
32  {
34  (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.tail(size-i)
35  + (alpha * numext::conj(v.coeff(i))) * u.tail(size-i);
36  }
37  }
38 };
39 
40 template<typename Scalar, typename Index, typename UType, typename VType>
42 {
43  static void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha)
44  {
45  const Index size = u.size();
46  for (Index i=0; i<size; ++i)
47  Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) +=
48  (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.head(i+1)
49  + (alpha * numext::conj(v.coeff(i))) * u.head(i+1);
50  }
51 };
52 
53 template<bool Cond, typename T> struct conj_expr_if
54  : conditional<!Cond, const T&,
55  CwiseUnaryOp<scalar_conjugate_op<typename traits<T>::Scalar>,T> > {};
56 
57 } // end namespace internal
58 
59 template<typename MatrixType, unsigned int UpLo>
60 template<typename DerivedU, typename DerivedV>
63 {
64  typedef internal::blas_traits<DerivedU> UBlasTraits;
65  typedef typename UBlasTraits::DirectLinearAccessType ActualUType;
66  typedef typename internal::remove_all<ActualUType>::type _ActualUType;
67  typename internal::add_const_on_value_type<ActualUType>::type actualU = UBlasTraits::extract(u.derived());
68 
69  typedef internal::blas_traits<DerivedV> VBlasTraits;
70  typedef typename VBlasTraits::DirectLinearAccessType ActualVType;
71  typedef typename internal::remove_all<ActualVType>::type _ActualVType;
72  typename internal::add_const_on_value_type<ActualVType>::type actualV = VBlasTraits::extract(v.derived());
73 
74  // If MatrixType is row major, then we use the routine for lower triangular in the upper triangular case and
75  // vice versa, and take the complex conjugate of all coefficients and vector entries.
76 
77  enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
78  Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived())
79  * numext::conj(VBlasTraits::extractScalarFactor(v.derived()));
80  if (IsRowMajor)
81  actualAlpha = numext::conj(actualAlpha);
82 
83  typedef typename internal::remove_all<typename internal::conj_expr_if<int(IsRowMajor) ^ int(UBlasTraits::NeedToConjugate), _ActualUType>::type>::type UType;
84  typedef typename internal::remove_all<typename internal::conj_expr_if<int(IsRowMajor) ^ int(VBlasTraits::NeedToConjugate), _ActualVType>::type>::type VType;
86  (IsRowMajor ? int(UpLo==Upper ? Lower : Upper) : UpLo)>
87  ::run(_expression().const_cast_derived().data(),_expression().outerStride(),UType(actualU),VType(actualV),actualAlpha);
88 
89  return *this;
90 }
91 
92 } // end namespace Eigen
93 
94 #endif // EIGEN_SELFADJOINTRANK2UPTADE_H
gtsam.examples.DogLegOptimizerExample.int
int
Definition: DogLegOptimizerExample.py:111
EIGEN_DEVICE_FUNC
#define EIGEN_DEVICE_FUNC
Definition: Macros.h:976
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
Eigen::internal::remove_all
Definition: Meta.h:126
Eigen::internal::selfadjoint_rank2_update_selector< Scalar, Index, UType, VType, Lower >::run
static EIGEN_DEVICE_FUNC void run(Scalar *mat, Index stride, const UType &u, const VType &v, const Scalar &alpha)
Definition: SelfadjointRank2Update.h:28
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:66
Eigen::Upper
@ Upper
Definition: Constants.h:211
type
Definition: pytypes.h:1525
Eigen::SelfAdjointView
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
Definition: SelfAdjointView.h:49
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
Eigen::SelfAdjointView< const Derived, UpLo >::Scalar
internal::traits< SelfAdjointView >::Scalar Scalar
The type of coefficients in this matrix.
Definition: SelfAdjointView.h:61
Eigen::internal::selfadjoint_rank2_update_selector< Scalar, Index, UType, VType, Upper >::run
static void run(Scalar *mat, Index stride, const UType &u, const VType &v, const Scalar &alpha)
Definition: SelfadjointRank2Update.h:43
gtsam.examples.DogLegOptimizerExample.run
def run(args)
Definition: DogLegOptimizerExample.py:21
Eigen::Triplet< double >
conj
AnnoyingScalar conj(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:104
Eigen::Lower
@ Lower
Definition: Constants.h:209
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
Eigen::internal::conditional
Definition: Meta.h:109
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
Eigen::internal::blas_traits
Definition: BlasUtil.h:402
Eigen::Matrix< Scalar, Dynamic, 1 >
Eigen::internal::selfadjoint_rank2_update_selector
Definition: SelfadjointRank2Update.h:22
internal
Definition: BandTriangularSolver.h:13
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Eigen::internal::size
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
Eigen::internal::conj_expr_if
Definition: SelfadjointRank2Update.h:53
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::SelfAdjointView::rankUpdate
EIGEN_DEVICE_FUNC SelfAdjointView & rankUpdate(const MatrixBase< DerivedU > &u, const MatrixBase< DerivedV > &v, const Scalar &alpha=Scalar(1))
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Eigen::internal::add_const_on_value_type
Definition: Meta.h:214


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:35:05