00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_SELFADJOINT_PRODUCT_H
00026 #define EIGEN_SELFADJOINT_PRODUCT_H
00027
00028
00029
00030
00031
00032
00033
00034 namespace Eigen {
00035
00036 template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs>
00037 struct selfadjoint_rank1_update;
00038
00039 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
00040 struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs>
00041 {
00042 static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha)
00043 {
00044 internal::conj_if<ConjRhs> cj;
00045 typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
00046 typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjRhsType;
00047 for (Index i=0; i<size; ++i)
00048 {
00049 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1)))
00050 += (alpha * cj(vec[i])) * ConjRhsType(OtherMap(vec+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1)));
00051 }
00052 }
00053 };
00054
00055 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
00056 struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs>
00057 {
00058 static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha)
00059 {
00060 selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vec,alpha);
00061 }
00062 };
00063
00064 template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime>
00065 struct selfadjoint_product_selector;
00066
00067 template<typename MatrixType, typename OtherType, int UpLo>
00068 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
00069 {
00070 static void run(MatrixType& mat, const OtherType& other, typename MatrixType::Scalar alpha)
00071 {
00072 typedef typename MatrixType::Scalar Scalar;
00073 typedef typename MatrixType::Index Index;
00074 typedef internal::blas_traits<OtherType> OtherBlasTraits;
00075 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
00076 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
00077 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
00078
00079 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
00080
00081 enum {
00082 StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor,
00083 UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
00084 };
00085 internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
00086
00087 ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
00088 (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
00089
00090 if(!UseOtherDirectly)
00091 Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
00092
00093 selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
00094 OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
00095 (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
00096 ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha);
00097 }
00098 };
00099
00100 template<typename MatrixType, typename OtherType, int UpLo>
00101 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
00102 {
00103 static void run(MatrixType& mat, const OtherType& other, typename MatrixType::Scalar alpha)
00104 {
00105 typedef typename MatrixType::Scalar Scalar;
00106 typedef typename MatrixType::Index Index;
00107 typedef internal::blas_traits<OtherType> OtherBlasTraits;
00108 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
00109 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
00110 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
00111
00112 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
00113
00114 enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
00115
00116 internal::general_matrix_matrix_triangular_product<Index,
00117 Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
00118 Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
00119 MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo>
00120 ::run(mat.cols(), actualOther.cols(),
00121 &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
00122 mat.data(), mat.outerStride(), actualAlpha);
00123 }
00124 };
00125
00126
00127
00128 template<typename MatrixType, unsigned int UpLo>
00129 template<typename DerivedU>
00130 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
00131 ::rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha)
00132 {
00133 selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
00134
00135 return *this;
00136 }
00137
00138 }
00139
00140 #endif // EIGEN_SELFADJOINT_PRODUCT_H