From b83654b5d01155c4ea875090f0779b99241a91a4 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 2 Sep 2009 15:04:10 +0200 Subject: [PATCH] * rename JacobiRotation => PlanarRotation * move the makeJacobi and make_givens_* to PlanarRotation * rename applyJacobi* => apply* --- Eigen/src/Core/MathFunctions.h | 4 + Eigen/src/Core/MatrixBase.h | 9 +- Eigen/src/Core/util/ForwardDeclarations.h | 2 +- Eigen/src/Jacobi/Jacobi.h | 205 +++++++++++++++------- Eigen/src/QR/ComplexSchur.h | 59 ++----- Eigen/src/QR/SelfAdjointEigenSolver.h | 45 ++--- Eigen/src/SVD/JacobiSVD.h | 28 +-- Eigen/src/SVD/SVD.h | 10 +- 8 files changed, 194 insertions(+), 168 deletions(-) diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 1570f01e0..40edf4a3c 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -116,6 +116,7 @@ inline float ei_imag(float) { return 0.f; } inline float ei_conj(float x) { return x; } inline float ei_abs(float x) { return std::abs(x); } inline float ei_abs2(float x) { return x*x; } +inline float ei_norm1(float x) { return ei_abs(x); } inline float ei_sqrt(float x) { return std::sqrt(x); } inline float ei_exp(float x) { return std::exp(x); } inline float ei_log(float x) { return std::log(x); } @@ -164,6 +165,7 @@ inline double ei_imag(double) { return 0.; } inline double ei_conj(double x) { return x; } inline double ei_abs(double x) { return std::abs(x); } inline double ei_abs2(double x) { return x*x; } +inline double ei_norm1(double x) { return ei_abs(x); } inline double ei_sqrt(double x) { return std::sqrt(x); } inline double ei_exp(double x) { return std::exp(x); } inline double ei_log(double x) { return std::log(x); } @@ -212,6 +214,7 @@ inline float& ei_imag_ref(std::complex& x) { return reinterpret_cast ei_conj(const std::complex& x) { return std::conj(x); } inline float ei_abs(const std::complex& x) { return std::abs(x); } inline float ei_abs2(const std::complex& x) { return std::norm(x); } +inline float ei_norm1(const std::complex &x) { return(ei_abs(x.real()) + ei_abs(x.imag())); } inline std::complex ei_exp(std::complex x) { return std::exp(x); } inline std::complex ei_sin(std::complex x) { return std::sin(x); } inline std::complex ei_cos(std::complex x) { return std::cos(x); } @@ -248,6 +251,7 @@ inline double& ei_imag_ref(std::complex& x) { return reinterpret_cast ei_conj(const std::complex& x) { return std::conj(x); } inline double ei_abs(const std::complex& x) { return std::abs(x); } inline double ei_abs2(const std::complex& x) { return std::norm(x); } +inline double ei_norm1(const std::complex &x) { return(ei_abs(x.real()) + ei_abs(x.imag())); } inline std::complex ei_exp(std::complex x) { return std::exp(x); } inline std::complex ei_sin(std::complex x) { return std::sin(x); } inline std::complex ei_cos(std::complex x) { return std::cos(x); } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index c2945ab46..9ac964168 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -803,11 +803,10 @@ template class MatrixBase ///////// Jacobi module ///////// - template - void applyJacobiOnTheLeft(int p, int q, const JacobiRotation& j); - template - void applyJacobiOnTheRight(int p, int q, const JacobiRotation& j); - bool makeJacobi(int p, int q, JacobiRotation *j) const; + template + void applyOnTheLeft(int p, int q, const PlanarRotation& j); + template + void applyOnTheRight(int p, int q, const PlanarRotation& j); #ifdef EIGEN_MATRIXBASE_PLUGIN #include EIGEN_MATRIXBASE_PLUGIN diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 18d3af7c5..c5f27d80b 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -123,7 +123,7 @@ template class SVD; template class JacobiSVD; template class LLT; template class LDLT; -template class JacobiRotation; +template class PlanarRotation; // Geometry module: template class RotationBase; diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 24fb7e782..c38d4583f 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -27,97 +27,72 @@ #define EIGEN_JACOBI_H /** \ingroup Jacobi - * \class JacobiRotation + * \class PlanarRotation * \brief Represents a rotation in the plane from a cosine-sine pair. * - * This class represents a Jacobi rotation which is also known as a Givens rotation. + * This class represents a Jacobi or Givens rotation. * This is a 2D clock-wise rotation in the plane \c J of angle \f$ \theta \f$ defined by * its cosine \c c and sine \c s as follow: * \f$ J = \left ( \begin{array}{cc} c & \overline s \\ -s & \overline c \end{array} \right ) \f$ * - * \sa MatrixBase::makeJacobi(), MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight() + * \sa MatrixBase::makeJacobi(), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ -template class JacobiRotation +template class PlanarRotation { public: - /** Default constructor without any initialization. */ - JacobiRotation() {} + typedef typename NumTraits::Real RealScalar; - /** Construct a Jacobi rotation from a cosine-sine pair (\a c, \c s). */ - JacobiRotation(const Scalar& c, const Scalar& s) : m_c(c), m_s(s) {} + /** Default constructor without any initialization. */ + PlanarRotation() {} + + /** Construct a planar rotation from a cosine-sine pair (\a c, \c s). */ + PlanarRotation(const Scalar& c, const Scalar& s) : m_c(c), m_s(s) {} Scalar& c() { return m_c; } Scalar c() const { return m_c; } Scalar& s() { return m_s; } Scalar s() const { return m_s; } - /** Concatenates two Jacobi rotation */ - JacobiRotation operator*(const JacobiRotation& other) + /** Concatenates two planar rotation */ + PlanarRotation operator*(const PlanarRotation& other) { - return JacobiRotation(m_c * other.m_c - ei_conj(m_s) * other.m_s, + return PlanarRotation(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))); } /** Returns the transposed transformation */ - JacobiRotation transpose() const { return JacobiRotation(m_c, -ei_conj(m_s)); } + PlanarRotation transpose() const { return PlanarRotation(m_c, -ei_conj(m_s)); } /** Returns the adjoint transformation */ - JacobiRotation adjoint() const { return JacobiRotation(ei_conj(m_c), -m_s); } + PlanarRotation adjoint() const { return PlanarRotation(ei_conj(m_c), -m_s); } + + template + bool makeJacobi(const MatrixBase&, int p, int q); + bool makeJacobi(RealScalar x, Scalar y, RealScalar z); + + 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); + Scalar m_c, m_s; }; -/** Applies the clock wise 2D rotation \a j to the set of 2D vectors of cordinates \a x and \a y: - * \f$ \left ( \begin{array}{cc} x \\ y \end{array} \right ) = J \left ( \begin{array}{cc} x \\ y \end{array} \right ) \f$ - * - * \sa MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight() - */ -template -void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation& j); - -/** 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::applyJacobiOnTheRight(), ei_apply_rotation_in_the_plane() - */ -template -template -inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, const JacobiRotation& j) -{ - RowXpr x(row(p)); - RowXpr y(row(q)); - ei_apply_rotation_in_the_plane(x, y, j); -} - -/** 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::applyJacobiOnTheLeft(), ei_apply_rotation_in_the_plane() - */ -template -template -inline void MatrixBase::applyJacobiOnTheRight(int p, int q, const JacobiRotation& j) -{ - ColXpr x(col(p)); - ColXpr y(col(q)); - ei_apply_rotation_in_the_plane(x, y, j.transpose()); -} - -/** Computes the Jacobi rotation \a J such that applying \a J on both the right and left sides of the 2x2 matrix +/** Makes \c *this as a Jacobi rotation \a J such that applying \a J on both the right and left sides of the 2x2 matrix * \f$ B = \left ( \begin{array}{cc} x & y \\ * & z \end{array} \right )\f$ yields * a diagonal matrix \f$ A = J^* B J \f$ * - * \sa MatrixBase::makeJacobi(), MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight() + * \sa MatrixBase::makeJacobi(), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ template -bool ei_makeJacobi(typename NumTraits::Real x, Scalar y, typename NumTraits::Real z, JacobiRotation *j) +bool PlanarRotation::makeJacobi(RealScalar x, Scalar y, RealScalar z) { typedef typename NumTraits::Real RealScalar; if(y == Scalar(0)) { - j->c() = Scalar(1); - j->s() = Scalar(0); + m_c = Scalar(1); + m_s = Scalar(0); return false; } else @@ -135,26 +110,132 @@ bool ei_makeJacobi(typename NumTraits::Real x, Scalar y, typename NumTra } RealScalar sign_t = t > 0 ? 1 : -1; RealScalar n = RealScalar(1) / ei_sqrt(ei_abs2(t)+1); - j->s() = - sign_t * (ei_conj(y) / ei_abs(y)) * ei_abs(t) * n; - j->c() = n; + m_s = - sign_t * (ei_conj(y) / ei_abs(y)) * ei_abs(t) * n; + m_c = n; return true; } } -/** Computes the Jacobi rotation \a J such that applying \a J on both the right and left sides of the 2x2 matrix +/** Makes \c *this as a Jacobi rotation \c J such that applying \a J on both the right and left sides of the 2x2 matrix * \f$ B = \left ( \begin{array}{cc} \text{this}_{pp} & \text{this}_{pq} \\ * & \text{this}_{qq} \end{array} \right )\f$ yields * a diagonal matrix \f$ A = J^* B J \f$ * - * \sa MatrixBase::ei_make_jacobi(), MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight() + * \sa PlanarRotation::makeJacobi(RealScalar, Scalar, RealScalar), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ +template template -inline bool MatrixBase::makeJacobi(int p, int q, JacobiRotation *j) const +inline bool PlanarRotation::makeJacobi(const MatrixBase& m, int p, int q) { - return ei_makeJacobi(ei_real(coeff(p,p)), coeff(p,q), ei_real(coeff(q,q)), j); + return makeJacobi(ei_real(m.coeff(p,p)), m.coeff(p,q), ei_real(m.coeff(q,q))); } -template -void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation& j) +/** Makes \c *this as a Givens rotation \c G such that applying \f$ G^* \f$ to the left of the vector + * \f$ V = \left ( \begin{array}{c} p \\ q \end{array} \right )\f$ yields: + * \f$ G^* V = \left ( \begin{array}{c} z \\ 0 \end{array} \right )\f$. + * + * The value of \a z is returned if \a z is not null (the default is null). + * Also note that G is built such that the cosine is always real. + * + * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() + */ +template +void PlanarRotation::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()); +} + + +// specialization for complexes +template +void PlanarRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar* z, ei_meta_true) +{ + RealScalar scale, absx, absxy; + if(q==Scalar(0)) + { + // return identity + m_c = Scalar(1); + m_s = Scalar(0); + if(z) *z = p; + } + else + { + scale = ei_norm1(p); + absx = scale * ei_sqrt(ei_abs2(p/scale)); + scale = ei_abs(scale) + ei_norm1(q); + absxy = scale * ei_sqrt((absx/scale)*(absx/scale) + ei_abs2(q/scale)); + m_c = Scalar(absx / absxy); + Scalar np = p/absx; + m_s = -ei_conj(np) * q / absxy; + if(z) *z = np * absxy; + } +} + +// specialization for reals +template +void PlanarRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar* z, ei_meta_false) +{ + // from Golub's "Matrix Computations", algorithm 5.1.3 + if(q==0) + { + m_c = 1; m_s = 0; + } + else if(ei_abs(q)>ei_abs(p)) + { + Scalar t = -p/q; + m_s = Scalar(1)/ei_sqrt(1+t*t); + m_c = m_s * t; + } + else + { + Scalar t = -q/p; + m_c = Scalar(1)/ei_sqrt(1+t*t); + m_s = m_c * t; + } +} + +/**************************************************************************************** +* Implementation of MatrixBase methods +/***************************************************************************************/ + +/** Applies the clock wise 2D rotation \a j to the set of 2D vectors of cordinates \a x and \a y: + * \f$ \left ( \begin{array}{cc} x \\ y \end{array} \right ) = J \left ( \begin{array}{cc} x \\ y \end{array} \right ) \f$ + * + * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() + */ +template +void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const PlanarRotation& j); + +/** 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 PlanarRotation, MatrixBase::applyOnTheRight(), ei_apply_rotation_in_the_plane() + */ +template +template +inline void MatrixBase::applyOnTheLeft(int p, int q, const PlanarRotation& j) +{ + RowXpr x(row(p)); + RowXpr y(row(q)); + ei_apply_rotation_in_the_plane(x, y, j); +} + +/** 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 PlanarRotation, MatrixBase::applyOnTheLeft(), ei_apply_rotation_in_the_plane() + */ +template +template +inline void MatrixBase::applyOnTheRight(int p, int q, const PlanarRotation& j) +{ + ColXpr x(col(p)); + ColXpr y(col(q)); + ei_apply_rotation_in_the_plane(x, y, j.transpose()); +} + + +template +void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const PlanarRotation& j) { typedef typename VectorX::Scalar Scalar; ei_assert(_x.size() == _y.size()); diff --git a/Eigen/src/QR/ComplexSchur.h b/Eigen/src/QR/ComplexSchur.h index a0b954d49..153826811 100644 --- a/Eigen/src/QR/ComplexSchur.h +++ b/Eigen/src/QR/ComplexSchur.h @@ -80,34 +80,6 @@ template class ComplexSchur bool m_isInitialized; }; -// computes the plane rotation G such that G' x |p| = | c s' |' |p| = |z| -// |q| |-s c' | |q| |0| -// and returns z if requested. Note that the returned c is real. -template void ei_make_givens(const std::complex& p, const std::complex& q, - JacobiRotation >& rot, std::complex* z=0) -{ - typedef std::complex Complex; - T scale, absx, absxy; - if(p==Complex(0)) - { - // return identity - rot.c() = Complex(1,0); - rot.s() = Complex(0,0); - if(z) *z = p; - } - else - { - scale = cnorm1(p); - absx = scale * ei_sqrt(ei_abs2(p/scale)); - scale = ei_abs(scale) + cnorm1(q); - absxy = scale * ei_sqrt((absx/scale)*(absx/scale) + ei_abs2(q/scale)); - rot.c() = Complex(absx / absxy); - Complex np = p/absx; - rot.s() = -ei_conj(np) * q / absxy; - if(z) *z = np * absxy; - } -} - template void ComplexSchur::compute(const MatrixType& matrix) { @@ -133,8 +105,8 @@ void ComplexSchur::compute(const MatrixType& matrix) //locate the range in which to iterate while(iu > 0) { - d = cnorm1(m_matT.coeffRef(iu,iu)) + cnorm1(m_matT.coeffRef(iu-1,iu-1)); - sd = cnorm1(m_matT.coeffRef(iu,iu-1)); + d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1)); + sd = ei_norm1(m_matT.coeffRef(iu,iu-1)); if(sd >= eps * d) break; // FIXME : precision criterion ?? @@ -156,8 +128,8 @@ void ComplexSchur::compute(const MatrixType& matrix) while( il > 0 ) { // check if the current 2x2 block on the diagonal is upper triangular - d = cnorm1(m_matT.coeffRef(il,il)) + cnorm1(m_matT.coeffRef(il-1,il-1)); - sd = cnorm1(m_matT.coeffRef(il,il-1)); + d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1)); + sd = ei_norm1(m_matT.coeffRef(il,il-1)); if(sd < eps * d) break; // FIXME : precision criterion ?? @@ -179,32 +151,32 @@ void ComplexSchur::compute(const MatrixType& matrix) r1 = (b+disc)/RealScalar(2); r2 = (b-disc)/RealScalar(2); - if(cnorm1(r1) > cnorm1(r2)) + if(ei_norm1(r1) > ei_norm1(r2)) r2 = c/r1; else r1 = c/r2; - if(cnorm1(r1-t.coeff(1,1)) < cnorm1(r2-t.coeff(1,1))) + if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1))) kappa = sf * r1; else kappa = sf * r2; // perform the QR step using Givens rotations - JacobiRotation rot; - ei_make_givens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il), rot); + PlanarRotation rot; + rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il)); for(int i=il ; i::compute(const MatrixType& matrix) m_isInitialized = true; } -// norm1 of complex numbers -template -T cnorm1(const std::complex &Z) -{ - return(ei_abs(Z.real()) + ei_abs(Z.imag())); -} - /** * Computes the principal value of the square root of the complex \a z. */ diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index dd1736e28..580b042f6 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -135,28 +135,6 @@ template class SelfAdjointEigenSolver #ifndef EIGEN_HIDE_HEAVY_CODE -// from Golub's "Matrix Computations", algorithm 5.1.3 -template -static void ei_givens_rotation(Scalar a, Scalar b, Scalar& c, Scalar& s) -{ - if (b==0) - { - c = 1; s = 0; - } - else if (ei_abs(b)>ei_abs(a)) - { - Scalar t = -a/b; - s = Scalar(1)/ei_sqrt(1+t*t); - c = s * t; - } - else - { - Scalar t = -b/a; - c = Scalar(1)/ei_sqrt(1+t*t); - s = c * t; - } -} - /** \internal * * \qr_module @@ -353,34 +331,33 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st for (int k = start; k < end; ++k) { - RealScalar c, s; - ei_givens_rotation(x, z, c, s); + PlanarRotation rot; + rot.makeGivens(x, z); // do T = G' T G - RealScalar sdk = s * diag[k] + c * subdiag[k]; - RealScalar dkp1 = s * subdiag[k] + c * diag[k+1]; + RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k]; + RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1]; - diag[k] = c * (c * diag[k] - s * subdiag[k]) - s * (c * subdiag[k] - s * diag[k+1]); - diag[k+1] = s * sdk + c * dkp1; - subdiag[k] = c * sdk - s * dkp1; + diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]); + diag[k+1] = rot.s() * sdk + rot.c() * dkp1; + subdiag[k] = rot.c() * sdk - rot.s() * dkp1; if (k > start) - subdiag[k - 1] = c * subdiag[k-1] - s * z; + subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; x = subdiag[k]; if (k < end - 1) { - z = -s * subdiag[k+1]; - subdiag[k + 1] = c * subdiag[k+1]; + z = -rot.s() * subdiag[k+1]; + subdiag[k + 1] = rot.c() * subdiag[k+1]; } // apply the givens rotation to the unit matrix Q = Q * G - // G only modifies the two columns k and k+1 if (matrixQ) { Map > q(matrixQ,n,n); - q.applyJacobiOnTheRight(k,k+1,JacobiRotation(c,s)); + q.applyOnTheRight(k,k+1,rot); } } } diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 654f8981a..ca359832b 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -125,7 +125,7 @@ struct ei_svd_precondition_2x2_block_to_be_real static void run(MatrixType& work_matrix, JacobiSVD& svd, int p, int q) { Scalar z; - JacobiRotation rot; + PlanarRotation rot; RealScalar n = ei_sqrt(ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(q,p))); if(n==0) { @@ -140,8 +140,8 @@ struct ei_svd_precondition_2x2_block_to_be_real { rot.c() = ei_conj(work_matrix.coeff(p,p)) / n; rot.s() = work_matrix.coeff(q,p) / n; - work_matrix.applyJacobiOnTheLeft(p,q,rot); - if(ComputeU) svd.m_matrixU.applyJacobiOnTheRight(p,q,rot.adjoint()); + work_matrix.applyOnTheLeft(p,q,rot); + if(ComputeU) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); if(work_matrix.coeff(p,q) != Scalar(0)) { Scalar z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); @@ -160,13 +160,13 @@ struct ei_svd_precondition_2x2_block_to_be_real template void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q, - JacobiRotation *j_left, - JacobiRotation *j_right) + PlanarRotation *j_left, + PlanarRotation *j_right) { Matrix m; m << ei_real(matrix.coeff(p,p)), ei_real(matrix.coeff(p,q)), ei_real(matrix.coeff(q,p)), ei_real(matrix.coeff(q,q)); - JacobiRotation rot1; + PlanarRotation rot1; RealScalar t = m.coeff(0,0) + m.coeff(1,1); RealScalar d = m.coeff(1,0) - m.coeff(0,1); if(t == RealScalar(0)) @@ -180,8 +180,8 @@ void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q, rot1.c() = RealScalar(1) / ei_sqrt(1 + ei_abs2(u)); rot1.s() = rot1.c() * u; } - m.applyJacobiOnTheLeft(0,1,rot1); - m.makeJacobi(0,1,j_right); + m.applyOnTheLeft(0,1,rot1); + j_right->makeJacobi(m,0,1); *j_left = rot1 * j_right->transpose(); } @@ -214,7 +214,7 @@ JacobiSVD& JacobiSVD::compute(const Ma if(ComputeU) for(int i = 0; i < rows; i++) m_matrixU.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); - + } else { @@ -232,14 +232,14 @@ sweep_again: { ei_svd_precondition_2x2_block_to_be_real::run(work_matrix, *this, p, q); - JacobiRotation j_left, j_right; + PlanarRotation j_left, j_right; ei_real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right); - work_matrix.applyJacobiOnTheLeft(p,q,j_left); - if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,j_left.transpose()); + work_matrix.applyOnTheLeft(p,q,j_left); + if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); - work_matrix.applyJacobiOnTheRight(p,q,j_right); - if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,j_right); + work_matrix.applyOnTheRight(p,q,j_right); + if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right); } } } diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 095e18b3e..da01cf396 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -309,7 +309,7 @@ SVD& SVD::compute(const MatrixType& matrix) h = Scalar(1.0)/h; c = g*h; s = -f*h; - V.applyJacobiOnTheRight(i,nm,JacobiRotation(c,s)); + V.applyOnTheRight(i,nm,PlanarRotation(c,s)); } } z = W[k]; @@ -342,7 +342,7 @@ SVD& SVD::compute(const MatrixType& matrix) y = W[i]; h = s*g; g = c*g; - + z = pythag(f,h); rv1[j] = z; c = f/z; @@ -351,8 +351,8 @@ SVD& SVD::compute(const MatrixType& matrix) g = g*c - x*s; h = y*s; y *= c; - V.applyJacobiOnTheRight(i,j,JacobiRotation(c,s)); - + V.applyOnTheRight(i,j,PlanarRotation(c,s)); + z = pythag(f,h); W[j] = z; // Rotation can be arbitrary if z = 0. @@ -364,7 +364,7 @@ SVD& SVD::compute(const MatrixType& matrix) } f = c*g + s*y; x = c*y - s*g; - A.applyJacobiOnTheRight(i,j,JacobiRotation(c,s)); + A.applyOnTheRight(i,j,PlanarRotation(c,s)); } rv1[l] = 0.0; rv1[k] = f;