10 #ifndef EIGEN_SELFADJOINT_PRODUCT_H 11 #define EIGEN_SELFADJOINT_PRODUCT_H 22 template<
typename Scalar,
typename Index,
int UpLo,
bool ConjLhs,
bool ConjRhs>
23 struct selfadjoint_rank1_update<Scalar,
Index,
ColMajor,UpLo,ConjLhs,ConjRhs>
25 static void run(
Index size, Scalar* mat,
Index stride,
const Scalar* vecX,
const Scalar* vecY,
const Scalar& alpha)
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)
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)));
38 template<
typename Scalar,
typename Index,
int UpLo,
bool ConjLhs,
bool ConjRhs>
39 struct selfadjoint_rank1_update<Scalar,
Index,
RowMajor,UpLo,ConjLhs,ConjRhs>
41 static void run(
Index size, Scalar* mat,
Index stride,
const Scalar* vecX,
const Scalar* vecY,
const Scalar& alpha)
43 selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vecY,vecX,alpha);
47 template<
typename MatrixType,
typename OtherType,
int UpLo,
bool OtherIsVector = OtherType::IsVectorAtCompileTime>
48 struct selfadjoint_product_selector;
50 template<
typename MatrixType,
typename OtherType,
int UpLo>
51 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
53 static void run(MatrixType& mat,
const OtherType& other,
const typename MatrixType::Scalar& alpha)
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());
61 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
65 UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
67 internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
69 ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
70 (UseOtherDirectly ?
const_cast<Scalar*
>(actualOther.data()) : static_other.data()));
73 Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
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);
82 template<
typename MatrixType,
typename OtherType,
int UpLo>
83 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
85 static void run(MatrixType& mat,
const OtherType& other,
const typename MatrixType::Scalar& alpha)
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());
93 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
96 IsRowMajor = (internal::traits<MatrixType>::Flags&
RowMajorBit) ? 1 : 0,
97 OtherIsRowMajor = _ActualOtherType::Flags&
RowMajorBit ? 1 : 0
100 Index size = mat.cols();
101 Index depth = actualOther.cols();
103 typedef internal::gemm_blocking_space<IsRowMajor ?
RowMajor :
ColMajor,Scalar,Scalar,
104 MatrixType::MaxColsAtCompileTime, MatrixType::MaxColsAtCompileTime, _ActualOtherType::MaxColsAtCompileTime> BlockingType;
106 BlockingType blocking(size, size, depth, 1,
false);
109 internal::general_matrix_matrix_triangular_product<
Index,
110 Scalar, OtherIsRowMajor ?
RowMajor :
ColMajor, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
111 Scalar, OtherIsRowMajor ?
ColMajor :
RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
114 &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
115 mat.data(), mat.outerStride(), actualAlpha, blocking);
121 template<
typename MatrixType,
unsigned int UpLo>
122 template<
typename DerivedU>
124 ::rankUpdate(
const MatrixBase<DerivedU>& u,
const Scalar& alpha)
126 selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
133 #endif // EIGEN_SELFADJOINT_PRODUCT_H SelfAdjointView & rankUpdate(const MatrixBase< DerivedU > &u, const MatrixBase< DerivedV > &v, const Scalar &alpha=Scalar(1))
Definition: Constants.h:320
Namespace containing all symbols from the Eigen library.
Definition: Core:287
const unsigned int RowMajorBit
Definition: Constants.h:61
Definition: Constants.h:204
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Definition: Constants.h:322