25 #ifndef EIGEN_SELFADJOINT_PRODUCT_H
26 #define EIGEN_SELFADJOINT_PRODUCT_H
36 template<
typename Scalar,
typename Index,
int StorageOrder,
int UpLo,
bool ConjLhs,
bool ConjRhs>
37 struct selfadjoint_rank1_update;
39 template<
typename Scalar,
typename Index,
int UpLo,
bool ConjLhs,
bool ConjRhs>
40 struct selfadjoint_rank1_update<Scalar,Index,
ColMajor,UpLo,ConjLhs,ConjRhs>
42 static void run(Index size, Scalar* mat, Index stride,
const Scalar* vec, Scalar alpha)
44 internal::conj_if<ConjRhs> cj;
46 typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjRhsType;
47 for (Index i=0; i<size; ++i)
50 += (alpha * cj(vec[i])) * ConjRhsType(OtherMap(vec+(UpLo==
Lower ? i : 0),UpLo==
Lower ? size-i : (i+1)));
55 template<
typename Scalar,
typename Index,
int UpLo,
bool ConjLhs,
bool ConjRhs>
56 struct selfadjoint_rank1_update<Scalar,Index,
RowMajor,UpLo,ConjLhs,ConjRhs>
58 static void run(Index size, Scalar* mat, Index stride,
const Scalar* vec, Scalar alpha)
60 selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vec,alpha);
64 template<
typename MatrixType,
typename OtherType,
int UpLo,
bool OtherIsVector = OtherType::IsVectorAtCompileTime>
65 struct selfadjoint_product_selector;
67 template<
typename MatrixType,
typename OtherType,
int UpLo>
68 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
70 static void run(MatrixType& mat,
const OtherType& other,
typename MatrixType::Scalar alpha)
72 typedef typename MatrixType::Scalar Scalar;
73 typedef typename MatrixType::Index Index;
74 typedef internal::blas_traits<OtherType> OtherBlasTraits;
75 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
76 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
77 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
79 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
83 UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
85 internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
88 (UseOtherDirectly ?
const_cast<Scalar*
>(actualOther.data()) : static_other.data()));
93 selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
95 (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
96 ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha);
100 template<
typename MatrixType,
typename OtherType,
int UpLo>
101 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
103 static void run(MatrixType& mat,
const OtherType& other,
typename MatrixType::Scalar alpha)
105 typedef typename MatrixType::Scalar Scalar;
106 typedef typename MatrixType::Index Index;
107 typedef internal::blas_traits<OtherType> OtherBlasTraits;
108 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
109 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
110 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
112 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
114 enum { IsRowMajor = (internal::traits<MatrixType>::Flags&
RowMajorBit) ? 1 : 0 };
116 internal::general_matrix_matrix_triangular_product<Index,
118 Scalar, _ActualOtherType::Flags&
RowMajorBit ?
ColMajor :
RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
120 ::run(mat.cols(), actualOther.cols(),
121 &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
122 mat.data(), mat.outerStride(), actualAlpha);
128 template<
typename MatrixType,
unsigned int UpLo>
129 template<
typename DerivedU>
133 selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
140 #endif // EIGEN_SELFADJOINT_PRODUCT_H