GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-13 00:21:49 +08:00
Code Issues Packages Projects Releases Wiki Activity

selfadjoint: do not reference the imaginary part of the diagonal

Browse Source
This commit is contained in:
Gael Guennebaud 2010-03-02 12:43:55 +01:00
parent abfed301cb
commit 7fd6458fec
1 changed files with 21 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
Eigen/src/Core/products/SelfadjointMatrixMatrix.h
View File

@ -43,7 +43,10 @@ struct ei_symm_pack_lhs
{ {
for(int w=0; w<h; w++) for(int w=0; w<h; w++)
blockA[count++] = ei_conj(lhs(k, i+w)); // transposed 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 blockA[count++] = lhs(i+w, k); // normal
++h; ++h;
} }
@ -71,8 +74,11 @@ struct ei_symm_pack_lhs
// do the same with mr==1 // do the same with mr==1
for(int i=peeled_mc; i<rows; i++) 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++] = lhs(i, k); // normal
blockA[count++] = ei_real(lhs(i, i)); // real (diagonal)
for(int k=i+1; k<cols; k++) for(int k=i+1; k<cols; k++)
blockA[count++] = ei_conj(lhs(k, i)); // transposed blockA[count++] = ei_conj(lhs(k, i)); // transposed
} }
@ -129,8 +135,11 @@ struct ei_symm_pack_rhs
// normal // normal
for (int w=0 ; w<h; ++w) for (int w=0 ; w<h; ++w)
ei_pstore(&blockB[count+w*PacketSize], ei_pset1(alpha*rhs(k,j2+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 // 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)))); ei_pstore(&blockB[count+w*PacketSize], ei_pset1(alpha*ei_conj(rhs(j2+w,k))));
count += nr*PacketSize; count += nr*PacketSize;
++h; ++h;
@ -175,8 +184,15 @@ struct ei_symm_pack_rhs
ei_pstore(&blockB[count], ei_pset1(alpha*ei_conj(rhs(j2,k)))); ei_pstore(&blockB[count], ei_pset1(alpha*ei_conj(rhs(j2,k))));
count += PacketSize; count += PacketSize;
} }
if(half==j2)
{
ei_pstore(&blockB[count], ei_pset1(alpha*ei_real(rhs(j2,j2))));
count += PacketSize;
}
// normal // 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))); ei_pstore(&blockB[count], ei_pset1(alpha*rhs(k,j2)));
count += PacketSize; count += PacketSize;
@ -385,7 +401,7 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>
* RhsBlasTraits::extractScalarFactor(m_rhs); * RhsBlasTraits::extractScalarFactor(m_rhs);
ei_product_selfadjoint_matrix<Scalar, ei_product_selfadjoint_matrix<Scalar,
EIGEN_LOGICAL_XOR(LhsIsUpper, EIGEN_LOGICAL_XOR(LhsIsUpper,
ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint, ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint,
NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsIsUpper,bool(LhsBlasTraits::NeedToConjugate)), NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsIsUpper,bool(LhsBlasTraits::NeedToConjugate)),
EIGEN_LOGICAL_XOR(RhsIsUpper, EIGEN_LOGICAL_XOR(RhsIsUpper,