SelfadjointProduct.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba 
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
19namespace Eigen {
20
21
22template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
23struct 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  {
27    internal::conj_if<ConjRhs> cj;
28    typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
29    typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjLhsType;
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
38template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
39struct 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  {
43    selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vecY,vecX,alpha);
44  }
45};
46
47template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime>
48struct selfadjoint_product_selector;
49
50template<typename MatrixType, typename OtherType, int UpLo>
51struct 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 typename MatrixType::Index Index;
57    typedef internal::blas_traits<OtherType> OtherBlasTraits;
58    typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
59    typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
60    typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
61
62    Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
63
64    enum {
65      StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor,
66      UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
67    };
68    internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
69
70    ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
71      (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
72
73    if(!UseOtherDirectly)
74      Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
75
76    selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
77                              OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
78                            (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
79          ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualOtherPtr, actualAlpha);
80  }
81};
82
83template<typename MatrixType, typename OtherType, int UpLo>
84struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
85{
86  static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha)
87  {
88    typedef typename MatrixType::Scalar Scalar;
89    typedef typename MatrixType::Index Index;
90    typedef internal::blas_traits<OtherType> OtherBlasTraits;
91    typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
92    typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
93    typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
94
95    Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
96
97    enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
98
99    internal::general_matrix_matrix_triangular_product<Index,
100      Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor,   OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
101      Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
102      MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo>
103      ::run(mat.cols(), actualOther.cols(),
104            &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
105            mat.data(), mat.outerStride(), actualAlpha);
106  }
107};
108
109// high level API
110
111template<typename MatrixType, unsigned int UpLo>
112template<typename DerivedU>
113SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
114::rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha)
115{
116  selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
117
118  return *this;
119}
120
121} // end namespace Eigen
122
123#endif // EIGEN_SELFADJOINT_PRODUCT_H
124