diff options
| author |  Gael Guennebaud <g.gael@free.fr> | 2010-03-02 12:43:55 +0100 |
|---|---|---|
| committer | Gael Guennebaud <g.gael@free.fr> | 2010-03-02 12:43:55 +0100 |
| commit | 7fd6458fec694f213323d6dd0718d315513adbb5 (patch) | |
| tree | 2e379158a5f510df5ebdb956869ed58744298ff7 /Eigen/src/Core/products/SelfadjointMatrixMatrix.h | |
| parent | abfed301cb474c27fbb76a41cc459602db2b145f (diff) | |
selfadjoint: do not reference the imaginary part of the diagonal
Diffstat (limited to 'Eigen/src/Core/products/SelfadjointMatrixMatrix.h')
| -rw-r--r-- | Eigen/src/Core/products/SelfadjointMatrixMatrix.h | 26 |
1 files changed, 21 insertions, 5 deletions
diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h index 89cbc3ac0..84d056c5d 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h +++ b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h @@ -43,7 +43,10 @@ struct ei_symm_pack_lhs { for(int w=0; w<h; w++) blockA[count++] = ei_conj(lhs(k, i+w)); // transposed - for(int w=h; w<BlockRows; w++) + + blockA[count++] = ei_real(lhs(k,k)); // real (diagonal) + + for(int w=h+1; w<BlockRows; w++) blockA[count++] = lhs(i+w, k); // normal ++h; } @@ -71,8 +74,11 @@ struct ei_symm_pack_lhs // do the same with mr==1 for(int i=peeled_mc; i<rows; i++) { - for(int k=0; k<=i; k++) + for(int k=0; k<i; k++) blockA[count++] = lhs(i, k); // normal + + blockA[count++] = ei_real(lhs(i, i)); // real (diagonal) + for(int k=i+1; k<cols; k++) blockA[count++] = ei_conj(lhs(k, i)); // transposed } @@ -129,8 +135,11 @@ struct ei_symm_pack_rhs // normal for (int w=0 ; w<h; ++w) ei_pstore(&blockB[count+w*PacketSize], ei_pset1(alpha*rhs(k,j2+w))); + + ei_pstore(&blockB[count+h*PacketSize], ei_pset1(alpha*ei_real(rhs(k,k)))); + // transpose - for (int w=h ; w<nr; ++w) + for (int w=h+1 ; w<nr; ++w) ei_pstore(&blockB[count+w*PacketSize], ei_pset1(alpha*ei_conj(rhs(j2+w,k)))); count += nr*PacketSize; ++h; @@ -175,8 +184,15 @@ struct ei_symm_pack_rhs ei_pstore(&blockB[count], ei_pset1(alpha*ei_conj(rhs(j2,k)))); count += PacketSize; } + + if(half==j2) + { + ei_pstore(&blockB[count], ei_pset1(alpha*ei_real(rhs(j2,j2)))); + count += PacketSize; + } + // normal - for(int k=half; k<k2+rows; k++) + for(int k=half+1; k<k2+rows; k++) { ei_pstore(&blockB[count], ei_pset1(alpha*rhs(k,j2))); count += PacketSize; @@ -385,7 +401,7 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> * RhsBlasTraits::extractScalarFactor(m_rhs); ei_product_selfadjoint_matrix<Scalar, - EIGEN_LOGICAL_XOR(LhsIsUpper, + EIGEN_LOGICAL_XOR(LhsIsUpper, ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint, NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsIsUpper,bool(LhsBlasTraits::NeedToConjugate)), EIGEN_LOGICAL_XOR(RhsIsUpper, |