diff options
author | Gael Guennebaud <g.gael@free.fr> | 2012-07-24 18:19:56 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2012-07-24 18:19:56 +0200 |
commit | 7b34b5f6f98115efbda9d127378aea5d403a4905 (patch) | |
tree | 2f19605b569a48b5e233723201e2f5eb148fe242 /Eigen/src/Jacobi | |
parent | e7c07de549ace45f11c000c15f6a9fcb1c9b8170 (diff) |
do not apply plane rotation when it is exactly the identity
Diffstat (limited to 'Eigen/src/Jacobi')
-rw-r--r-- | Eigen/src/Jacobi/Jacobi.h | 26 |
1 files changed, 15 insertions, 11 deletions
diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index a9c17dcdf..7eb19a023 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -207,7 +207,6 @@ void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar template<typename Scalar> void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::false_type) { - if(q==Scalar(0)) { m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1); @@ -303,6 +302,11 @@ void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0); Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0); + + OtherScalar c = j.c(); + OtherScalar s = j.s(); + if (c==OtherScalar(1) && s==OtherScalar(0)) + return; /*** dynamic-size vectorized paths ***/ @@ -316,16 +320,16 @@ void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, Index alignedStart = internal::first_aligned(y, size); Index alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize; - const Packet pc = pset1<Packet>(j.c()); - const Packet ps = pset1<Packet>(j.s()); + const Packet pc = pset1<Packet>(c); + const Packet ps = pset1<Packet>(s); conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex,false> pcj; for(Index i=0; i<alignedStart; ++i) { Scalar xi = x[i]; Scalar yi = y[i]; - x[i] = j.c() * xi + conj(j.s()) * yi; - y[i] = -j.s() * xi + conj(j.c()) * yi; + x[i] = c * xi + conj(s) * yi; + y[i] = -s * xi + conj(c) * yi; } Scalar* EIGEN_RESTRICT px = x + alignedStart; @@ -372,8 +376,8 @@ void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, { Scalar xi = x[i]; Scalar yi = y[i]; - x[i] = j.c() * xi + conj(j.s()) * yi; - y[i] = -j.s() * xi + conj(j.c()) * yi; + x[i] = c * xi + conj(s) * yi; + y[i] = -s * xi + conj(c) * yi; } } @@ -382,8 +386,8 @@ void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, (VectorX::Flags & VectorY::Flags & PacketAccessBit) && (VectorX::Flags & VectorY::Flags & AlignedBit)) { - const Packet pc = pset1<Packet>(j.c()); - const Packet ps = pset1<Packet>(j.s()); + const Packet pc = pset1<Packet>(c); + const Packet ps = pset1<Packet>(s); conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex,false> pcj; Scalar* EIGEN_RESTRICT px = x; Scalar* EIGEN_RESTRICT py = y; @@ -405,8 +409,8 @@ void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, { Scalar xi = *x; Scalar yi = *y; - *x = j.c() * xi + conj(j.s()) * yi; - *y = -j.s() * xi + conj(j.c()) * yi; + *x = c * xi + conj(s) * yi; + *y = -s * xi + conj(c) * yi; x += incrx; y += incry; } |