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_SELFADJOINTRANK2UPTADE_H
00026 #define EIGEN_SELFADJOINTRANK2UPTADE_H
00027
00028 namespace Eigen {
00029
00030 namespace internal {
00031
00032
00033
00034
00035
00036 template<typename Scalar, typename Index, typename UType, typename VType, int UpLo>
00037 struct selfadjoint_rank2_update_selector;
00038
00039 template<typename Scalar, typename Index, typename UType, typename VType>
00040 struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Lower>
00041 {
00042 static void run(Scalar* mat, Index stride, const UType& u, const VType& v, Scalar alpha)
00043 {
00044 const Index size = u.size();
00045 for (Index i=0; i<size; ++i)
00046 {
00047 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) +=
00048 (conj(alpha) * conj(u.coeff(i))) * v.tail(size-i)
00049 + (alpha * conj(v.coeff(i))) * u.tail(size-i);
00050 }
00051 }
00052 };
00053
00054 template<typename Scalar, typename Index, typename UType, typename VType>
00055 struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Upper>
00056 {
00057 static void run(Scalar* mat, Index stride, const UType& u, const VType& v, Scalar alpha)
00058 {
00059 const Index size = u.size();
00060 for (Index i=0; i<size; ++i)
00061 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) +=
00062 (conj(alpha) * conj(u.coeff(i))) * v.head(i+1)
00063 + (alpha * conj(v.coeff(i))) * u.head(i+1);
00064 }
00065 };
00066
00067 template<bool Cond, typename T> struct conj_expr_if
00068 : conditional<!Cond, const T&,
00069 CwiseUnaryOp<scalar_conjugate_op<typename traits<T>::Scalar>,T> > {};
00070
00071 }
00072
00073 template<typename MatrixType, unsigned int UpLo>
00074 template<typename DerivedU, typename DerivedV>
00075 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
00076 ::rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha)
00077 {
00078 typedef internal::blas_traits<DerivedU> UBlasTraits;
00079 typedef typename UBlasTraits::DirectLinearAccessType ActualUType;
00080 typedef typename internal::remove_all<ActualUType>::type _ActualUType;
00081 typename internal::add_const_on_value_type<ActualUType>::type actualU = UBlasTraits::extract(u.derived());
00082
00083 typedef internal::blas_traits<DerivedV> VBlasTraits;
00084 typedef typename VBlasTraits::DirectLinearAccessType ActualVType;
00085 typedef typename internal::remove_all<ActualVType>::type _ActualVType;
00086 typename internal::add_const_on_value_type<ActualVType>::type actualV = VBlasTraits::extract(v.derived());
00087
00088
00089
00090
00091 enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
00092 Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived())
00093 * internal::conj(VBlasTraits::extractScalarFactor(v.derived()));
00094 if (IsRowMajor)
00095 actualAlpha = internal::conj(actualAlpha);
00096
00097 internal::selfadjoint_rank2_update_selector<Scalar, Index,
00098 typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::type>::type,
00099 typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::type>::type,
00100 (IsRowMajor ? int(UpLo==Upper ? Lower : Upper) : UpLo)>
00101 ::run(_expression().const_cast_derived().data(),_expression().outerStride(),actualU,actualV,actualAlpha);
00102
00103 return *this;
00104 }
00105
00106 }
00107
00108 #endif // EIGEN_SELFADJOINTRANK2UPTADE_H