// 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