From 4716040703be1ee906439385d20475dcddad5ce3 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 25 Oct 2010 10:15:22 -0400 Subject: bug #86 : use internal:: namespace instead of ei_ prefix --- Eigen/src/Jacobi/Jacobi.h | 163 +++++++++++++++++++++++----------------------- 1 file changed, 83 insertions(+), 80 deletions(-) (limited to 'Eigen/src/Jacobi') diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 39420bf42..58f6c98fb 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -29,7 +29,7 @@ /** \ingroup Jacobi_Module * \jacobi_module * \class JacobiRotation - * \brief Represents a rotation given by a cosine-sine pair. + * \brief Rotation given by a cosine-sine pair. * * This class represents a Jacobi or Givens rotation. * This is a 2D rotation in the plane \c J of angle \f$ \theta \f$ defined by @@ -63,15 +63,15 @@ template class JacobiRotation /** Concatenates two planar rotation */ JacobiRotation operator*(const JacobiRotation& other) { - return JacobiRotation(m_c * other.m_c - ei_conj(m_s) * other.m_s, - ei_conj(m_c * ei_conj(other.m_s) + ei_conj(m_s) * ei_conj(other.m_c))); + 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))); } /** Returns the transposed transformation */ - JacobiRotation transpose() const { return JacobiRotation(m_c, -ei_conj(m_s)); } + JacobiRotation transpose() const { return JacobiRotation(m_c, -internal::conj(m_s)); } /** Returns the adjoint transformation */ - JacobiRotation adjoint() const { return JacobiRotation(ei_conj(m_c), -m_s); } + JacobiRotation adjoint() const { return JacobiRotation(internal::conj(m_c), -m_s); } template bool makeJacobi(const MatrixBase&, typename Derived::Index p, typename Derived::Index q); @@ -80,8 +80,8 @@ template class JacobiRotation void makeGivens(const Scalar& p, const Scalar& q, Scalar* z=0); protected: - void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, ei_meta_true); - void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, ei_meta_false); + void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, internal::meta_true); + void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, internal::meta_false); Scalar m_c, m_s; }; @@ -103,8 +103,8 @@ bool JacobiRotation::makeJacobi(RealScalar x, Scalar y, RealScalar z) } else { - RealScalar tau = (x-z)/(RealScalar(2)*ei_abs(y)); - RealScalar w = ei_sqrt(ei_abs2(tau) + 1); + RealScalar tau = (x-z)/(RealScalar(2)*internal::abs(y)); + RealScalar w = internal::sqrt(internal::abs2(tau) + 1); RealScalar t; if(tau>0) { @@ -115,8 +115,8 @@ bool JacobiRotation::makeJacobi(RealScalar x, Scalar y, RealScalar z) t = RealScalar(1) / (tau - w); } RealScalar sign_t = t > 0 ? 1 : -1; - RealScalar n = RealScalar(1) / ei_sqrt(ei_abs2(t)+1); - m_s = - sign_t * (ei_conj(y) / ei_abs(y)) * ei_abs(t) * n; + RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+1); + m_s = - sign_t * (internal::conj(y) / internal::abs(y)) * internal::abs(t) * n; m_c = n; return true; } @@ -135,7 +135,7 @@ template template inline bool JacobiRotation::makeJacobi(const MatrixBase& m, typename Derived::Index p, typename Derived::Index q) { - return makeJacobi(ei_real(m.coeff(p,p)), m.coeff(p,q), ei_real(m.coeff(q,q))); + return makeJacobi(internal::real(m.coeff(p,p)), m.coeff(p,q), internal::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 @@ -157,60 +157,60 @@ inline bool JacobiRotation::makeJacobi(const MatrixBase& m, typ template void JacobiRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar* z) { - makeGivens(p, q, z, typename ei_meta_if::IsComplex, ei_meta_true, ei_meta_false>::ret()); + makeGivens(p, q, z, typename internal::meta_if::IsComplex, internal::meta_true, internal::meta_false>::ret()); } // specialization for complexes template -void JacobiRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, ei_meta_true) +void JacobiRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::meta_true) { if(q==Scalar(0)) { - m_c = ei_real(p)<0 ? Scalar(-1) : Scalar(1); + m_c = internal::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/ei_abs(q); - if(r) *r = ei_abs(q); + m_s = -q/internal::abs(q); + if(r) *r = internal::abs(q); } else { - RealScalar p1 = ei_norm1(p); - RealScalar q1 = ei_norm1(q); + RealScalar p1 = internal::norm1(p); + RealScalar q1 = internal::norm1(q); if(p1>=q1) { Scalar ps = p / p1; - RealScalar p2 = ei_abs2(ps); + RealScalar p2 = internal::abs2(ps); Scalar qs = q / p1; - RealScalar q2 = ei_abs2(qs); + RealScalar q2 = internal::abs2(qs); - RealScalar u = ei_sqrt(RealScalar(1) + q2/p2); - if(ei_real(p)::makeGivens(const Scalar& p, const Scalar& q, Scalar // specialization for reals template -void JacobiRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, ei_meta_false) +void JacobiRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::meta_false) { if(q==0) { m_c = p ei_abs(q)) + else if(internal::abs(p) > internal::abs(q)) { Scalar t = q/p; - Scalar u = ei_sqrt(Scalar(1) + ei_abs2(t)); + Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t)); if(p::makeGivens(const Scalar& p, const Scalar& q, Scalar else { Scalar t = p/q; - Scalar u = ei_sqrt(Scalar(1) + ei_abs2(t)); + Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t)); if(q::makeGivens(const Scalar& p, const Scalar& q, Scalar * * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ +namespace internal { template -void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation& j); +void apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation& j); +} /** \jacobi_module * Applies the rotation in the plane \a j to the rows \a p and \a q of \c *this, i.e., it computes B = J * B, * with \f$ B = \left ( \begin{array}{cc} \text{*this.row}(p) \\ \text{*this.row}(q) \end{array} \right ) \f$. * - * \sa class JacobiRotation, MatrixBase::applyOnTheRight(), ei_apply_rotation_in_the_plane() + * \sa class JacobiRotation, MatrixBase::applyOnTheRight(), internal::apply_rotation_in_the_plane() */ template template @@ -281,14 +283,14 @@ inline void MatrixBase::applyOnTheLeft(Index p, Index q, const JacobiRo { RowXpr x(this->row(p)); RowXpr y(this->row(q)); - ei_apply_rotation_in_the_plane(x, y, j); + internal::apply_rotation_in_the_plane(x, y, j); } /** \ingroup Jacobi_Module * Applies the rotation in the plane \a j to the columns \a p and \a q of \c *this, i.e., it computes B = B * J * with \f$ B = \left ( \begin{array}{cc} \text{*this.col}(p) & \text{*this.col}(q) \end{array} \right ) \f$. * - * \sa class JacobiRotation, MatrixBase::applyOnTheLeft(), ei_apply_rotation_in_the_plane() + * \sa class JacobiRotation, MatrixBase::applyOnTheLeft(), internal::apply_rotation_in_the_plane() */ template template @@ -296,18 +298,18 @@ inline void MatrixBase::applyOnTheRight(Index p, Index q, const JacobiR { ColXpr x(this->col(p)); ColXpr y(this->col(q)); - ei_apply_rotation_in_the_plane(x, y, j.transpose()); + internal::apply_rotation_in_the_plane(x, y, j.transpose()); } - +namespace internal { template -void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation& j) +void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation& j) { typedef typename VectorX::Index Index; typedef typename VectorX::Scalar Scalar; - enum { PacketSize = ei_packet_traits::size }; - typedef typename ei_packet_traits::type Packet; - ei_assert(_x.size() == _y.size()); + enum { PacketSize = packet_traits::size }; + typedef typename packet_traits::type Packet; + eigen_assert(_x.size() == _y.size()); Index size = _x.size(); Index incrx = _x.innerStride(); Index incry = _y.innerStride(); @@ -324,32 +326,32 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& // both vectors are sequentially stored in memory => vectorization enum { Peeling = 2 }; - Index alignedStart = ei_first_aligned(y, size); + Index alignedStart = first_aligned(y, size); Index alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize; - const Packet pc = ei_pset1(j.c()); - const Packet ps = ei_pset1(j.s()); - ei_conj_helper::IsComplex,false> pcj; + const Packet pc = pset1(j.c()); + const Packet ps = pset1(j.s()); + conj_helper::IsComplex,false> pcj; for(Index i=0; i(px); - Packet yi = ei_pload(py); - ei_pstore(px, ei_padd(ei_pmul(pc,xi),pcj.pmul(ps,yi))); - ei_pstore(py, ei_psub(pcj.pmul(pc,yi),ei_pmul(ps,xi))); + Packet xi = pload(px); + Packet yi = pload(py); + pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); + pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); px += PacketSize; py += PacketSize; } @@ -359,23 +361,23 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& Index peelingEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize); for(Index i=alignedStart; i(px); - Packet xi1 = ei_ploadu(px+PacketSize); - Packet yi = ei_pload (py); - Packet yi1 = ei_pload (py+PacketSize); - ei_pstoreu(px, ei_padd(ei_pmul(pc,xi),pcj.pmul(ps,yi))); - ei_pstoreu(px+PacketSize, ei_padd(ei_pmul(pc,xi1),pcj.pmul(ps,yi1))); - ei_pstore (py, ei_psub(pcj.pmul(pc,yi),ei_pmul(ps,xi))); - ei_pstore (py+PacketSize, ei_psub(pcj.pmul(pc,yi1),ei_pmul(ps,xi1))); + Packet xi = ploadu(px); + Packet xi1 = ploadu(px+PacketSize); + Packet yi = pload (py); + Packet yi1 = pload (py+PacketSize); + pstoreu(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); + pstoreu(px+PacketSize, padd(pmul(pc,xi1),pcj.pmul(ps,yi1))); + pstore (py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); + pstore (py+PacketSize, psub(pcj.pmul(pc,yi1),pmul(ps,xi1))); px += Peeling*PacketSize; py += Peeling*PacketSize; } if(alignedEnd!=peelingEnd) { - Packet xi = ei_ploadu(x+peelingEnd); - Packet yi = ei_pload (y+peelingEnd); - ei_pstoreu(x+peelingEnd, ei_padd(ei_pmul(pc,xi),pcj.pmul(ps,yi))); - ei_pstore (y+peelingEnd, ei_psub(pcj.pmul(pc,yi),ei_pmul(ps,xi))); + Packet xi = ploadu(x+peelingEnd); + Packet yi = pload (y+peelingEnd); + pstoreu(x+peelingEnd, padd(pmul(pc,xi),pcj.pmul(ps,yi))); + pstore (y+peelingEnd, psub(pcj.pmul(pc,yi),pmul(ps,xi))); } } @@ -383,8 +385,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Scalar xi = x[i]; Scalar yi = y[i]; - x[i] = j.c() * xi + ei_conj(j.s()) * yi; - y[i] = -j.s() * xi + ei_conj(j.c()) * yi; + x[i] = j.c() * xi + conj(j.s()) * yi; + y[i] = -j.s() * xi + conj(j.c()) * yi; } } @@ -393,17 +395,17 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& (VectorX::Flags & VectorY::Flags & PacketAccessBit) && (VectorX::Flags & VectorY::Flags & AlignedBit)) { - const Packet pc = ei_pset1(j.c()); - const Packet ps = ei_pset1(j.s()); - ei_conj_helper::IsComplex,false> pcj; + const Packet pc = pset1(j.c()); + const Packet ps = pset1(j.s()); + conj_helper::IsComplex,false> pcj; Scalar* EIGEN_RESTRICT px = x; Scalar* EIGEN_RESTRICT py = y; for(Index i=0; i(px); - Packet yi = ei_pload(py); - ei_pstore(px, ei_padd(ei_pmul(pc,xi),pcj.pmul(ps,yi))); - ei_pstore(py, ei_psub(pcj.pmul(pc,yi),ei_pmul(ps,xi))); + Packet xi = pload(px); + Packet yi = pload(py); + pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); + pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); px += PacketSize; py += PacketSize; } @@ -416,12 +418,13 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Scalar xi = *x; Scalar yi = *y; - *x = j.c() * xi + ei_conj(j.s()) * yi; - *y = -j.s() * xi + ei_conj(j.c()) * yi; + *x = j.c() * xi + conj(j.s()) * yi; + *y = -j.s() * xi + conj(j.c()) * yi; x += incrx; y += incry; } } } +} #endif // EIGEN_JACOBI_H -- cgit v1.2.3