diff options
| author |  Gael Guennebaud <g.gael@free.fr> | 2013-02-24 23:05:42 +0100 |
|---|---|---|
| committer | Gael Guennebaud <g.gael@free.fr> | 2013-02-24 23:05:42 +0100 |
| commit | 04367447ac295d9818713f54b0a539efef7f0caa (patch) | |
| tree | c9d2d53d94ecaef1f6c8a3ea923bb6a605a58104 /Eigen/src/Core/products/SelfadjointProduct.h | |
| parent | 08388cc71209ac29ad768c628f1b014dd460d819 (diff) | |
Fix bug #496: generalize internal rank1_update implementation to accept uplo(A) += v * w and make A.triangularView() += v * w uses it.
Update unit tests and blas interface respectively.
Diffstat (limited to 'Eigen/src/Core/products/SelfadjointProduct.h')
| -rw-r--r-- | Eigen/src/Core/products/SelfadjointProduct.h | 14 |
1 files changed, 6 insertions, 8 deletions
diff --git a/Eigen/src/Core/products/SelfadjointProduct.h b/Eigen/src/Core/products/SelfadjointProduct.h index 6a55f3d77..302e0d841 100644 --- a/Eigen/src/Core/products/SelfadjointProduct.h +++ b/Eigen/src/Core/products/SelfadjointProduct.h @@ -18,21 +18,19 @@ namespace Eigen { -template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs> -struct selfadjoint_rank1_update; template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs> { - static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha) + static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, Scalar alpha) { internal::conj_if<ConjRhs> cj; typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; - typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjRhsType; + typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjLhsType; for (Index i=0; i<size; ++i) { Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1))) - += (alpha * cj(vec[i])) * ConjRhsType(OtherMap(vec+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1))); + += (alpha * cj(vecY[i])) * ConjLhsType(OtherMap(vecX+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1))); } } }; @@ -40,9 +38,9 @@ struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs> template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs> { - static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha) + static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, Scalar alpha) { - selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vec,alpha); + selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vecY,vecX,alpha); } }; @@ -78,7 +76,7 @@ struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true> selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex> - ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha); + ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualOtherPtr, actualAlpha); } }; |