diff options
Diffstat (limited to 'eigen/unsupported/Eigen/src/SVD/BDCSVD.h')
| -rw-r--r-- | eigen/unsupported/Eigen/src/SVD/BDCSVD.h | 748 |
1 files changed, 0 insertions, 748 deletions
diff --git a/eigen/unsupported/Eigen/src/SVD/BDCSVD.h b/eigen/unsupported/Eigen/src/SVD/BDCSVD.h deleted file mode 100644 index 11d4882..0000000 --- a/eigen/unsupported/Eigen/src/SVD/BDCSVD.h +++ /dev/null @@ -1,748 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// We used the "A Divide-And-Conquer Algorithm for the Bidiagonal SVD" -// research report written by Ming Gu and Stanley C.Eisenstat -// The code variable names correspond to the names they used in their -// report -// -// Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com> -// Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr> -// Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr> -// Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr> -// -// Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#ifndef EIGEN_BDCSVD_H -#define EIGEN_BDCSVD_H - -#define EPSILON 0.0000000000000001 - -#define ALGOSWAP 32 - -namespace Eigen { -/** \ingroup SVD_Module - * - * - * \class BDCSVD - * - * \brief class Bidiagonal Divide and Conquer SVD - * - * \param MatrixType the type of the matrix of which we are computing the SVD decomposition - * We plan to have a very similar interface to JacobiSVD on this class. - * It should be used to speed up the calcul of SVD for big matrices. - */ -template<typename _MatrixType> -class BDCSVD : public SVDBase<_MatrixType> -{ - typedef SVDBase<_MatrixType> Base; - -public: - using Base::rows; - using Base::cols; - - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - typedef typename MatrixType::Index Index; - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime, - DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime), - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime, MaxColsAtCompileTime), - MatrixOptions = MatrixType::Options - }; - - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, - MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> - MatrixUType; - typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, - MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> - MatrixVType; - typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; - typedef typename internal::plain_row_type<MatrixType>::type RowType; - typedef typename internal::plain_col_type<MatrixType>::type ColType; - typedef Matrix<Scalar, Dynamic, Dynamic> MatrixX; - typedef Matrix<RealScalar, Dynamic, Dynamic> MatrixXr; - typedef Matrix<RealScalar, Dynamic, 1> VectorType; - - /** \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via BDCSVD::compute(const MatrixType&). - */ - BDCSVD() - : SVDBase<_MatrixType>::SVDBase(), - algoswap(ALGOSWAP) - {} - - - /** \brief Default Constructor with memory preallocation - * - * Like the default constructor but with preallocation of the internal data - * according to the specified problem size. - * \sa BDCSVD() - */ - BDCSVD(Index rows, Index cols, unsigned int computationOptions = 0) - : SVDBase<_MatrixType>::SVDBase(), - algoswap(ALGOSWAP) - { - allocate(rows, cols, computationOptions); - } - - /** \brief Constructor performing the decomposition of given matrix. - * - * \param matrix the matrix to decompose - * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. - * By default, none is computed. This is a bit - field, the possible bits are #ComputeFullU, #ComputeThinU, - * #ComputeFullV, #ComputeThinV. - * - * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not - * available with the (non - default) FullPivHouseholderQR preconditioner. - */ - BDCSVD(const MatrixType& matrix, unsigned int computationOptions = 0) - : SVDBase<_MatrixType>::SVDBase(), - algoswap(ALGOSWAP) - { - compute(matrix, computationOptions); - } - - ~BDCSVD() - { - } - /** \brief Method performing the decomposition of given matrix using custom options. - * - * \param matrix the matrix to decompose - * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. - * By default, none is computed. This is a bit - field, the possible bits are #ComputeFullU, #ComputeThinU, - * #ComputeFullV, #ComputeThinV. - * - * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not - * available with the (non - default) FullPivHouseholderQR preconditioner. - */ - SVDBase<MatrixType>& compute(const MatrixType& matrix, unsigned int computationOptions); - - /** \brief Method performing the decomposition of given matrix using current options. - * - * \param matrix the matrix to decompose - * - * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int). - */ - SVDBase<MatrixType>& compute(const MatrixType& matrix) - { - return compute(matrix, this->m_computationOptions); - } - - void setSwitchSize(int s) - { - eigen_assert(s>3 && "BDCSVD the size of the algo switch has to be greater than 4"); - algoswap = s; - } - - - /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. - * - * \param b the right - hand - side of the equation to solve. - * - * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V. - * - * \note SVD solving is implicitly least - squares. Thus, this method serves both purposes of exact solving and least - squares solving. - * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. - */ - template<typename Rhs> - inline const internal::solve_retval<BDCSVD, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(this->m_isInitialized && "BDCSVD is not initialized."); - eigen_assert(SVDBase<_MatrixType>::computeU() && SVDBase<_MatrixType>::computeV() && - "BDCSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); - return internal::solve_retval<BDCSVD, Rhs>(*this, b.derived()); - } - - - const MatrixUType& matrixU() const - { - eigen_assert(this->m_isInitialized && "SVD is not initialized."); - if (isTranspose){ - eigen_assert(this->computeV() && "This SVD decomposition didn't compute U. Did you ask for it?"); - return this->m_matrixV; - } - else - { - eigen_assert(this->computeU() && "This SVD decomposition didn't compute U. Did you ask for it?"); - return this->m_matrixU; - } - - } - - - const MatrixVType& matrixV() const - { - eigen_assert(this->m_isInitialized && "SVD is not initialized."); - if (isTranspose){ - eigen_assert(this->computeU() && "This SVD decomposition didn't compute V. Did you ask for it?"); - return this->m_matrixU; - } - else - { - eigen_assert(this->computeV() && "This SVD decomposition didn't compute V. Did you ask for it?"); - return this->m_matrixV; - } - } - -private: - void allocate(Index rows, Index cols, unsigned int computationOptions); - void divide (Index firstCol, Index lastCol, Index firstRowW, - Index firstColW, Index shift); - void deflation43(Index firstCol, Index shift, Index i, Index size); - void deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size); - void deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift); - void copyUV(MatrixXr naiveU, MatrixXr naiveV, MatrixX householderU, MatrixX houseHolderV); - -protected: - MatrixXr m_naiveU, m_naiveV; - MatrixXr m_computed; - Index nRec; - int algoswap; - bool isTranspose, compU, compV; - -}; //end class BDCSVD - - -// Methode to allocate ans initialize matrix and attributs -template<typename MatrixType> -void BDCSVD<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions) -{ - isTranspose = (cols > rows); - if (SVDBase<MatrixType>::allocate(rows, cols, computationOptions)) return; - m_computed = MatrixXr::Zero(this->m_diagSize + 1, this->m_diagSize ); - if (isTranspose){ - compU = this->computeU(); - compV = this->computeV(); - } - else - { - compV = this->computeU(); - compU = this->computeV(); - } - if (compU) m_naiveU = MatrixXr::Zero(this->m_diagSize + 1, this->m_diagSize + 1 ); - else m_naiveU = MatrixXr::Zero(2, this->m_diagSize + 1 ); - - if (compV) m_naiveV = MatrixXr::Zero(this->m_diagSize, this->m_diagSize); - - - //should be changed for a cleaner implementation - if (isTranspose){ - bool aux; - if (this->computeU()||this->computeV()){ - aux = this->m_computeFullU; - this->m_computeFullU = this->m_computeFullV; - this->m_computeFullV = aux; - aux = this->m_computeThinU; - this->m_computeThinU = this->m_computeThinV; - this->m_computeThinV = aux; - } - } -}// end allocate - -// Methode which compute the BDCSVD for the int -template<> -SVDBase<Matrix<int, Dynamic, Dynamic> >& -BDCSVD<Matrix<int, Dynamic, Dynamic> >::compute(const MatrixType& matrix, unsigned int computationOptions) { - allocate(matrix.rows(), matrix.cols(), computationOptions); - this->m_nonzeroSingularValues = 0; - m_computed = Matrix<int, Dynamic, Dynamic>::Zero(rows(), cols()); - for (int i=0; i<this->m_diagSize; i++) { - this->m_singularValues.coeffRef(i) = 0; - } - if (this->m_computeFullU) this->m_matrixU = Matrix<int, Dynamic, Dynamic>::Zero(rows(), rows()); - if (this->m_computeFullV) this->m_matrixV = Matrix<int, Dynamic, Dynamic>::Zero(cols(), cols()); - this->m_isInitialized = true; - return *this; -} - - -// Methode which compute the BDCSVD -template<typename MatrixType> -SVDBase<MatrixType>& -BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions) -{ - allocate(matrix.rows(), matrix.cols(), computationOptions); - using std::abs; - - //**** step 1 Bidiagonalization isTranspose = (matrix.cols()>matrix.rows()) ; - MatrixType copy; - if (isTranspose) copy = matrix.adjoint(); - else copy = matrix; - - internal::UpperBidiagonalization<MatrixX > bid(copy); - - //**** step 2 Divide - // this is ugly and has to be redone (care of complex cast) - MatrixXr temp; - temp = bid.bidiagonal().toDenseMatrix().transpose(); - m_computed.setZero(); - for (int i=0; i<this->m_diagSize - 1; i++) { - m_computed(i, i) = temp(i, i); - m_computed(i + 1, i) = temp(i + 1, i); - } - m_computed(this->m_diagSize - 1, this->m_diagSize - 1) = temp(this->m_diagSize - 1, this->m_diagSize - 1); - divide(0, this->m_diagSize - 1, 0, 0, 0); - - //**** step 3 copy - for (int i=0; i<this->m_diagSize; i++) { - RealScalar a = abs(m_computed.coeff(i, i)); - this->m_singularValues.coeffRef(i) = a; - if (a == 0){ - this->m_nonzeroSingularValues = i; - break; - } - else if (i == this->m_diagSize - 1) - { - this->m_nonzeroSingularValues = i + 1; - break; - } - } - copyUV(m_naiveV, m_naiveU, bid.householderU(), bid.householderV()); - this->m_isInitialized = true; - return *this; -}// end compute - - -template<typename MatrixType> -void BDCSVD<MatrixType>::copyUV(MatrixXr naiveU, MatrixXr naiveV, MatrixX householderU, MatrixX householderV){ - if (this->computeU()){ - MatrixX temp = MatrixX::Zero(naiveU.rows(), naiveU.cols()); - temp.real() = naiveU; - if (this->m_computeThinU){ - this->m_matrixU = MatrixX::Identity(householderU.cols(), this->m_nonzeroSingularValues ); - this->m_matrixU.block(0, 0, this->m_diagSize, this->m_nonzeroSingularValues) = - temp.block(0, 0, this->m_diagSize, this->m_nonzeroSingularValues); - this->m_matrixU = householderU * this->m_matrixU ; - } - else - { - this->m_matrixU = MatrixX::Identity(householderU.cols(), householderU.cols()); - this->m_matrixU.block(0, 0, this->m_diagSize, this->m_diagSize) = temp.block(0, 0, this->m_diagSize, this->m_diagSize); - this->m_matrixU = householderU * this->m_matrixU ; - } - } - if (this->computeV()){ - MatrixX temp = MatrixX::Zero(naiveV.rows(), naiveV.cols()); - temp.real() = naiveV; - if (this->m_computeThinV){ - this->m_matrixV = MatrixX::Identity(householderV.cols(),this->m_nonzeroSingularValues ); - this->m_matrixV.block(0, 0, this->m_nonzeroSingularValues, this->m_nonzeroSingularValues) = - temp.block(0, 0, this->m_nonzeroSingularValues, this->m_nonzeroSingularValues); - this->m_matrixV = householderV * this->m_matrixV ; - } - else - { - this->m_matrixV = MatrixX::Identity(householderV.cols(), householderV.cols()); - this->m_matrixV.block(0, 0, this->m_diagSize, this->m_diagSize) = temp.block(0, 0, this->m_diagSize, this->m_diagSize); - this->m_matrixV = householderV * this->m_matrixV; - } - } -} - -// The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods takes as argument the -// place of the submatrix we are currently working on. - -//@param firstCol : The Index of the first column of the submatrix of m_computed and for m_naiveU; -//@param lastCol : The Index of the last column of the submatrix of m_computed and for m_naiveU; -// lastCol + 1 - firstCol is the size of the submatrix. -//@param firstRowW : The Index of the first row of the matrix W that we are to change. (see the reference paper section 1 for more information on W) -//@param firstRowW : Same as firstRowW with the column. -//@param shift : Each time one takes the left submatrix, one must add 1 to the shift. Why? Because! We actually want the last column of the U submatrix -// to become the first column (*coeff) and to shift all the other columns to the right. There are more details on the reference paper. -template<typename MatrixType> -void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW, - Index firstColW, Index shift) -{ - // requires nbRows = nbCols + 1; - using std::pow; - using std::sqrt; - using std::abs; - const Index n = lastCol - firstCol + 1; - const Index k = n/2; - RealScalar alphaK; - RealScalar betaK; - RealScalar r0; - RealScalar lambda, phi, c0, s0; - MatrixXr l, f; - // We use the other algorithm which is more efficient for small - // matrices. - if (n < algoswap){ - JacobiSVD<MatrixXr> b(m_computed.block(firstCol, firstCol, n + 1, n), - ComputeFullU | (ComputeFullV * compV)) ; - if (compU) m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() << b.matrixU(); - else - { - m_naiveU.row(0).segment(firstCol, n + 1).real() << b.matrixU().row(0); - m_naiveU.row(1).segment(firstCol, n + 1).real() << b.matrixU().row(n); - } - if (compV) m_naiveV.block(firstRowW, firstColW, n, n).real() << b.matrixV(); - m_computed.block(firstCol + shift, firstCol + shift, n + 1, n).setZero(); - for (int i=0; i<n; i++) - { - m_computed(firstCol + shift + i, firstCol + shift +i) = b.singularValues().coeffRef(i); - } - return; - } - // We use the divide and conquer algorithm - alphaK = m_computed(firstCol + k, firstCol + k); - betaK = m_computed(firstCol + k + 1, firstCol + k); - // The divide must be done in that order in order to have good results. Divide change the data inside the submatrices - // and the divide of the right submatrice reads one column of the left submatrice. That's why we need to treat the - // right submatrix before the left one. - divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift); - divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1); - if (compU) - { - lambda = m_naiveU(firstCol + k, firstCol + k); - phi = m_naiveU(firstCol + k + 1, lastCol + 1); - } - else - { - lambda = m_naiveU(1, firstCol + k); - phi = m_naiveU(0, lastCol + 1); - } - r0 = sqrt((abs(alphaK * lambda) * abs(alphaK * lambda)) - + abs(betaK * phi) * abs(betaK * phi)); - if (compU) - { - l = m_naiveU.row(firstCol + k).segment(firstCol, k); - f = m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1, n - k - 1); - } - else - { - l = m_naiveU.row(1).segment(firstCol, k); - f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1); - } - if (compV) m_naiveV(firstRowW+k, firstColW) = 1; - if (r0 == 0) - { - c0 = 1; - s0 = 0; - } - else - { - c0 = alphaK * lambda / r0; - s0 = betaK * phi / r0; - } - if (compU) - { - MatrixXr q1 (m_naiveU.col(firstCol + k).segment(firstCol, k + 1)); - // we shiftW Q1 to the right - for (Index i = firstCol + k - 1; i >= firstCol; i--) - { - m_naiveU.col(i + 1).segment(firstCol, k + 1) << m_naiveU.col(i).segment(firstCol, k + 1); - } - // we shift q1 at the left with a factor c0 - m_naiveU.col(firstCol).segment( firstCol, k + 1) << (q1 * c0); - // last column = q1 * - s0 - m_naiveU.col(lastCol + 1).segment(firstCol, k + 1) << (q1 * ( - s0)); - // first column = q2 * s0 - m_naiveU.col(firstCol).segment(firstCol + k + 1, n - k) << - m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) *s0; - // q2 *= c0 - m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) *= c0; - } - else - { - RealScalar q1 = (m_naiveU(0, firstCol + k)); - // we shift Q1 to the right - for (Index i = firstCol + k - 1; i >= firstCol; i--) - { - m_naiveU(0, i + 1) = m_naiveU(0, i); - } - // we shift q1 at the left with a factor c0 - m_naiveU(0, firstCol) = (q1 * c0); - // last column = q1 * - s0 - m_naiveU(0, lastCol + 1) = (q1 * ( - s0)); - // first column = q2 * s0 - m_naiveU(1, firstCol) = m_naiveU(1, lastCol + 1) *s0; - // q2 *= c0 - m_naiveU(1, lastCol + 1) *= c0; - m_naiveU.row(1).segment(firstCol + 1, k).setZero(); - m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1).setZero(); - } - m_computed(firstCol + shift, firstCol + shift) = r0; - m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) << alphaK * l.transpose().real(); - m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) << betaK * f.transpose().real(); - - - // the line below do the deflation of the matrix for the third part of the algorithm - // Here the deflation is commented because the third part of the algorithm is not implemented - // the third part of the algorithm is a fast SVD on the matrix m_computed which works thanks to the deflation - - deflation(firstCol, lastCol, k, firstRowW, firstColW, shift); - - // Third part of the algorithm, since the real third part of the algorithm is not implemeted we use a JacobiSVD - JacobiSVD<MatrixXr> res= JacobiSVD<MatrixXr>(m_computed.block(firstCol + shift, firstCol +shift, n + 1, n), - ComputeFullU | (ComputeFullV * compV)) ; - if (compU) m_naiveU.block(firstCol, firstCol, n + 1, n + 1) *= res.matrixU(); - else m_naiveU.block(0, firstCol, 2, n + 1) *= res.matrixU(); - - if (compV) m_naiveV.block(firstRowW, firstColW, n, n) *= res.matrixV(); - m_computed.block(firstCol + shift, firstCol + shift, n, n) << MatrixXr::Zero(n, n); - for (int i=0; i<n; i++) - m_computed(firstCol + shift + i, firstCol + shift +i) = res.singularValues().coeffRef(i); - // end of the third part - - -}// end divide - - -// page 12_13 -// i >= 1, di almost null and zi non null. -// We use a rotation to zero out zi applied to the left of M -template <typename MatrixType> -void BDCSVD<MatrixType>::deflation43(Index firstCol, Index shift, Index i, Index size){ - using std::abs; - using std::sqrt; - using std::pow; - RealScalar c = m_computed(firstCol + shift, firstCol + shift); - RealScalar s = m_computed(i, firstCol + shift); - RealScalar r = sqrt(pow(abs(c), 2) + pow(abs(s), 2)); - if (r == 0){ - m_computed(i, i)=0; - return; - } - c/=r; - s/=r; - m_computed(firstCol + shift, firstCol + shift) = r; - m_computed(i, firstCol + shift) = 0; - m_computed(i, i) = 0; - if (compU){ - m_naiveU.col(firstCol).segment(firstCol,size) = - c * m_naiveU.col(firstCol).segment(firstCol, size) - - s * m_naiveU.col(i).segment(firstCol, size) ; - - m_naiveU.col(i).segment(firstCol, size) = - (c + s*s/c) * m_naiveU.col(i).segment(firstCol, size) + - (s/c) * m_naiveU.col(firstCol).segment(firstCol,size); - } -}// end deflation 43 - - -// page 13 -// i,j >= 1, i != j and |di - dj| < epsilon * norm2(M) -// We apply two rotations to have zj = 0; -template <typename MatrixType> -void BDCSVD<MatrixType>::deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size){ - using std::abs; - using std::sqrt; - using std::conj; - using std::pow; - RealScalar c = m_computed(firstColm, firstColm + j - 1); - RealScalar s = m_computed(firstColm, firstColm + i - 1); - RealScalar r = sqrt(pow(abs(c), 2) + pow(abs(s), 2)); - if (r==0){ - m_computed(firstColm + i, firstColm + i) = m_computed(firstColm + j, firstColm + j); - return; - } - c/=r; - s/=r; - m_computed(firstColm + i, firstColm) = r; - m_computed(firstColm + i, firstColm + i) = m_computed(firstColm + j, firstColm + j); - m_computed(firstColm + j, firstColm) = 0; - if (compU){ - m_naiveU.col(firstColu + i).segment(firstColu, size) = - c * m_naiveU.col(firstColu + i).segment(firstColu, size) - - s * m_naiveU.col(firstColu + j).segment(firstColu, size) ; - - m_naiveU.col(firstColu + j).segment(firstColu, size) = - (c + s*s/c) * m_naiveU.col(firstColu + j).segment(firstColu, size) + - (s/c) * m_naiveU.col(firstColu + i).segment(firstColu, size); - } - if (compV){ - m_naiveV.col(firstColW + i).segment(firstRowW, size - 1) = - c * m_naiveV.col(firstColW + i).segment(firstRowW, size - 1) + - s * m_naiveV.col(firstColW + j).segment(firstRowW, size - 1) ; - - m_naiveV.col(firstColW + j).segment(firstRowW, size - 1) = - (c + s*s/c) * m_naiveV.col(firstColW + j).segment(firstRowW, size - 1) - - (s/c) * m_naiveV.col(firstColW + i).segment(firstRowW, size - 1); - } -}// end deflation 44 - - - -template <typename MatrixType> -void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift){ - //condition 4.1 - RealScalar EPS = EPSILON * (std::max<RealScalar>(m_computed(firstCol + shift + 1, firstCol + shift + 1), m_computed(firstCol + k, firstCol + k))); - const Index length = lastCol + 1 - firstCol; - if (m_computed(firstCol + shift, firstCol + shift) < EPS){ - m_computed(firstCol + shift, firstCol + shift) = EPS; - } - //condition 4.2 - for (Index i=firstCol + shift + 1;i<=lastCol + shift;i++){ - if (std::abs(m_computed(i, firstCol + shift)) < EPS){ - m_computed(i, firstCol + shift) = 0; - } - } - - //condition 4.3 - for (Index i=firstCol + shift + 1;i<=lastCol + shift; i++){ - if (m_computed(i, i) < EPS){ - deflation43(firstCol, shift, i, length); - } - } - - //condition 4.4 - - Index i=firstCol + shift + 1, j=firstCol + shift + k + 1; - //we stock the final place of each line - Index *permutation = new Index[length]; - - for (Index p =1; p < length; p++) { - if (i> firstCol + shift + k){ - permutation[p] = j; - j++; - } else if (j> lastCol + shift) - { - permutation[p] = i; - i++; - } - else - { - if (m_computed(i, i) < m_computed(j, j)){ - permutation[p] = j; - j++; - } - else - { - permutation[p] = i; - i++; - } - } - } - //we do the permutation - RealScalar aux; - //we stock the current index of each col - //and the column of each index - Index *realInd = new Index[length]; - Index *realCol = new Index[length]; - for (int pos = 0; pos< length; pos++){ - realCol[pos] = pos + firstCol + shift; - realInd[pos] = pos; - } - const Index Zero = firstCol + shift; - VectorType temp; - for (int i = 1; i < length - 1; i++){ - const Index I = i + Zero; - const Index realI = realInd[i]; - const Index j = permutation[length - i] - Zero; - const Index J = realCol[j]; - - //diag displace - aux = m_computed(I, I); - m_computed(I, I) = m_computed(J, J); - m_computed(J, J) = aux; - - //firstrow displace - aux = m_computed(I, Zero); - m_computed(I, Zero) = m_computed(J, Zero); - m_computed(J, Zero) = aux; - - // change columns - if (compU) { - temp = m_naiveU.col(I - shift).segment(firstCol, length + 1); - m_naiveU.col(I - shift).segment(firstCol, length + 1) << - m_naiveU.col(J - shift).segment(firstCol, length + 1); - m_naiveU.col(J - shift).segment(firstCol, length + 1) << temp; - } - else - { - temp = m_naiveU.col(I - shift).segment(0, 2); - m_naiveU.col(I - shift).segment(0, 2) << - m_naiveU.col(J - shift).segment(0, 2); - m_naiveU.col(J - shift).segment(0, 2) << temp; - } - if (compV) { - const Index CWI = I + firstColW - Zero; - const Index CWJ = J + firstColW - Zero; - temp = m_naiveV.col(CWI).segment(firstRowW, length); - m_naiveV.col(CWI).segment(firstRowW, length) << m_naiveV.col(CWJ).segment(firstRowW, length); - m_naiveV.col(CWJ).segment(firstRowW, length) << temp; - } - - //update real pos - realCol[realI] = J; - realCol[j] = I; - realInd[J - Zero] = realI; - realInd[I - Zero] = j; - } - for (Index i = firstCol + shift + 1; i<lastCol + shift;i++){ - if ((m_computed(i + 1, i + 1) - m_computed(i, i)) < EPS){ - deflation44(firstCol , - firstCol + shift, - firstRowW, - firstColW, - i - Zero, - i + 1 - Zero, - length); - } - } - delete [] permutation; - delete [] realInd; - delete [] realCol; - -}//end deflation - - -namespace internal{ - -template<typename _MatrixType, typename Rhs> -struct solve_retval<BDCSVD<_MatrixType>, Rhs> - : solve_retval_base<BDCSVD<_MatrixType>, Rhs> -{ - typedef BDCSVD<_MatrixType> BDCSVDType; - EIGEN_MAKE_SOLVE_HELPERS(BDCSVDType, Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - eigen_assert(rhs().rows() == dec().rows()); - // A = U S V^* - // So A^{ - 1} = V S^{ - 1} U^* - Index diagSize = (std::min)(dec().rows(), dec().cols()); - typename BDCSVDType::SingularValuesType invertedSingVals(diagSize); - Index nonzeroSingVals = dec().nonzeroSingularValues(); - invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse(); - invertedSingVals.tail(diagSize - nonzeroSingVals).setZero(); - - dst = dec().matrixV().leftCols(diagSize) - * invertedSingVals.asDiagonal() - * dec().matrixU().leftCols(diagSize).adjoint() - * rhs(); - return; - } -}; - -} //end namespace internal - - /** \svd_module - * - * \return the singular value decomposition of \c *this computed by - * BDC Algorithm - * - * \sa class BDCSVD - */ -/* -template<typename Derived> -BDCSVD<typename MatrixBase<Derived>::PlainObject> -MatrixBase<Derived>::bdcSvd(unsigned int computationOptions) const -{ - return BDCSVD<PlainObject>(*this, computationOptions); -} -*/ - -} // end namespace Eigen - -#endif |