src/SVD/BDCSVD.h

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