mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-01 18:26:24 +08:00
Code Issues Packages Projects Releases Wiki Activity

this product still badly read the imaginary part on the diagonal

Browse Source
This commit is contained in:
Gael Guennebaud 2010-11-22 18:00:47 +01:00
parent a3f214ade9
commit 7213dd1e6b
1 changed files with 9 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
Eigen/src/Core/products/SelfadjointMatrixVector.h
View File

@ -40,6 +40,7 @@ static EIGEN_DONT_INLINE void product_selfadjoint_vector(
Scalar* res, Scalar alpha)
{
typedef typename packet_traits<Scalar>::type Packet;
typedef typename NumTraits<Scalar>::Real RealScalar;
const Index PacketSize = sizeof(Packet)/sizeof(Scalar);
enum {
@ -50,6 +51,7 @@ static EIGEN_DONT_INLINE void product_selfadjoint_vector(
conj_helper<Scalar,Scalar,NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, IsRowMajor), ConjugateRhs> cj0;
conj_helper<Scalar,Scalar,NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, !IsRowMajor), ConjugateRhs> cj1;
conj_helper<Scalar,Scalar,NumTraits<Scalar>::IsComplex, ConjugateRhs> cjd;
conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, IsRowMajor), ConjugateRhs> pcj0;
conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, !IsRowMajor), ConjugateRhs> pcj1;
@ -90,20 +92,20 @@ static EIGEN_DONT_INLINE void product_selfadjoint_vector(
size_t starti = FirstTriangular ? 0 : j+2;
size_t endi = FirstTriangular ? j : size;
size_t alignedEnd = starti;
size_t alignedStart = (starti) + first_aligned(&res[starti], endi-starti);
alignedEnd = alignedStart + ((endi-alignedStart)/(PacketSize))*(PacketSize);
size_t alignedEnd = alignedStart + ((endi-alignedStart)/(PacketSize))*(PacketSize);
res[j] += internal::real(A0[j]) * t0;
// TODO make sure this product is a real * complex and that the rhs is properly conjugated if needed
res[j] += cjd.pmul(internal::real(A0[j]), t0);
res[j+1] += cjd.pmul(internal::real(A1[j+1]), t1);
if(FirstTriangular)
{
res[j+1] += cj0.pmul(A1[j+1], t1);
res[j] += cj0.pmul(A1[j], t1);
t3 += cj1.pmul(A1[j], rhs[j]);
}
else
{
res[j+1] += cj0.pmul(A0[j+1],t0) + cj0.pmul(A1[j+1],t1);
res[j+1] += cj0.pmul(A0[j+1],t0);
t2 += cj1.pmul(A0[j+1], rhs[j+1]);
}
@ -147,7 +149,8 @@ static EIGEN_DONT_INLINE void product_selfadjoint_vector(
Scalar t1 = cjAlpha * rhs[j];
Scalar t2 = 0;
res[j] += internal::real(A0[j]) * t1;
// TODO make sure this product is a real * complex and that the rhs is properly conjugated if needed
res[j] += cjd.pmul(internal::real(A0[j]), t1);
for (Index i=FirstTriangular ? 0 : j+1; i<(FirstTriangular ? j : size); i++)
{
res[i] += cj0.pmul(A0[i], t1);