Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/eigen3/Eigen/src/Jacobi/Jacobi.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010 
0011 #ifndef EIGEN_JACOBI_H
0012 #define EIGEN_JACOBI_H
0013 
0014 namespace Eigen {
0015 
0016 /** \ingroup Jacobi_Module
0017   * \jacobi_module
0018   * \class JacobiRotation
0019   * \brief Rotation given by a cosine-sine pair.
0020   *
0021   * This class represents a Jacobi or Givens rotation.
0022   * This is a 2D rotation in the plane \c J of angle \f$ \theta \f$ defined by
0023   * its cosine \c c and sine \c s as follow:
0024   * \f$ J = \left ( \begin{array}{cc} c & \overline s \\ -s  & \overline c \end{array} \right ) \f$
0025   *
0026   * You can apply the respective counter-clockwise rotation to a column vector \c v by
0027   * applying its adjoint on the left: \f$ v = J^* v \f$ that translates to the following Eigen code:
0028   * \code
0029   * v.applyOnTheLeft(J.adjoint());
0030   * \endcode
0031   *
0032   * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
0033   */
0034 template<typename Scalar> class JacobiRotation
0035 {
0036   public:
0037     typedef typename NumTraits<Scalar>::Real RealScalar;
0038 
0039     /** Default constructor without any initialization. */
0040     EIGEN_DEVICE_FUNC
0041     JacobiRotation() {}
0042 
0043     /** Construct a planar rotation from a cosine-sine pair (\a c, \c s). */
0044     EIGEN_DEVICE_FUNC
0045     JacobiRotation(const Scalar& c, const Scalar& s) : m_c(c), m_s(s) {}
0046 
0047     EIGEN_DEVICE_FUNC Scalar& c() { return m_c; }
0048     EIGEN_DEVICE_FUNC Scalar c() const { return m_c; }
0049     EIGEN_DEVICE_FUNC Scalar& s() { return m_s; }
0050     EIGEN_DEVICE_FUNC Scalar s() const { return m_s; }
0051 
0052     /** Concatenates two planar rotation */
0053     EIGEN_DEVICE_FUNC
0054     JacobiRotation operator*(const JacobiRotation& other)
0055     {
0056       using numext::conj;
0057       return JacobiRotation(m_c * other.m_c - conj(m_s) * other.m_s,
0058                             conj(m_c * conj(other.m_s) + conj(m_s) * conj(other.m_c)));
0059     }
0060 
0061     /** Returns the transposed transformation */
0062     EIGEN_DEVICE_FUNC
0063     JacobiRotation transpose() const { using numext::conj; return JacobiRotation(m_c, -conj(m_s)); }
0064 
0065     /** Returns the adjoint transformation */
0066     EIGEN_DEVICE_FUNC
0067     JacobiRotation adjoint() const { using numext::conj; return JacobiRotation(conj(m_c), -m_s); }
0068 
0069     template<typename Derived>
0070     EIGEN_DEVICE_FUNC
0071     bool makeJacobi(const MatrixBase<Derived>&, Index p, Index q);
0072     EIGEN_DEVICE_FUNC
0073     bool makeJacobi(const RealScalar& x, const Scalar& y, const RealScalar& z);
0074 
0075     EIGEN_DEVICE_FUNC
0076     void makeGivens(const Scalar& p, const Scalar& q, Scalar* r=0);
0077 
0078   protected:
0079     EIGEN_DEVICE_FUNC
0080     void makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::true_type);
0081     EIGEN_DEVICE_FUNC
0082     void makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::false_type);
0083 
0084     Scalar m_c, m_s;
0085 };
0086 
0087 /** Makes \c *this as a Jacobi rotation \a J such that applying \a J on both the right and left sides of the selfadjoint 2x2 matrix
0088   * \f$ B = \left ( \begin{array}{cc} x & y \\ \overline y & z \end{array} \right )\f$ yields a diagonal matrix \f$ A = J^* B J \f$
0089   *
0090   * \sa MatrixBase::makeJacobi(const MatrixBase<Derived>&, Index, Index), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
0091   */
0092 template<typename Scalar>
0093 EIGEN_DEVICE_FUNC
0094 bool JacobiRotation<Scalar>::makeJacobi(const RealScalar& x, const Scalar& y, const RealScalar& z)
0095 {
0096   using std::sqrt;
0097   using std::abs;
0098 
0099   RealScalar deno = RealScalar(2)*abs(y);
0100   if(deno < (std::numeric_limits<RealScalar>::min)())
0101   {
0102     m_c = Scalar(1);
0103     m_s = Scalar(0);
0104     return false;
0105   }
0106   else
0107   {
0108     RealScalar tau = (x-z)/deno;
0109     RealScalar w = sqrt(numext::abs2(tau) + RealScalar(1));
0110     RealScalar t;
0111     if(tau>RealScalar(0))
0112     {
0113       t = RealScalar(1) / (tau + w);
0114     }
0115     else
0116     {
0117       t = RealScalar(1) / (tau - w);
0118     }
0119     RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
0120     RealScalar n = RealScalar(1) / sqrt(numext::abs2(t)+RealScalar(1));
0121     m_s = - sign_t * (numext::conj(y) / abs(y)) * abs(t) * n;
0122     m_c = n;
0123     return true;
0124   }
0125 }
0126 
0127 /** Makes \c *this as a Jacobi rotation \c J such that applying \a J on both the right and left sides of the 2x2 selfadjoint matrix
0128   * \f$ B = \left ( \begin{array}{cc} \text{this}_{pp} & \text{this}_{pq} \\ (\text{this}_{pq})^* & \text{this}_{qq} \end{array} \right )\f$ yields
0129   * a diagonal matrix \f$ A = J^* B J \f$
0130   *
0131   * Example: \include Jacobi_makeJacobi.cpp
0132   * Output: \verbinclude Jacobi_makeJacobi.out
0133   *
0134   * \sa JacobiRotation::makeJacobi(RealScalar, Scalar, RealScalar), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
0135   */
0136 template<typename Scalar>
0137 template<typename Derived>
0138 EIGEN_DEVICE_FUNC
0139 inline bool JacobiRotation<Scalar>::makeJacobi(const MatrixBase<Derived>& m, Index p, Index q)
0140 {
0141   return makeJacobi(numext::real(m.coeff(p,p)), m.coeff(p,q), numext::real(m.coeff(q,q)));
0142 }
0143 
0144 /** Makes \c *this as a Givens rotation \c G such that applying \f$ G^* \f$ to the left of the vector
0145   * \f$ V = \left ( \begin{array}{c} p \\ q \end{array} \right )\f$ yields:
0146   * \f$ G^* V = \left ( \begin{array}{c} r \\ 0 \end{array} \right )\f$.
0147   *
0148   * The value of \a r is returned if \a r is not null (the default is null).
0149   * Also note that G is built such that the cosine is always real.
0150   *
0151   * Example: \include Jacobi_makeGivens.cpp
0152   * Output: \verbinclude Jacobi_makeGivens.out
0153   *
0154   * This function implements the continuous Givens rotation generation algorithm
0155   * found in Anderson (2000), Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem.
0156   * LAPACK Working Note 150, University of Tennessee, UT-CS-00-454, December 4, 2000.
0157   *
0158   * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
0159   */
0160 template<typename Scalar>
0161 EIGEN_DEVICE_FUNC
0162 void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r)
0163 {
0164   makeGivens(p, q, r, typename internal::conditional<NumTraits<Scalar>::IsComplex, internal::true_type, internal::false_type>::type());
0165 }
0166 
0167 
0168 // specialization for complexes
0169 template<typename Scalar>
0170 EIGEN_DEVICE_FUNC
0171 void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::true_type)
0172 {
0173   using std::sqrt;
0174   using std::abs;
0175   using numext::conj;
0176 
0177   if(q==Scalar(0))
0178   {
0179     m_c = numext::real(p)<0 ? Scalar(-1) : Scalar(1);
0180     m_s = 0;
0181     if(r) *r = m_c * p;
0182   }
0183   else if(p==Scalar(0))
0184   {
0185     m_c = 0;
0186     m_s = -q/abs(q);
0187     if(r) *r = abs(q);
0188   }
0189   else
0190   {
0191     RealScalar p1 = numext::norm1(p);
0192     RealScalar q1 = numext::norm1(q);
0193     if(p1>=q1)
0194     {
0195       Scalar ps = p / p1;
0196       RealScalar p2 = numext::abs2(ps);
0197       Scalar qs = q / p1;
0198       RealScalar q2 = numext::abs2(qs);
0199 
0200       RealScalar u = sqrt(RealScalar(1) + q2/p2);
0201       if(numext::real(p)<RealScalar(0))
0202         u = -u;
0203 
0204       m_c = Scalar(1)/u;
0205       m_s = -qs*conj(ps)*(m_c/p2);
0206       if(r) *r = p * u;
0207     }
0208     else
0209     {
0210       Scalar ps = p / q1;
0211       RealScalar p2 = numext::abs2(ps);
0212       Scalar qs = q / q1;
0213       RealScalar q2 = numext::abs2(qs);
0214 
0215       RealScalar u = q1 * sqrt(p2 + q2);
0216       if(numext::real(p)<RealScalar(0))
0217         u = -u;
0218 
0219       p1 = abs(p);
0220       ps = p/p1;
0221       m_c = p1/u;
0222       m_s = -conj(ps) * (q/u);
0223       if(r) *r = ps * u;
0224     }
0225   }
0226 }
0227 
0228 // specialization for reals
0229 template<typename Scalar>
0230 EIGEN_DEVICE_FUNC
0231 void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::false_type)
0232 {
0233   using std::sqrt;
0234   using std::abs;
0235   if(q==Scalar(0))
0236   {
0237     m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
0238     m_s = Scalar(0);
0239     if(r) *r = abs(p);
0240   }
0241   else if(p==Scalar(0))
0242   {
0243     m_c = Scalar(0);
0244     m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
0245     if(r) *r = abs(q);
0246   }
0247   else if(abs(p) > abs(q))
0248   {
0249     Scalar t = q/p;
0250     Scalar u = sqrt(Scalar(1) + numext::abs2(t));
0251     if(p<Scalar(0))
0252       u = -u;
0253     m_c = Scalar(1)/u;
0254     m_s = -t * m_c;
0255     if(r) *r = p * u;
0256   }
0257   else
0258   {
0259     Scalar t = p/q;
0260     Scalar u = sqrt(Scalar(1) + numext::abs2(t));
0261     if(q<Scalar(0))
0262       u = -u;
0263     m_s = -Scalar(1)/u;
0264     m_c = -t * m_s;
0265     if(r) *r = q * u;
0266   }
0267 
0268 }
0269 
0270 /****************************************************************************************
0271 *   Implementation of MatrixBase methods
0272 ****************************************************************************************/
0273 
0274 namespace internal {
0275 /** \jacobi_module
0276   * Applies the clock wise 2D rotation \a j to the set of 2D vectors of coordinates \a x and \a y:
0277   * \f$ \left ( \begin{array}{cc} x \\ y \end{array} \right )  =  J \left ( \begin{array}{cc} x \\ y \end{array} \right ) \f$
0278   *
0279   * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
0280   */
0281 template<typename VectorX, typename VectorY, typename OtherScalar>
0282 EIGEN_DEVICE_FUNC
0283 void apply_rotation_in_the_plane(DenseBase<VectorX>& xpr_x, DenseBase<VectorY>& xpr_y, const JacobiRotation<OtherScalar>& j);
0284 }
0285 
0286 /** \jacobi_module
0287   * 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,
0288   * with \f$ B = \left ( \begin{array}{cc} \text{*this.row}(p) \\ \text{*this.row}(q) \end{array} \right ) \f$.
0289   *
0290   * \sa class JacobiRotation, MatrixBase::applyOnTheRight(), internal::apply_rotation_in_the_plane()
0291   */
0292 template<typename Derived>
0293 template<typename OtherScalar>
0294 EIGEN_DEVICE_FUNC
0295 inline void MatrixBase<Derived>::applyOnTheLeft(Index p, Index q, const JacobiRotation<OtherScalar>& j)
0296 {
0297   RowXpr x(this->row(p));
0298   RowXpr y(this->row(q));
0299   internal::apply_rotation_in_the_plane(x, y, j);
0300 }
0301 
0302 /** \ingroup Jacobi_Module
0303   * 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
0304   * with \f$ B = \left ( \begin{array}{cc} \text{*this.col}(p) & \text{*this.col}(q) \end{array} \right ) \f$.
0305   *
0306   * \sa class JacobiRotation, MatrixBase::applyOnTheLeft(), internal::apply_rotation_in_the_plane()
0307   */
0308 template<typename Derived>
0309 template<typename OtherScalar>
0310 EIGEN_DEVICE_FUNC
0311 inline void MatrixBase<Derived>::applyOnTheRight(Index p, Index q, const JacobiRotation<OtherScalar>& j)
0312 {
0313   ColXpr x(this->col(p));
0314   ColXpr y(this->col(q));
0315   internal::apply_rotation_in_the_plane(x, y, j.transpose());
0316 }
0317 
0318 namespace internal {
0319 
0320 template<typename Scalar, typename OtherScalar,
0321          int SizeAtCompileTime, int MinAlignment, bool Vectorizable>
0322 struct apply_rotation_in_the_plane_selector
0323 {
0324   static EIGEN_DEVICE_FUNC
0325   inline void run(Scalar *x, Index incrx, Scalar *y, Index incry, Index size, OtherScalar c, OtherScalar s)
0326   {
0327     for(Index i=0; i<size; ++i)
0328     {
0329       Scalar xi = *x;
0330       Scalar yi = *y;
0331       *x =  c * xi + numext::conj(s) * yi;
0332       *y = -s * xi + numext::conj(c) * yi;
0333       x += incrx;
0334       y += incry;
0335     }
0336   }
0337 };
0338 
0339 template<typename Scalar, typename OtherScalar,
0340          int SizeAtCompileTime, int MinAlignment>
0341 struct apply_rotation_in_the_plane_selector<Scalar,OtherScalar,SizeAtCompileTime,MinAlignment,true /* vectorizable */>
0342 {
0343   static inline void run(Scalar *x, Index incrx, Scalar *y, Index incry, Index size, OtherScalar c, OtherScalar s)
0344   {
0345     enum {
0346       PacketSize = packet_traits<Scalar>::size,
0347       OtherPacketSize = packet_traits<OtherScalar>::size
0348     };
0349     typedef typename packet_traits<Scalar>::type Packet;
0350     typedef typename packet_traits<OtherScalar>::type OtherPacket;
0351 
0352     /*** dynamic-size vectorized paths ***/
0353     if(SizeAtCompileTime == Dynamic && ((incrx==1 && incry==1) || PacketSize == 1))
0354     {
0355       // both vectors are sequentially stored in memory => vectorization
0356       enum { Peeling = 2 };
0357 
0358       Index alignedStart = internal::first_default_aligned(y, size);
0359       Index alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize;
0360 
0361       const OtherPacket pc = pset1<OtherPacket>(c);
0362       const OtherPacket ps = pset1<OtherPacket>(s);
0363       conj_helper<OtherPacket,Packet,NumTraits<OtherScalar>::IsComplex,false> pcj;
0364       conj_helper<OtherPacket,Packet,false,false> pm;
0365 
0366       for(Index i=0; i<alignedStart; ++i)
0367       {
0368         Scalar xi = x[i];
0369         Scalar yi = y[i];
0370         x[i] =  c * xi + numext::conj(s) * yi;
0371         y[i] = -s * xi + numext::conj(c) * yi;
0372       }
0373 
0374       Scalar* EIGEN_RESTRICT px = x + alignedStart;
0375       Scalar* EIGEN_RESTRICT py = y + alignedStart;
0376 
0377       if(internal::first_default_aligned(x, size)==alignedStart)
0378       {
0379         for(Index i=alignedStart; i<alignedEnd; i+=PacketSize)
0380         {
0381           Packet xi = pload<Packet>(px);
0382           Packet yi = pload<Packet>(py);
0383           pstore(px, padd(pm.pmul(pc,xi),pcj.pmul(ps,yi)));
0384           pstore(py, psub(pcj.pmul(pc,yi),pm.pmul(ps,xi)));
0385           px += PacketSize;
0386           py += PacketSize;
0387         }
0388       }
0389       else
0390       {
0391         Index peelingEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize);
0392         for(Index i=alignedStart; i<peelingEnd; i+=Peeling*PacketSize)
0393         {
0394           Packet xi   = ploadu<Packet>(px);
0395           Packet xi1  = ploadu<Packet>(px+PacketSize);
0396           Packet yi   = pload <Packet>(py);
0397           Packet yi1  = pload <Packet>(py+PacketSize);
0398           pstoreu(px, padd(pm.pmul(pc,xi),pcj.pmul(ps,yi)));
0399           pstoreu(px+PacketSize, padd(pm.pmul(pc,xi1),pcj.pmul(ps,yi1)));
0400           pstore (py, psub(pcj.pmul(pc,yi),pm.pmul(ps,xi)));
0401           pstore (py+PacketSize, psub(pcj.pmul(pc,yi1),pm.pmul(ps,xi1)));
0402           px += Peeling*PacketSize;
0403           py += Peeling*PacketSize;
0404         }
0405         if(alignedEnd!=peelingEnd)
0406         {
0407           Packet xi = ploadu<Packet>(x+peelingEnd);
0408           Packet yi = pload <Packet>(y+peelingEnd);
0409           pstoreu(x+peelingEnd, padd(pm.pmul(pc,xi),pcj.pmul(ps,yi)));
0410           pstore (y+peelingEnd, psub(pcj.pmul(pc,yi),pm.pmul(ps,xi)));
0411         }
0412       }
0413 
0414       for(Index i=alignedEnd; i<size; ++i)
0415       {
0416         Scalar xi = x[i];
0417         Scalar yi = y[i];
0418         x[i] =  c * xi + numext::conj(s) * yi;
0419         y[i] = -s * xi + numext::conj(c) * yi;
0420       }
0421     }
0422 
0423     /*** fixed-size vectorized path ***/
0424     else if(SizeAtCompileTime != Dynamic && MinAlignment>0) // FIXME should be compared to the required alignment
0425     {
0426       const OtherPacket pc = pset1<OtherPacket>(c);
0427       const OtherPacket ps = pset1<OtherPacket>(s);
0428       conj_helper<OtherPacket,Packet,NumTraits<OtherPacket>::IsComplex,false> pcj;
0429       conj_helper<OtherPacket,Packet,false,false> pm;
0430       Scalar* EIGEN_RESTRICT px = x;
0431       Scalar* EIGEN_RESTRICT py = y;
0432       for(Index i=0; i<size; i+=PacketSize)
0433       {
0434         Packet xi = pload<Packet>(px);
0435         Packet yi = pload<Packet>(py);
0436         pstore(px, padd(pm.pmul(pc,xi),pcj.pmul(ps,yi)));
0437         pstore(py, psub(pcj.pmul(pc,yi),pm.pmul(ps,xi)));
0438         px += PacketSize;
0439         py += PacketSize;
0440       }
0441     }
0442 
0443     /*** non-vectorized path ***/
0444     else
0445     {
0446       apply_rotation_in_the_plane_selector<Scalar,OtherScalar,SizeAtCompileTime,MinAlignment,false>::run(x,incrx,y,incry,size,c,s);
0447     }
0448   }
0449 };
0450 
0451 template<typename VectorX, typename VectorY, typename OtherScalar>
0452 EIGEN_DEVICE_FUNC
0453 void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(DenseBase<VectorX>& xpr_x, DenseBase<VectorY>& xpr_y, const JacobiRotation<OtherScalar>& j)
0454 {
0455   typedef typename VectorX::Scalar Scalar;
0456   const bool Vectorizable =    (int(VectorX::Flags) & int(VectorY::Flags) & PacketAccessBit)
0457                             && (int(packet_traits<Scalar>::size) == int(packet_traits<OtherScalar>::size));
0458 
0459   eigen_assert(xpr_x.size() == xpr_y.size());
0460   Index size = xpr_x.size();
0461   Index incrx = xpr_x.derived().innerStride();
0462   Index incry = xpr_y.derived().innerStride();
0463 
0464   Scalar* EIGEN_RESTRICT x = &xpr_x.derived().coeffRef(0);
0465   Scalar* EIGEN_RESTRICT y = &xpr_y.derived().coeffRef(0);
0466 
0467   OtherScalar c = j.c();
0468   OtherScalar s = j.s();
0469   if (c==OtherScalar(1) && s==OtherScalar(0))
0470     return;
0471 
0472   apply_rotation_in_the_plane_selector<
0473     Scalar,OtherScalar,
0474     VectorX::SizeAtCompileTime,
0475     EIGEN_PLAIN_ENUM_MIN(evaluator<VectorX>::Alignment, evaluator<VectorY>::Alignment),
0476     Vectorizable>::run(x,incrx,y,incry,size,c,s);
0477 }
0478 
0479 } // end namespace internal
0480 
0481 } // end namespace Eigen
0482 
0483 #endif // EIGEN_JACOBI_H
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!