10 #ifndef EIGEN_SELFADJOINT_PRODUCT_H
11 #define EIGEN_SELFADJOINT_PRODUCT_H
21 template<
typename Scalar,
typename Index,
int StorageOrder,
int UpLo,
bool ConjLhs,
bool ConjRhs>
22 struct selfadjoint_rank1_update;
24 template<
typename Scalar,
typename Index,
int UpLo,
bool ConjLhs,
bool ConjRhs>
25 struct selfadjoint_rank1_update<Scalar,Index,
ColMajor,UpLo,ConjLhs,ConjRhs>
27 static void run(Index size, Scalar* mat, Index stride,
const Scalar* vec, Scalar alpha)
29 internal::conj_if<ConjRhs> cj;
30 typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
31 typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjRhsType;
32 for (Index i=0; i<size; ++i)
34 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==
Lower ? i : 0), (UpLo==
Lower ? size-i : (i+1)))
35 += (alpha * cj(vec[i])) * ConjRhsType(OtherMap(vec+(UpLo==
Lower ? i : 0),UpLo==
Lower ? size-i : (i+1)));
40 template<
typename Scalar,
typename Index,
int UpLo,
bool ConjLhs,
bool ConjRhs>
41 struct selfadjoint_rank1_update<Scalar,Index,
RowMajor,UpLo,ConjLhs,ConjRhs>
43 static void run(Index size, Scalar* mat, Index stride,
const Scalar* vec, Scalar alpha)
45 selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vec,alpha);
49 template<
typename MatrixType,
typename OtherType,
int UpLo,
bool OtherIsVector = OtherType::IsVectorAtCompileTime>
50 struct selfadjoint_product_selector;
52 template<
typename MatrixType,
typename OtherType,
int UpLo>
53 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
55 static void run(MatrixType& mat,
const OtherType& other,
typename MatrixType::Scalar alpha)
57 typedef typename MatrixType::Scalar Scalar;
58 typedef typename MatrixType::Index Index;
59 typedef internal::blas_traits<OtherType> OtherBlasTraits;
60 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
61 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
62 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
64 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
68 UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
70 internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
72 ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
73 (UseOtherDirectly ?
const_cast<Scalar*
>(actualOther.data()) : static_other.data()));
76 Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
78 selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
79 OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
80 (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
81 ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha);
85 template<
typename MatrixType,
typename OtherType,
int UpLo>
86 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
88 static void run(MatrixType& mat,
const OtherType& other,
typename MatrixType::Scalar alpha)
90 typedef typename MatrixType::Scalar Scalar;
91 typedef typename MatrixType::Index Index;
92 typedef internal::blas_traits<OtherType> OtherBlasTraits;
93 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
94 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
95 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
97 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
99 enum { IsRowMajor = (internal::traits<MatrixType>::Flags&
RowMajorBit) ? 1 : 0 };
101 internal::general_matrix_matrix_triangular_product<Index,
102 Scalar, _ActualOtherType::Flags&
RowMajorBit ?
RowMajor :
ColMajor, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
103 Scalar, _ActualOtherType::Flags&
RowMajorBit ? ColMajor :
RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
105 ::run(mat.cols(), actualOther.cols(),
106 &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
107 mat.data(), mat.outerStride(), actualAlpha);
113 template<
typename MatrixType,
unsigned int UpLo>
114 template<
typename DerivedU>
118 selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
125 #endif // EIGEN_SELFADJOINT_PRODUCT_H