diff options
| -rw-r--r-- | Eigen/SVD | 2 | ||||
| -rw-r--r-- | Eigen/src/Core/MatrixBase.h | 8 | ||||
| -rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 1 | ||||
| -rw-r--r-- | Eigen/src/Jacobi/Jacobi.h | 113 | ||||
| -rw-r--r-- | Eigen/src/SVD/JacobiSVD.h | 258 | ||||
| -rw-r--r-- | Eigen/src/SVD/JacobiSquareSVD.h | 169 | ||||
| -rw-r--r-- | test/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | test/jacobisvd.cpp | 105 |
8 files changed, 416 insertions, 241 deletions
@@ -23,7 +23,7 @@ namespace Eigen { */ #include "src/SVD/SVD.h" -#include "src/SVD/JacobiSquareSVD.h" +#include "src/SVD/JacobiSVD.h" } // namespace Eigen diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 7cb8853e6..9b1bf9e19 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -803,11 +803,11 @@ template<typename Derived> class MatrixBase ///////// Jacobi module ///////// - void applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s); - void applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s); + template<typename JacobiScalar> + void applyJacobiOnTheLeft(int p, int q, JacobiScalar c, JacobiScalar s); + template<typename JacobiScalar> + void applyJacobiOnTheRight(int p, int q, JacobiScalar c, JacobiScalar s); bool makeJacobi(int p, int q, Scalar *c, Scalar *s) const; - bool makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const; - bool makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const; #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 414aaceca..d7dc61e73 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -120,6 +120,7 @@ template<typename MatrixType> class HouseholderQR; template<typename MatrixType> class ColPivotingHouseholderQR; template<typename MatrixType> class FullPivotingHouseholderQR; template<typename MatrixType> class SVD; +template<typename MatrixType, unsigned int Options = 0> class JacobiSVD; template<typename MatrixType, int UpLo = LowerTriangular> class LLT; template<typename MatrixType> class LDLT; 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<typename VectorX, typename VectorY> -void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s); +template<typename VectorX, typename VectorY, typename JacobiScalar> +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<typename Derived> -inline void MatrixBase<Derived>::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) +template<typename JacobiScalar> +inline void MatrixBase<Derived>::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<Derived>::applyJacobiOnTheLeft(int p, int q, Scalar c, Sc * \sa MatrixBase::applyJacobiOnTheLeft(), ei_apply_rotation_in_the_plane() */ template<typename Derived> -inline void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) +template<typename JacobiScalar> +inline void MatrixBase<Derived>::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<typename Scalar> -bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) +bool ei_makeJacobi(typename NumTraits<Scalar>::Real x, Scalar y, typename NumTraits<Scalar>::Real z, Scalar *c, Scalar *s) { - if(y == 0) + typedef typename NumTraits<Scalar>::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<typename Derived> inline bool MatrixBase<Derived>::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<typename Derived> -inline bool MatrixBase<Derived>::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<typename Derived> -inline bool MatrixBase<Derived>::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<typename Scalar> -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<typename VectorX, typename VectorY> -void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s) +template<typename VectorX, typename VectorY, typename JacobiScalar> +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<Scalar>::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<NumTraits<Scalar>::IsComplex,false> cj; for(int i=0; i<alignedStart; ++i) { Scalar xi = x[i]; Scalar yi = y[i]; - x[i] = c * xi - ei_conj(s) * yi; - y[i] = s * xi + c * yi; + x[i] = c * xi + ei_conj(s) * yi; + y[i] = - s * xi + ei_conj(c) * yi; } Scalar* px = x + alignedStart; @@ -168,8 +147,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Packet xi = ei_pload(px); Packet yi = ei_pload(py); - ei_pstore(px, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi))); - ei_pstore(py, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi))); + ei_pstore(px, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi))); + ei_pstore(py, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi))); px += PacketSize; py += PacketSize; } @@ -183,10 +162,10 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& Packet xi1 = ei_ploadu(px+PacketSize); Packet yi = ei_pload (py); Packet yi1 = ei_pload (py+PacketSize); - ei_pstoreu(px, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi))); - ei_pstoreu(px+PacketSize, ei_psub(ei_pmul(pc,xi1),cj.pmul(ps,yi1))); - ei_pstore (py, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi))); - ei_pstore (py+PacketSize, ei_padd(ei_pmul(ps,xi1),ei_pmul(pc,yi1))); + ei_pstoreu(px, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi))); + ei_pstoreu(px+PacketSize, ei_padd(ei_pmul(pc,xi1),cj.pmul(ps,yi1))); + ei_pstore (py, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi))); + ei_pstore (py+PacketSize, ei_psub(ei_pmul(pc,yi1),ei_pmul(ps,xi1))); px += Peeling*PacketSize; py += Peeling*PacketSize; } @@ -194,8 +173,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Packet xi = ei_ploadu(x+peelingEnd); Packet yi = ei_pload (y+peelingEnd); - ei_pstoreu(x+peelingEnd, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi))); - ei_pstore (y+peelingEnd, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi))); + ei_pstoreu(x+peelingEnd, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi))); + ei_pstore (y+peelingEnd, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi))); } } @@ -203,8 +182,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Scalar xi = x[i]; Scalar yi = y[i]; - x[i] = c * xi - ei_conj(s) * yi; - y[i] = s * xi + c * yi; + x[i] = c * xi + ei_conj(s) * yi; + y[i] = -s * xi + ei_conj(c) * yi; } } else @@ -213,8 +192,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Scalar xi = *x; Scalar yi = *y; - *x = c * xi - ei_conj(s) * yi; - *y = s * xi + c * yi; + *x = c * xi + ei_conj(s) * yi; + *y = -s * xi + ei_conj(c) * yi; x += incrx; y += incry; } diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h new file mode 100644 index 000000000..e5e7fb6fc --- /dev/null +++ b/Eigen/src/SVD/JacobiSVD.h @@ -0,0 +1,258 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_JACOBISVD_H +#define EIGEN_JACOBISVD_H + +/** \ingroup SVD_Module + * \nonstableyet + * + * \class JacobiSVD + * + * \brief Jacobi SVD decomposition of a square matrix + * + * \param MatrixType the type of the matrix of which we are computing the SVD decomposition + * \param ComputeU whether the U matrix should be computed + * \param ComputeV whether the V matrix should be computed + * + * \sa MatrixBase::jacobiSvd() + */ +template<typename MatrixType, unsigned int Options> class JacobiSVD +{ + private: + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; + enum { + ComputeU = 1, + ComputeV = 1, + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + DiagSizeAtCompileTime = EIGEN_ENUM_MIN(RowsAtCompileTime,ColsAtCompileTime), + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + MaxDiagSizeAtCompileTime = EIGEN_ENUM_MIN(MaxRowsAtCompileTime,MaxColsAtCompileTime), + MatrixOptions = MatrixType::Options + }; + + typedef Matrix<Scalar, Dynamic, Dynamic, MatrixOptions> DummyMatrixType; + typedef typename ei_meta_if<ComputeU, + Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, + MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>, + DummyMatrixType>::ret MatrixUType; + typedef typename ei_meta_if<ComputeV, + Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, + MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>, + DummyMatrixType>::ret MatrixVType; + typedef Matrix<RealScalar, DiagSizeAtCompileTime, 1, + Options, MaxDiagSizeAtCompileTime, 1> SingularValuesType; + typedef Matrix<Scalar, 1, RowsAtCompileTime, MatrixOptions, 1, MaxRowsAtCompileTime> RowType; + typedef Matrix<Scalar, RowsAtCompileTime, 1, MatrixOptions, MaxRowsAtCompileTime, 1> ColType; + + public: + + JacobiSVD() : m_isInitialized(false) {} + + JacobiSVD(const MatrixType& matrix) : m_isInitialized(false) + { + compute(matrix); + } + + JacobiSVD& compute(const MatrixType& matrix); + + const MatrixUType& matrixU() const + { + ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + return m_matrixU; + } + + const SingularValuesType& singularValues() const + { + ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + return m_singularValues; + } + + const MatrixUType& matrixV() const + { + ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + return m_matrixV; + } + + protected: + MatrixUType m_matrixU; + MatrixVType m_matrixV; + SingularValuesType m_singularValues; + bool m_isInitialized; + + template<typename _MatrixType, unsigned int _Options, bool _IsComplex> + friend struct ei_svd_precondition_2x2_block_to_be_real; +}; + +template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> +struct ei_svd_precondition_2x2_block_to_be_real +{ + static void run(MatrixType&, JacobiSVD<MatrixType, Options>&, int, int) {} +}; + +template<typename MatrixType, unsigned int Options> +struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options, true> +{ + typedef JacobiSVD<MatrixType, Options> SVD; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + + enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV }; + static void run(MatrixType& work_matrix, JacobiSVD<MatrixType, Options>& svd, int p, int q) + { + Scalar c, s, z; + RealScalar n = ei_sqrt(ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(q,p))); + if(n==0) + { + z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); + work_matrix.row(p) *= z; + if(ComputeU) svd.m_matrixU.col(p) *= ei_conj(z); + z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + work_matrix.row(q) *= z; + if(ComputeU) svd.m_matrixU.col(q) *= ei_conj(z); + } + else + { + c = ei_conj(work_matrix.coeff(p,p)) / n; + s = work_matrix.coeff(q,p) / n; + work_matrix.applyJacobiOnTheLeft(p,q,c,s); + if(ComputeU) svd.m_matrixU.applyJacobiOnTheRight(p,q,ei_conj(c),-s); + if(work_matrix.coeff(p,q) != Scalar(0)) + { + Scalar z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); + work_matrix.col(q) *= z; + if(ComputeV) svd.m_matrixV.col(q) *= z; + } + if(work_matrix.coeff(q,q) != Scalar(0)) + { + z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + work_matrix.row(q) *= z; + if(ComputeU) svd.m_matrixU.col(q) *= ei_conj(z); + } + } + } +}; + +template<typename MatrixType, typename RealScalar> +void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q, + RealScalar *c_left, RealScalar *s_left, + RealScalar *c_right, RealScalar *s_right) +{ + Matrix<RealScalar,2,2> 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)); + RealScalar c1, s1; + RealScalar t = m.coeff(0,0) + m.coeff(1,1); + RealScalar d = m.coeff(1,0) - m.coeff(0,1); + if(t == RealScalar(0)) + { + c1 = 0; + s1 = d > 0 ? 1 : -1; + } + else + { + RealScalar u = d / t; + c1 = RealScalar(1) / ei_sqrt(1 + ei_abs2(u)); + s1 = c1 * u; + } + m.applyJacobiOnTheLeft(0,1,c1,s1); + RealScalar c2, s2; + m.makeJacobi(0,1,&c2,&s2); + *c_left = c1*c2 + s1*s2; + *s_left = s1*c2 - c1*s2; + *c_right = c2; + *s_right = s2; +} + +template<typename MatrixType, unsigned int Options> +JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const MatrixType& matrix) +{ + MatrixType work_matrix(matrix); + int size = matrix.rows(); + if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); + if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); + m_singularValues.resize(size); + const RealScalar precision = 2 * epsilon<Scalar>(); + +sweep_again: + for(int p = 1; p < size; ++p) + { + for(int q = 0; q < p; ++q) + { + if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) + > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision) + { + ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(work_matrix, *this, p, q); + + RealScalar c_left, s_left, c_right, s_right; + ei_real_2x2_jacobi_svd(work_matrix, p, q, &c_left, &s_left, &c_right, &s_right); + + work_matrix.applyJacobiOnTheLeft(p,q,c_left,s_left); + if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c_left,-s_left); + + work_matrix.applyJacobiOnTheRight(p,q,c_right,s_right); + if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c_right,s_right); + } + } + } + + RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff(); + RealScalar maxAllowedOffDiag = biggestOnDiag * precision; + for(int p = 0; p < size; ++p) + { + for(int q = 0; q < p; ++q) + if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) + goto sweep_again; + for(int q = p+1; q < size; ++q) + if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) + goto sweep_again; + } + + for(int i = 0; i < size; ++i) + { + RealScalar a = ei_abs(work_matrix.coeff(i,i)); + m_singularValues.coeffRef(i) = a; + if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a; + } + + for(int i = 0; i < size; i++) + { + int pos; + m_singularValues.end(size-i).maxCoeff(&pos); + if(pos) + { + pos += i; + std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); + if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i)); + if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); + } + } + + m_isInitialized = true; + return *this; +} +#endif // EIGEN_JACOBISVD_H diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h deleted file mode 100644 index f21c9edca..000000000 --- a/Eigen/src/SVD/JacobiSquareSVD.h +++ /dev/null @@ -1,169 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_JACOBISQUARESVD_H -#define EIGEN_JACOBISQUARESVD_H - -/** \ingroup SVD_Module - * \nonstableyet - * - * \class JacobiSquareSVD - * - * \brief Jacobi SVD decomposition of a square matrix - * - * \param MatrixType the type of the matrix of which we are computing the SVD decomposition - * \param ComputeU whether the U matrix should be computed - * \param ComputeV whether the V matrix should be computed - * - * \sa MatrixBase::jacobiSvd() - */ -template<typename MatrixType, bool ComputeU, bool ComputeV> class JacobiSquareSVD -{ - private: - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - Options = MatrixType::Options - }; - - typedef Matrix<Scalar, Dynamic, Dynamic, Options> DummyMatrixType; - typedef typename ei_meta_if<ComputeU, - Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, - Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime>, - DummyMatrixType>::ret MatrixUType; - typedef typename Diagonal<MatrixType,0>::PlainMatrixType SingularValuesType; - typedef Matrix<Scalar, 1, RowsAtCompileTime, Options, 1, MaxRowsAtCompileTime> RowType; - typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> ColType; - - public: - - JacobiSquareSVD() : m_isInitialized(false) {} - - JacobiSquareSVD(const MatrixType& matrix) : m_isInitialized(false) - { - compute(matrix); - } - - JacobiSquareSVD& compute(const MatrixType& matrix); - - const MatrixUType& matrixU() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matrixU; - } - - const SingularValuesType& singularValues() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_singularValues; - } - - const MatrixUType& matrixV() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matrixV; - } - - protected: - MatrixUType m_matrixU; - MatrixUType m_matrixV; - SingularValuesType m_singularValues; - bool m_isInitialized; -}; - -template<typename MatrixType, bool ComputeU, bool ComputeV> -JacobiSquareSVD<MatrixType, ComputeU, ComputeV>& JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType& matrix) -{ - MatrixType work_matrix(matrix); - int size = matrix.rows(); - if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); - if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); - m_singularValues.resize(size); - const RealScalar precision = 2 * epsilon<Scalar>(); - -sweep_again: - for(int p = 1; p < size; ++p) - { - for(int q = 0; q < p; ++q) - { - Scalar c, s; - while(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) - > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision) - { - if(work_matrix.makeJacobiForAtA(p,q,&c,&s)) - { - work_matrix.applyJacobiOnTheRight(p,q,c,s); - if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s); - } - if(work_matrix.makeJacobiForAAt(p,q,&c,&s)) - { - ei_normalizeJacobi(&c, &s, work_matrix.coeff(p,p), work_matrix.coeff(q,p)), - work_matrix.applyJacobiOnTheLeft(p,q,c,s); - if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c,s); - } - } - } - } - - RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff(); - RealScalar maxAllowedOffDiag = biggestOnDiag * precision; - for(int p = 0; p < size; ++p) - { - for(int q = 0; q < p; ++q) - if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) - goto sweep_again; - for(int q = p+1; q < size; ++q) - if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) - goto sweep_again; - } - - m_singularValues = work_matrix.diagonal().cwise().abs(); - RealScalar biggestSingularValue = m_singularValues.maxCoeff(); - - for(int i = 0; i < size; ++i) - { - RealScalar a = ei_abs(work_matrix.coeff(i,i)); - m_singularValues.coeffRef(i) = a; - if(ComputeU && !ei_isMuchSmallerThan(a, biggestSingularValue)) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a; - } - - for(int i = 0; i < size; i++) - { - int pos; - m_singularValues.end(size-i).maxCoeff(&pos); - if(pos) - { - pos += i; - std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); - if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i)); - if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); - } - } - - m_isInitialized = true; - return *this; -} -#endif // EIGEN_JACOBISQUARESVD_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fc43bbb1d..896392188 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -126,6 +126,7 @@ ei_add_test(qr_fullpivoting) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") ei_add_test(svd) +ei_add_test(jacobisvd ${EI_OFLAG}) ei_add_test(geo_orthomethods) ei_add_test(geo_homogeneous) ei_add_test(geo_quaternion) diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp new file mode 100644 index 000000000..9a4d79e45 --- /dev/null +++ b/test/jacobisvd.cpp @@ -0,0 +1,105 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/SVD> +#include <Eigen/LU> + +template<typename MatrixType> void svd(const MatrixType& m, bool pickrandom = true) +{ + int rows = m.rows(); + int cols = m.cols(); + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; + typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; + typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; + typedef Matrix<Scalar, ColsAtCompileTime, 1> InputVectorType; + + MatrixType a; + if(pickrandom) a = MatrixType::Random(rows,cols); + else a = m; + + JacobiSVD<MatrixType> svd(a); + MatrixType sigma = MatrixType::Zero(rows,cols); + sigma.diagonal() = svd.singularValues().template cast<Scalar>(); + MatrixUType u = svd.matrixU(); + MatrixVType v = svd.matrixV(); + + VERIFY_IS_APPROX(a, u * sigma * v.adjoint()); + VERIFY_IS_UNITARY(u); + VERIFY_IS_UNITARY(v); +} + +template<typename MatrixType> void svd_verify_assert() +{ + MatrixType tmp; + + SVD<MatrixType> svd; + //VERIFY_RAISES_ASSERT(svd.solve(tmp, &tmp)) + VERIFY_RAISES_ASSERT(svd.matrixU()) + VERIFY_RAISES_ASSERT(svd.singularValues()) + VERIFY_RAISES_ASSERT(svd.matrixV()) + /*VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp))*/ +} + +void test_jacobisvd() +{ + for(int i = 0; i < g_repeat; i++) { + Matrix2cd m; + m << 0, 1, + 0, 1; + CALL_SUBTEST( svd(m, false) ); + m << 1, 0, + 1, 0; + CALL_SUBTEST( svd(m, false) ); + Matrix2d n; + n << 1, 1, + 1, -1; + CALL_SUBTEST( svd(n, false) ); + CALL_SUBTEST( svd(Matrix3f()) ); + CALL_SUBTEST( svd(Matrix4d()) ); + CALL_SUBTEST( svd(MatrixXf(50,50)) ); +// CALL_SUBTEST( svd(MatrixXd(14,7)) ); + CALL_SUBTEST( svd(MatrixXcf(3,3)) ); + CALL_SUBTEST( svd(MatrixXd(30,30)) ); + } + CALL_SUBTEST( svd(MatrixXf(200,200)) ); + CALL_SUBTEST( svd(MatrixXcd(100,100)) ); + + CALL_SUBTEST( svd_verify_assert<Matrix3f>() ); + CALL_SUBTEST( svd_verify_assert<Matrix3d>() ); + CALL_SUBTEST( svd_verify_assert<MatrixXf>() ); + CALL_SUBTEST( svd_verify_assert<MatrixXd>() ); +} |