diff options
Diffstat (limited to 'extern/Eigen3/Eigen/src/Jacobi')
-rw-r--r-- | extern/Eigen3/Eigen/src/Jacobi/Jacobi.h | 99 |
1 files changed, 56 insertions, 43 deletions
diff --git a/extern/Eigen3/Eigen/src/Jacobi/Jacobi.h b/extern/Eigen3/Eigen/src/Jacobi/Jacobi.h index a9c17dcdf19..956f72d5701 100644 --- a/extern/Eigen3/Eigen/src/Jacobi/Jacobi.h +++ b/extern/Eigen3/Eigen/src/Jacobi/Jacobi.h @@ -50,19 +50,20 @@ template<typename Scalar> class JacobiRotation /** Concatenates two planar rotation */ JacobiRotation operator*(const JacobiRotation& other) { - return JacobiRotation(m_c * other.m_c - internal::conj(m_s) * other.m_s, - internal::conj(m_c * internal::conj(other.m_s) + internal::conj(m_s) * internal::conj(other.m_c))); + using numext::conj; + return JacobiRotation(m_c * other.m_c - conj(m_s) * other.m_s, + conj(m_c * conj(other.m_s) + conj(m_s) * conj(other.m_c))); } /** Returns the transposed transformation */ - JacobiRotation transpose() const { return JacobiRotation(m_c, -internal::conj(m_s)); } + JacobiRotation transpose() const { using numext::conj; return JacobiRotation(m_c, -conj(m_s)); } /** Returns the adjoint transformation */ - JacobiRotation adjoint() const { return JacobiRotation(internal::conj(m_c), -m_s); } + JacobiRotation adjoint() const { using numext::conj; return JacobiRotation(conj(m_c), -m_s); } template<typename Derived> bool makeJacobi(const MatrixBase<Derived>&, typename Derived::Index p, typename Derived::Index q); - bool makeJacobi(RealScalar x, Scalar y, RealScalar z); + bool makeJacobi(const RealScalar& x, const Scalar& y, const RealScalar& z); void makeGivens(const Scalar& p, const Scalar& q, Scalar* z=0); @@ -79,8 +80,10 @@ template<typename Scalar> class JacobiRotation * \sa MatrixBase::makeJacobi(const MatrixBase<Derived>&, Index, Index), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ template<typename Scalar> -bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z) +bool JacobiRotation<Scalar>::makeJacobi(const RealScalar& x, const Scalar& y, const RealScalar& z) { + using std::sqrt; + using std::abs; typedef typename NumTraits<Scalar>::Real RealScalar; if(y == Scalar(0)) { @@ -90,8 +93,8 @@ bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z) } else { - RealScalar tau = (x-z)/(RealScalar(2)*internal::abs(y)); - RealScalar w = internal::sqrt(internal::abs2(tau) + RealScalar(1)); + RealScalar tau = (x-z)/(RealScalar(2)*abs(y)); + RealScalar w = sqrt(numext::abs2(tau) + RealScalar(1)); RealScalar t; if(tau>RealScalar(0)) { @@ -102,8 +105,8 @@ bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z) t = RealScalar(1) / (tau - w); } RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1); - RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+RealScalar(1)); - m_s = - sign_t * (internal::conj(y) / internal::abs(y)) * internal::abs(t) * n; + RealScalar n = RealScalar(1) / sqrt(numext::abs2(t)+RealScalar(1)); + m_s = - sign_t * (numext::conj(y) / abs(y)) * abs(t) * n; m_c = n; return true; } @@ -122,7 +125,7 @@ template<typename Scalar> template<typename Derived> inline bool JacobiRotation<Scalar>::makeJacobi(const MatrixBase<Derived>& m, typename Derived::Index p, typename Derived::Index q) { - return makeJacobi(internal::real(m.coeff(p,p)), m.coeff(p,q), internal::real(m.coeff(q,q))); + return makeJacobi(numext::real(m.coeff(p,p)), m.coeff(p,q), numext::real(m.coeff(q,q))); } /** Makes \c *this as a Givens rotation \c G such that applying \f$ G^* \f$ to the left of the vector @@ -152,52 +155,56 @@ 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::true_type) { + using std::sqrt; + using std::abs; + using numext::conj; + if(q==Scalar(0)) { - m_c = internal::real(p)<0 ? Scalar(-1) : Scalar(1); + m_c = numext::real(p)<0 ? Scalar(-1) : Scalar(1); m_s = 0; if(r) *r = m_c * p; } else if(p==Scalar(0)) { m_c = 0; - m_s = -q/internal::abs(q); - if(r) *r = internal::abs(q); + m_s = -q/abs(q); + if(r) *r = abs(q); } else { - RealScalar p1 = internal::norm1(p); - RealScalar q1 = internal::norm1(q); + RealScalar p1 = numext::norm1(p); + RealScalar q1 = numext::norm1(q); if(p1>=q1) { Scalar ps = p / p1; - RealScalar p2 = internal::abs2(ps); + RealScalar p2 = numext::abs2(ps); Scalar qs = q / p1; - RealScalar q2 = internal::abs2(qs); + RealScalar q2 = numext::abs2(qs); - RealScalar u = internal::sqrt(RealScalar(1) + q2/p2); - if(internal::real(p)<RealScalar(0)) + RealScalar u = sqrt(RealScalar(1) + q2/p2); + if(numext::real(p)<RealScalar(0)) u = -u; m_c = Scalar(1)/u; - m_s = -qs*internal::conj(ps)*(m_c/p2); + m_s = -qs*conj(ps)*(m_c/p2); if(r) *r = p * u; } else { Scalar ps = p / q1; - RealScalar p2 = internal::abs2(ps); + RealScalar p2 = numext::abs2(ps); Scalar qs = q / q1; - RealScalar q2 = internal::abs2(qs); + RealScalar q2 = numext::abs2(qs); - RealScalar u = q1 * internal::sqrt(p2 + q2); - if(internal::real(p)<RealScalar(0)) + RealScalar u = q1 * sqrt(p2 + q2); + if(numext::real(p)<RealScalar(0)) u = -u; - p1 = internal::abs(p); + p1 = abs(p); ps = p/p1; m_c = p1/u; - m_s = -internal::conj(ps) * (q/u); + m_s = -conj(ps) * (q/u); if(r) *r = ps * u; } } @@ -207,23 +214,24 @@ 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) { - + using std::sqrt; + using std::abs; if(q==Scalar(0)) { m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1); m_s = Scalar(0); - if(r) *r = internal::abs(p); + if(r) *r = abs(p); } else if(p==Scalar(0)) { m_c = Scalar(0); m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1); - if(r) *r = internal::abs(q); + if(r) *r = abs(q); } - else if(internal::abs(p) > internal::abs(q)) + else if(abs(p) > abs(q)) { Scalar t = q/p; - Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t)); + Scalar u = sqrt(Scalar(1) + numext::abs2(t)); if(p<Scalar(0)) u = -u; m_c = Scalar(1)/u; @@ -233,7 +241,7 @@ void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar else { Scalar t = p/q; - Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t)); + Scalar u = sqrt(Scalar(1) + numext::abs2(t)); if(q<Scalar(0)) u = -u; m_s = -Scalar(1)/u; @@ -303,6 +311,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 +329,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 + numext::conj(s) * yi; + y[i] = -s * xi + numext::conj(c) * yi; } Scalar* EIGEN_RESTRICT px = x + alignedStart; @@ -372,8 +385,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 + numext::conj(s) * yi; + y[i] = -s * xi + numext::conj(c) * yi; } } @@ -382,8 +395,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 +418,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 + numext::conj(s) * yi; + *y = -s * xi + numext::conj(c) * yi; x += incrx; y += incry; } |