From 6e4e94ff3266f8d85adbfe6556e612a2ff9a7e68 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 31 Aug 2009 22:26:15 -0400 Subject: * JacobiSVD: - support complex numbers - big rewrite of the 2x2 kernel, much more robust * Jacobi: - fix weirdness in initial design, e.g. applyJacobiOnTheRight actually did the inverse transformation - fully support complex numbers - fix logic to decide whether to vectorize - remove several clumsy methods fix for complex numbers --- Eigen/src/Jacobi/Jacobi.h | 113 +++++++++++++++++++--------------------------- 1 file changed, 46 insertions(+), 67 deletions(-) (limited to 'Eigen/src/Jacobi') diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 76f0800fe..96f08d54a 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -33,19 +33,20 @@ * * \sa MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight() */ -template -void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s); +template +void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, JacobiScalar c, JacobiScalar s); /** Applies a rotation in the plane defined by \a c, \a s to the rows \a p and \a q of \c *this. * More precisely, it computes B = J' * B, with J = [c s ; -s' c] and B = [ *this.row(p) ; *this.row(q) ] * \sa MatrixBase::applyJacobiOnTheRight(), ei_apply_rotation_in_the_plane() */ template -inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) +template +inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, JacobiScalar c, JacobiScalar s) { RowXpr x(row(p)); RowXpr y(row(q)); - ei_apply_rotation_in_the_plane(x, y, ei_conj(c), ei_conj(s)); + ei_apply_rotation_in_the_plane(x, y, c, s); } /** Applies a rotation in the plane defined by \a c, \a s to the columns \a p and \a q of \c *this. @@ -53,23 +54,25 @@ inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Sc * \sa MatrixBase::applyJacobiOnTheLeft(), ei_apply_rotation_in_the_plane() */ template -inline void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) +template +inline void MatrixBase::applyJacobiOnTheRight(int p, int q, JacobiScalar c, JacobiScalar s) { ColXpr x(col(p)); ColXpr y(col(q)); - ei_apply_rotation_in_the_plane(x, y, c, s); + ei_apply_rotation_in_the_plane(x, y, c, -ei_conj(s)); } /** Computes the cosine-sine pair (\a c, \a s) such that its associated - * rotation \f$ J = ( \begin{array}{cc} c & s \\ -s' c \end{array} )\f$ + * rotation \f$ J = ( \begin{array}{cc} c & \overline s \\ -s & \overline c \end{array} )\f$ * applied to both the right and left of the 2x2 matrix * \f$ B = ( \begin{array}{cc} x & y \\ * & z \end{array} )\f$ yields - * a diagonal matrix A: \f$ A = J' B J \f$ + * a diagonal matrix A: \f$ A = J^* B J \f$ */ template -bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) +bool ei_makeJacobi(typename NumTraits::Real x, Scalar y, typename NumTraits::Real z, Scalar *c, Scalar *s) { - if(y == 0) + typedef typename NumTraits::Real RealScalar; + if(y == Scalar(0)) { *c = Scalar(1); *s = Scalar(0); @@ -77,15 +80,21 @@ bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) } else { - Scalar tau = (z - x) / (2 * y); - Scalar w = ei_sqrt(1 + ei_abs2(tau)); - Scalar t; + RealScalar tau = (x-z)/(RealScalar(2)*ei_abs(y)); + RealScalar w = ei_sqrt(ei_abs2(tau) + 1); + RealScalar t; if(tau>0) - t = Scalar(1) / (tau + w); + { + t = RealScalar(1) / (tau + w); + } else - t = Scalar(1) / (tau - w); - *c = Scalar(1) / ei_sqrt(1 + ei_abs2(t)); - *s = *c * t; + { + t = RealScalar(1) / (tau - w); + } + RealScalar sign_t = t > 0 ? 1 : -1; + RealScalar n = RealScalar(1) / ei_sqrt(ei_abs2(t)+1); + *s = - sign_t * (ei_conj(y) / ei_abs(y)) * ei_abs(t) * n; + *c = n; return true; } } @@ -93,41 +102,11 @@ bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) template inline bool MatrixBase::makeJacobi(int p, int q, Scalar *c, Scalar *s) const { - return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), c, s); -} - -template -inline bool MatrixBase::makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const -{ - return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(q,p)), - ei_conj(coeff(p,p))*coeff(p,q) + ei_conj(coeff(q,p))*coeff(q,q), - ei_abs2(coeff(p,q)) + ei_abs2(coeff(q,q)), - c,s); -} - -template -inline bool MatrixBase::makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const -{ - return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(p,q)), - ei_conj(coeff(q,p))*coeff(p,p) + ei_conj(coeff(q,q))*coeff(p,q), - ei_abs2(coeff(q,p)) + ei_abs2(coeff(q,q)), - c,s); -} - -template -inline void ei_normalizeJacobi(Scalar *c, Scalar *s, const Scalar& x, const Scalar& y) -{ - Scalar a = x * *c - y * *s; - Scalar b = x * *s + y * *c; - if(ei_abs(b)>ei_abs(a)) { - Scalar x = *c; - *c = -*s; - *s = x; - } + return ei_makeJacobi(ei_real(coeff(p,p)), coeff(p,q), ei_real(coeff(q,q)), c, s); } -template -void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s) +template +void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, JacobiScalar c, JacobiScalar s) { typedef typename VectorX::Scalar Scalar; ei_assert(_x.size() == _y.size()); @@ -138,7 +117,7 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0); Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0); - if (incrx==1 && incry==1) + if((VectorX::Flags & VectorY::Flags & PacketAccessBit) && incrx==1 && incry==1) { // both vectors are sequentially stored in memory => vectorization typedef typename ei_packet_traits::type Packet; @@ -147,16 +126,16 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& int alignedStart = ei_alignmentOffset(y, size); int alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize; - const Packet pc = ei_pset1(c); - const Packet ps = ei_pset1(s); + const Packet pc = ei_pset1(Scalar(c)); + const Packet ps = ei_pset1(Scalar(s)); ei_conj_helper::IsComplex,false> cj; for(int i=0; i