diff options
Diffstat (limited to 'blas/PackedSelfadjointProduct.h')
| -rw-r--r-- | blas/PackedSelfadjointProduct.h | 19 |
1 files changed, 7 insertions, 12 deletions
diff --git a/blas/PackedSelfadjointProduct.h b/blas/PackedSelfadjointProduct.h index 1ba67b9c1..f7c9b9341 100644 --- a/blas/PackedSelfadjointProduct.h +++ b/blas/PackedSelfadjointProduct.h @@ -14,12 +14,6 @@ namespace internal { /* Optimized matrix += alpha * uv' * The matrix is in packed form. - * - * FIXME I always fail tests for complex self-adjoint matrices. - * - * ******* FATAL ERROR - PARAMETER NUMBER 6 WAS CHANGED INCORRECTLY ******* - * ******* xHPR FAILED ON CALL NUMBER: - * 2: xHPR ('U', 1, 0.0, X, 1, AP) */ template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs> struct selfadjoint_packed_rank1_update; @@ -27,20 +21,20 @@ struct selfadjoint_packed_rank1_update; template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> struct selfadjoint_packed_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs> { - static void run(Index size, Scalar* mat, const Scalar* vec, Scalar alpha) + typedef typename NumTraits<Scalar>::Real RealScalar; + static void run(Index size, Scalar* mat, const Scalar* vec, RealScalar alpha) { typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; typedef typename conj_expr_if<ConjLhs,OtherMap>::type ConjRhsType; conj_if<ConjRhs> cj; - Index offset = 0; for (Index i=0; i<size; ++i) { - Map<Matrix<Scalar,Dynamic,1> >(mat+offset, UpLo==Lower ? size-i : (i+1)) + Map<Matrix<Scalar,Dynamic,1> >(mat, UpLo==Lower ? size-i : (i+1)) += alpha * cj(vec[i]) * ConjRhsType(OtherMap(vec+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1))); //FIXME This should be handled outside. - mat[offset+(UpLo==Lower ? 0 : i)] = real(mat[offset+(UpLo==Lower ? 0 : i)]); - offset += UpLo==Lower ? size-i : (i+1); + mat[UpLo==Lower ? 0 : i] = real(mat[UpLo==Lower ? 0 : i]); + mat += UpLo==Lower ? size-i : (i+1); } } }; @@ -48,7 +42,8 @@ struct selfadjoint_packed_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRh template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> struct selfadjoint_packed_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs> { - static void run(Index size, Scalar* mat, const Scalar* vec, Scalar alpha) + typedef typename NumTraits<Scalar>::Real RealScalar; + static void run(Index size, Scalar* mat, const Scalar* vec, RealScalar alpha) { selfadjoint_packed_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,vec,alpha); } |