diff --git a/test/svd_common.h b/test/svd_common.h
index e902d23..347ea80 100644
--- a/test/svd_common.h
+++ b/test/svd_common.h
@@ -38,7 +38,6 @@
   sigma.diagonal() = svd.singularValues().template cast<Scalar>();
   MatrixUType u = svd.matrixU();
   MatrixVType v = svd.matrixV();
-
   VERIFY_IS_APPROX(m, u * sigma * v.adjoint());
   VERIFY_IS_UNITARY(u);
   VERIFY_IS_UNITARY(v);
@@ -90,31 +89,31 @@
   
   SolutionType x = svd.solve(rhs);
   
+  // evaluate normal equation which works also for least-squares solutions
+  if(internal::is_same<RealScalar,double>::value || svd.rank()==m.diagonal().size())
+  {
+    // This test is not stable with single precision.
+    // This is probably because squaring m signicantly affects the precision.
+    VERIFY_IS_APPROX(m.adjoint()*(m*x),m.adjoint()*rhs);
+  }
+  
   RealScalar residual = (m*x-rhs).norm();
   // Check that there is no significantly better solution in the neighborhood of x
   if(!test_isMuchSmallerThan(residual,rhs.norm()))
   {
-    // If the residual is very small, then we have an exact solution, so we are already good.
-    for(int k=0;k<x.rows();++k)
+    // ^^^ If the residual is very small, then we have an exact solution, so we are already good.
+    for(Index k=0;k<x.rows();++k)
     {
       SolutionType y(x);
-      y.row(k).array() += 2*NumTraits<RealScalar>::epsilon();
+      y.row(k) = (1.+2*NumTraits<RealScalar>::epsilon())*x.row(k);
       RealScalar residual_y = (m*y-rhs).norm();
       VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
       
-      y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon();
+      y.row(k) = (1.-2*NumTraits<RealScalar>::epsilon())*x.row(k);
       residual_y = (m*y-rhs).norm();
       VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
     }
   }
-  
-  // evaluate normal equation which works also for least-squares solutions
-  if(internal::is_same<RealScalar,double>::value)
-  {
-    // This test is not stable with single precision.
-    // This is probably because squaring m signicantly affects the precision.
-    VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
-  }
 }
 
 // check minimal norm solutions, the inoput matrix m is only used to recover problem size
@@ -234,11 +233,49 @@
   Matrix<RealScalar,Dynamic,1> d =  Matrix<RealScalar,Dynamic,1>::Random(diagSize);
   for(Index k=0; k<diagSize; ++k)
     d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s));
-  m = Matrix<Scalar,Dynamic,Dynamic>::Random(m.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,m.cols());
+
+  bool dup     = internal::random<int>(0,10) < 3;
+  bool unit_uv = internal::random<int>(0,10) < (dup?7:3); // if we duplicate some diagonal entries, then increase the chance to preserve them using unitary U and V factors
+  
+  // duplicate some singular values
+  if(dup)
+  {
+    Index n = internal::random<Index>(0,d.size()-1);
+    for(Index i=0; i<n; ++i)
+      d(internal::random<Index>(0,d.size()-1)) = d(internal::random<Index>(0,d.size()-1));
+  }
+  
+  Matrix<Scalar,Dynamic,Dynamic> U(m.rows(),diagSize);
+  Matrix<Scalar,Dynamic,Dynamic> VT(diagSize,m.cols());
+  if(unit_uv)
+  {
+    // in very rare cases let's try with a pure diagonal matrix
+    if(internal::random<int>(0,10) < 1)
+    {
+      U.setIdentity();
+      VT.setIdentity();
+    }
+    else
+    {
+      createRandomPIMatrixOfRank(diagSize,U.rows(), U.cols(), U);
+      createRandomPIMatrixOfRank(diagSize,VT.rows(), VT.cols(), VT);
+    }
+  }
+  else
+  {
+    U.setRandom();
+    VT.setRandom();
+  }
+  
+  m = U * d.asDiagonal() * VT;
+  
   // cancel some coeffs
-  Index n  = internal::random<Index>(0,m.size()-1);
-  for(Index i=0; i<n; ++i)
-    m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0);
+  if(!(dup && unit_uv))
+  {
+    Index n  = internal::random<Index>(0,m.size()-1);
+    for(Index i=0; i<n; ++i)
+      m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0);
+  }
 }
 
 
diff --git a/unsupported/Eigen/src/BDCSVD/BDCSVD.h b/unsupported/Eigen/src/BDCSVD/BDCSVD.h
index d5e8140..e8e7955 100644
--- a/unsupported/Eigen/src/BDCSVD/BDCSVD.h
+++ b/unsupported/Eigen/src/BDCSVD/BDCSVD.h
@@ -23,6 +23,10 @@
 // #define EIGEN_BDCSVD_SANITY_CHECKS
 namespace Eigen {
 
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+IOFormat bdcsvdfmt(8, 0, ", ", "\n", "  [", "]");
+#endif
+  
 template<typename _MatrixType> class BDCSVD;
 
 namespace internal {
@@ -80,6 +84,7 @@
   typedef Matrix<RealScalar, Dynamic, Dynamic> MatrixXr;
   typedef Matrix<RealScalar, Dynamic, 1> VectorType;
   typedef Array<RealScalar, Dynamic, 1> ArrayXr;
+  typedef Array<Index,1,Dynamic> ArrayXi;
 
   /** \brief Default Constructor.
    *
@@ -155,19 +160,16 @@
   void allocate(Index rows, Index cols, unsigned int computationOptions);
   void divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift);
   void computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V);
-  void computeSingVals(const ArrayXr& col0, const ArrayXr& diag, VectorType& singVals,
-                       ArrayXr& shifts, ArrayXr& mus);
-  void perturbCol0(const ArrayXr& col0, const ArrayXr& diag, const VectorType& singVals,
-                   const ArrayXr& shifts, const ArrayXr& mus, ArrayXr& zhat);
-  void computeSingVecs(const ArrayXr& zhat, const ArrayXr& diag, const VectorType& singVals,
-                       const ArrayXr& shifts, const ArrayXr& mus, MatrixXr& U, MatrixXr& V);
+  void computeSingVals(const ArrayXr& col0, const ArrayXr& diag, const ArrayXi& perm, VectorType& singVals, ArrayXr& shifts, ArrayXr& mus);
+  void perturbCol0(const ArrayXr& col0, const ArrayXr& diag, const ArrayXi& perm, const VectorType& singVals, const ArrayXr& shifts, const ArrayXr& mus, ArrayXr& zhat);
+  void computeSingVecs(const ArrayXr& zhat, const ArrayXr& diag, const ArrayXi& perm, const VectorType& singVals, const ArrayXr& shifts, const ArrayXr& mus, MatrixXr& U, MatrixXr& V);
   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);
   template<typename HouseholderU, typename HouseholderV, typename NaiveU, typename NaiveV>
   void copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naivev);
   static void structured_update(Block<MatrixXr,Dynamic,Dynamic> A, const MatrixXr &B, Index n1);
-  static RealScalar secularEq(RealScalar x, const ArrayXr& col0, const ArrayXr& diag, const ArrayXr& diagShifted, RealScalar shift, Index n);
+  static RealScalar secularEq(RealScalar x, const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm, const ArrayXr& diagShifted, RealScalar shift);
 
 protected:
   MatrixXr m_naiveU, m_naiveV;
@@ -234,6 +236,9 @@
 template<typename MatrixType>
 BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions) 
 {
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  std::cout << "\n\n\n======================================================================================================================\n\n\n";
+#endif
   allocate(matrix.rows(), matrix.cols(), computationOptions);
   using std::abs;
   
@@ -478,9 +483,23 @@
   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();
 
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  ArrayXr tmp1 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues();
+#endif
   // Second part: try to deflate singular values in combined matrix
   deflation(firstCol, lastCol, k, firstRowW, firstColW, shift);
-
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  ArrayXr tmp2 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues();
+  std::cout << "\n\nj1 = " << tmp1.transpose().format(bdcsvdfmt) << "\n";
+  std::cout << "j2 = " << tmp2.transpose().format(bdcsvdfmt) << "\n\n";
+  std::cout << "err:      " << ((tmp1-tmp2).abs()>1e-12*tmp2.abs()).transpose() << "\n";
+  static int count = 0;
+  std::cout << "# " << ++count << "\n\n";
+  assert((tmp1-tmp2).matrix().norm() < 1e-14*tmp2.matrix().norm());
+//   assert(count<681);
+//   assert(((tmp1-tmp2).abs()<1e-13*tmp2.abs()).all());
+#endif
+  
   // Third part: compute SVD of combined matrix
   MatrixXr UofSVD, VofSVD;
   VectorType singVals;
@@ -522,13 +541,24 @@
   ArrayXr diag = m_computed.block(firstCol, firstCol, n, n).diagonal();
   diag(0) = 0;
 
-  // compute singular values and vectors (in decreasing order)
+  // Allocate space for singular values and vectors
   singVals.resize(n);
   U.resize(n+1, n+1);
   if (m_compV) V.resize(n, n);
 
   if (col0.hasNaN() || diag.hasNaN()) { std::cout << "\n\nHAS NAN\n\n"; return; }
-
+  
+  // Many singular values might have been deflated, the zero ones have been moved to the end,
+  // but others are interleaved and we must ignore them at this stage.
+  // To this end, let's compute a permutation skipping them:
+  Index actual_n = n;
+  while(actual_n>1 && diag(actual_n-1)==0) --actual_n;
+  Index m = 0; // size of the deflated problem
+  ArrayXi perm(actual_n);
+  for(Index k=0;k<actual_n;++k)
+    if(col0(k)!=0)
+      perm(m++) = k;
+  perm.conservativeResize(m);
   
   ArrayXr shifts(n), mus(n), zhat(n);
 
@@ -539,12 +569,23 @@
 #endif
   
   // Compute singVals, shifts, and mus
-  computeSingVals(col0, diag, singVals, shifts, mus);
+  computeSingVals(col0, diag, perm, singVals, shifts, mus);
   
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  std::cout << "  j:        " << (m_computed.block(firstCol, firstCol, n, n)).jacobiSvd().singularValues().transpose().reverse() << "\n\n";
   std::cout << "  sing-val: " << singVals.transpose() << "\n";
   std::cout << "  mu:       " << mus.transpose() << "\n";
   std::cout << "  shift:    " << shifts.transpose() << "\n";
+  
+  {
+    Index actual_n = n;
+    while(actual_n>1 && col0(actual_n-1)==0) --actual_n;
+    std::cout << "\n\n    mus:    " << mus.head(actual_n).transpose() << "\n\n";
+    std::cout << "    check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n";
+    std::cout << "    check2 (>0)      : " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n";
+    std::cout << "    check3 (>0)      : " << ((diag.segment(1,actual_n-1)-singVals.head(actual_n-1).array()) / singVals.head(actual_n-1).array()).transpose() << "\n\n\n";
+    std::cout << "    check4 (>0)      : " << ((singVals.segment(1,actual_n-1)-singVals.head(actual_n-1))).transpose() << "\n\n\n";
+  }
 #endif
   
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
@@ -554,24 +595,16 @@
 #endif
   
   // Compute zhat
-  perturbCol0(col0, diag, singVals, shifts, mus, zhat);
+  perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat);
 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
   std::cout << "  zhat: " << zhat.transpose() << "\n";
-  {
-    Index actual_n = n;
-    while(actual_n>1 && col0(actual_n-1)==0) --actual_n;
-    std::cout << "\n\n    mus:    " << mus.head(actual_n).transpose() << "\n\n";
-    std::cout << "    check1: " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n";
-    std::cout << "    check2: " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n";
-    std::cout << "    check3: " << ((diag.segment(1,actual_n-1)-singVals.head(actual_n-1).array()) / singVals.head(actual_n-1).array()).transpose() << "\n\n\n";
-  }
 #endif
-
+  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   assert(zhat.allFinite());
 #endif
   
-  computeSingVecs(zhat, diag, singVals, shifts, mus, U, V);
+  computeSingVecs(zhat, diag, perm, singVals, shifts, mus, U, V);
   
 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
   std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n";
@@ -586,23 +619,48 @@
   assert(m_computed.allFinite());
 #endif
   
+  // Because of deflation, the singular values might not be completely sorted.
+  // Fortunately, reordering them is a O(n) problem
+  for(Index i=0; i<actual_n-1; ++i)
+  {
+    if(singVals(i)>singVals(i+1))
+    {
+      using std::swap;
+      swap(singVals(i),singVals(i+1));
+      U.col(i).swap(U.col(i+1));
+      if(m_compV) V.col(i).swap(V.col(i+1));
+    }
+  }
+  
   // Reverse order so that singular values in increased order
   // Because of deflation, the zeros singular-values are already at the end
-  Index actual_n = n;
-  while(actual_n>1 && singVals(actual_n-1)==0) --actual_n;
   singVals.head(actual_n).reverseInPlace();
   U.leftCols(actual_n) = U.leftCols(actual_n).rowwise().reverse().eval();               // FIXME this requires a temporary
   if (m_compV) V.leftCols(actual_n) = V.leftCols(actual_n).rowwise().reverse().eval();  // FIXME this requires a temporary
+  
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  JacobiSVD<MatrixXr> jsvd(m_computed.block(firstCol, firstCol, n, n) );
+  std::cout << "  * j:        " << jsvd.singularValues().transpose() << "\n\n";
+  std::cout << "  * sing-val: " << singVals.transpose() << "\n";
+//   std::cout << "  * err:      " << ((jsvd.singularValues()-singVals)>1e-13*singVals.norm()).transpose() << "\n";
+#endif
 }
 
 template <typename MatrixType>
-typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar mu, const ArrayXr& col0, const ArrayXr& diag, const ArrayXr& diagShifted, RealScalar shift, Index n)
+typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar mu, const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm, const ArrayXr& diagShifted, RealScalar shift)
 {
-  return 1 + (col0.square() / ((diagShifted - mu) )/( (diag + shift + mu))).head(n).sum();
+  Index m = perm.size();
+  RealScalar res = 1;
+  for(Index i=0; i<m; ++i)
+  {
+    Index j = perm(i);
+    res += numext::abs2(col0(j)) / ((diagShifted(j) - mu) * (diag(j) + shift + mu));
+  }
+  return res;
 }
 
 template <typename MatrixType>
-void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& diag, 
+void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm,
                                          VectorType& singVals, ArrayXr& shifts, ArrayXr& mus)
 {
   using std::abs;
@@ -612,26 +670,16 @@
   Index n = col0.size();
   Index actual_n = n;
   while(actual_n>1 && col0(actual_n-1)==0) --actual_n;
-//   Index m = 0;
-//   Array<Index,1,Dynamic> perm(actual_n);
-//   {
-//     for(Index k=0;k<actual_n;++k)
-//     {
-//       if(col0(k)!=0)
-//         perm(m++) = k;
-//     }
-//   }
-//   perm.conservativeResize(m);
-  
+
   for (Index k = 0; k < n; ++k)
   {
     if (col0(k) == 0 || actual_n==1)
     {
       // if col0(k) == 0, then entry is deflated, so singular value is on diagonal
       // if actual_n==1, then the deflated problem is already diagonalized
-      singVals(k) = diag(k);
+      singVals(k) = k==0 ? col0(0) : diag(k);
       mus(k) = 0;
-      shifts(k) = diag(k);
+      shifts(k) = k==0 ? col0(0) : diag(k);
       continue;
     } 
 
@@ -650,8 +698,22 @@
 
     // first decide whether it's closer to the left end or the right end
     RealScalar mid = left + (right-left) / 2;
-    RealScalar fMid = 1 + (col0.square() / ((diag + mid) * (diag - mid))).head(actual_n).sum();
-
+    RealScalar fMid = secularEq(mid, col0, diag, perm, diag, 0);
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+    std::cout << right-left << "\n";
+    std::cout << "fMid = " << fMid << " " << secularEq(mid-left, col0, diag, perm, diag-left, left) << " " << secularEq(mid-right, col0, diag, perm, diag-right, right)   << "\n";
+    std::cout << "     = " << secularEq(0.1*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.2*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.3*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.4*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.49*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.5*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.51*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.6*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.7*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.8*(left+right), col0, diag, perm, diag, 0)
+              << " "       << secularEq(0.9*(left+right), col0, diag, perm, diag, 0) << "\n";
+#endif
     RealScalar shift = (k == actual_n-1 || fMid > 0) ? left : right;
     
     // measure everything relative to shift
@@ -671,8 +733,8 @@
       muCur = -(right - left) * 0.5;
     }
 
-    RealScalar fPrev = secularEq(muPrev, col0, diag, diagShifted, shift, actual_n);
-    RealScalar fCur = secularEq(muCur, col0, diag, diagShifted, shift, actual_n);
+    RealScalar fPrev = secularEq(muPrev, col0, diag, perm, diagShifted, shift);
+    RealScalar fCur = secularEq(muCur, col0, diag, perm, diagShifted, shift);
     if (abs(fPrev) < abs(fCur))
     {
       swap(fPrev, fCur);
@@ -682,20 +744,26 @@
     // rational interpolation: fit a function of the form a / mu + b through the two previous
     // iterates and use its zero to compute the next iterate
     bool useBisection = fPrev*fCur>0;
-    while (abs(muCur - muPrev) > 8 * NumTraits<RealScalar>::epsilon() * (max)(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
+    while (fCur!=0 && abs(muCur - muPrev) > 8 * NumTraits<RealScalar>::epsilon() * (max)(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
     {
       ++m_numIters;
 
+      // Find a and b such that the function f(mu) = a / mu + b matches the current and previous samples.
       RealScalar a = (fCur - fPrev) / (1/muCur - 1/muPrev);
       RealScalar b = fCur - a / muCur;
-
+      // And find mu such that f(mu)==0:
+      RealScalar muZero = -a/b;
+      RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift);
+      
       muPrev = muCur;
       fPrev = fCur;
-      muCur = -a / b;
-      fCur = secularEq(muCur, col0, diag, diagShifted, shift, actual_n);
+      muCur = muZero;
+      fCur = fZero;
+      
       
       if (shift == left  && (muCur < 0 || muCur > right - left)) useBisection = true;
       if (shift == right && (muCur < -(right - left) || muCur > 0)) useBisection = true;
+      if (abs(fCur)>abs(fPrev)) useBisection = true;
     }
 
     // fall back on bisection method if rational interpolation did not work
@@ -710,7 +778,7 @@
         leftShifted = RealScalar(1)/NumTraits<RealScalar>::highest();
         // I don't understand why the case k==0 would be special there:
         // if (k == 0) rightShifted = right - left; else 
-        rightShifted = (right - left) * 0.6; // theoretically we can take 0.5, but let's be safe
+        rightShifted = (k==actual_n-1) ? right : ((right - left) * 0.6); // theoretically we can take 0.5, but let's be safe
       }
       else
       {
@@ -718,11 +786,11 @@
         rightShifted = -RealScalar(1)/NumTraits<RealScalar>::highest();
       }
       
-      RealScalar fLeft = secularEq(leftShifted, col0, diag, diagShifted, shift, actual_n);
-      RealScalar fRight = secularEq(rightShifted, col0, diag, diagShifted, shift, actual_n);
+      RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
+      RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift);
 
 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
-      if(fLeft * fRight>=0)
+      if(!(fLeft * fRight<0))
         std::cout << k << " : " <<  fLeft << " * " << fRight << " == " << fLeft * fRight << "  ;  " << left << " - " << right << " -> " <<  leftShifted << " " << rightShifted << "   shift=" << shift << "\n";
 #endif
       eigen_internal_assert(fLeft * fRight < 0);
@@ -730,7 +798,7 @@
       while (rightShifted - leftShifted > 2 * NumTraits<RealScalar>::epsilon() * (max)(abs(leftShifted), abs(rightShifted)))
       {
         RealScalar midShifted = (leftShifted + rightShifted) / 2;
-        RealScalar fMid = secularEq(midShifted, col0, diag, diagShifted, shift, actual_n);
+        RealScalar fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
         if (fLeft * fMid < 0)
         {
           rightShifted = midShifted;
@@ -762,28 +830,13 @@
 // zhat is perturbation of col0 for which singular vectors can be computed stably (see Section 3.1)
 template <typename MatrixType>
 void BDCSVD<MatrixType>::perturbCol0
-   (const ArrayXr& col0, const ArrayXr& diag, const VectorType& singVals,
+   (const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm, const VectorType& singVals,
     const ArrayXr& shifts, const ArrayXr& mus, ArrayXr& zhat)
 {
   using std::sqrt;
   Index n = col0.size();
-  
-  // Ignore trailing zeros:
-  Index actual_n = n;
-  while(actual_n>1 && col0(actual_n-1)==0) --actual_n;
-  // Deflated non-zero singular values are kept in-place,
-  // we thus compute an indirection array to properly ignore all deflated entries.
-  // TODO compute it once!
-  Index m = 0; // size of the deflated problem
-  Array<Index,1,Dynamic> perm(actual_n);
-  {
-    for(Index k=0;k<actual_n;++k)
-    {
-      if(col0(k)!=0)
-        perm(m++) = k;
-    }
-  }
-  perm.conservativeResize(m);
+  Index m = perm.size();
+  Index last = perm(m-1);
   
   // The offset permits to skip deflated entries while computing zhat
   for (Index k = 0; k < n; ++k)
@@ -794,7 +847,7 @@
     {
       // see equation (3.6)
       RealScalar dk = diag(k);
-      RealScalar prod = (singVals(actual_n-1) + dk) * (mus(actual_n-1) + (shifts(actual_n-1) - dk));
+      RealScalar prod = (singVals(last) + dk) * (mus(last) + (shifts(last) - dk));
 
       for(Index l = 0; l<m; ++l)
       {
@@ -805,12 +858,13 @@
           prod *= ((singVals(j)+dk) / ((diag(i)+dk))) * ((mus(j)+(shifts(j)-dk)) / ((diag(i)-dk)));
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
           if(i!=k && std::abs(((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) - 1) > 0.9 )
-            std::cout << "     " << ((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) << " == (" << (singVals(j)+dk) << " * " << (mus(j)+(shifts(j)-dk)) << ") / (" << (diag(i)+dk) << " * " << (diag(i)-dk) << ")\n";
+            std::cout << "     " << ((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) << " == (" << (singVals(j)+dk) << " * " << (mus(j)+(shifts(j)-dk))
+                       << ") / (" << (diag(i)+dk) << " * " << (diag(i)-dk) << ")\n";
 #endif
         }
       }
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
-      std::cout << "zhat(" << k << ") =  sqrt( " << prod << ")  ;  " << (singVals(actual_n-1) + dk) << " * " << mus(actual_n-1) + shifts(actual_n-1) << " - " << dk << "\n";
+      std::cout << "zhat(" << k << ") =  sqrt( " << prod << ")  ;  " << (singVals(last) + dk) << " * " << mus(last) + shifts(last) << " - " << dk << "\n";
 #endif
       RealScalar tmp = sqrt(prod);
       zhat(k) = col0(k) > 0 ? tmp : -tmp;
@@ -821,26 +875,11 @@
 // compute singular vectors
 template <typename MatrixType>
 void BDCSVD<MatrixType>::computeSingVecs
-   (const ArrayXr& zhat, const ArrayXr& diag, const VectorType& singVals,
+   (const ArrayXr& zhat, const ArrayXr& diag, const ArrayXi &perm, const VectorType& singVals,
     const ArrayXr& shifts, const ArrayXr& mus, MatrixXr& U, MatrixXr& V)
 {
   Index n = zhat.size();
-  
-  // Deflated non-zero singular values are kept in-place,
-  // we thus compute an indirection array to properly ignore all deflated entries.
-  // TODO compute it once!
-  Index actual_n = n;
-  while(actual_n>1 && zhat(actual_n-1)==0) --actual_n;
-  Index m = 0;
-  Array<Index,1,Dynamic> perm(actual_n);
-  {
-    for(Index k=0;k<actual_n;++k)
-    {
-      if(zhat(k)!=0)
-        perm(m++) = k;
-    }
-  }
-  perm.conservativeResize(m);
+  Index m = perm.size();
   
   for (Index k = 0; k < n; ++k)
   {
@@ -863,8 +902,6 @@
       if (m_compV)
       {
         V.col(k).setZero();
-//         for(Index i=1;i<actual_n;++i)
-//           V(i,k) = diag(i) * zhat(i) / (((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
         for(Index l=1;l<m;++l)
         {
           Index i = perm(l);
@@ -908,7 +945,7 @@
 
 
 // page 13
-// i,j >= 1, i != j and |di - dj| < epsilon * norm2(M)
+// i,j >= 1, i!=j and |di - dj| < epsilon * norm2(M)
 // We apply two rotations to have zj = 0;
 // TODO deflation44 is still broken and not properly tested
 template <typename MatrixType>
@@ -940,18 +977,13 @@
   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 + j) = m_computed(firstColm + i, firstColm + i);
   m_computed(firstColm + j, firstColm) = 0;
 
-  JacobiRotation<RealScalar> J(c,s);
-  if (m_compU)
-  {
-    m_naiveU.middleRows(firstColu, size).applyOnTheRight(firstColu + i, firstColu + j, J);
-  } 
-  if (m_compV)
-  {
-    m_naiveU.middleRows(firstRowW, size-1).applyOnTheRight(firstColW + i, firstColW + j, J.transpose());
-  }
+  JacobiRotation<RealScalar> J(c,-s);
+  if (m_compU)  m_naiveU.middleRows(firstColu, size+1).applyOnTheRight(firstColu + i, firstColu + j, J);
+  else          m_naiveU.applyOnTheRight(firstColu+i, firstColu+j, J);
+  if (m_compV)  m_naiveV.middleRows(firstRowW, size).applyOnTheRight(firstColW + i, firstColW + j, J);
 }// end deflation 44
 
 
@@ -964,43 +996,49 @@
   using std::max;
   const Index length = lastCol + 1 - firstCol;
   
+  Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1);
+  Diagonal<MatrixXr> fulldiag(m_computed);
+  VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length);
+  
+  RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff();
+  RealScalar epsilon_strict = NumTraits<RealScalar>::epsilon() * maxDiag;
+  RealScalar epsilon_coarse = 8 * NumTraits<RealScalar>::epsilon() * (max)(col0.cwiseAbs().maxCoeff(), maxDiag);
+  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   assert(m_naiveU.allFinite());
   assert(m_naiveV.allFinite());
   assert(m_computed.allFinite());
 #endif
-  
-  Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1);
-  Diagonal<MatrixXr> fulldiag(m_computed);
-  VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length);
-  
-  RealScalar epsilon = 8 * NumTraits<RealScalar>::epsilon() * (max)(col0.cwiseAbs().maxCoeff(), diag.cwiseAbs().maxCoeff());
+
+#ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE  
+  std::cout << "\ndeflate:" << diag.head(k+1).transpose() << "  |  " << diag.segment(k+1,length-k-1).transpose() << "\n";
+#endif
   
   //condition 4.1
-  if (diag(0) < epsilon)
+  if (diag(0) < epsilon_coarse)
   { 
 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
-    std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon << "\n";
+    std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon_coarse << "\n";
 #endif
-    diag(0) = epsilon;
+    diag(0) = epsilon_coarse;
   }
 
   //condition 4.2
   for (Index i=1;i<length;++i)
-    if (abs(col0(i)) < epsilon)
+    if (abs(col0(i)) < epsilon_strict)
     {
 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
-      std::cout << "deflation 4.2, set z(" << i << ") to zero because " << abs(col0(i)) << " < " << epsilon << "  (diag(" << i << ")=" << diag(i) << ")\n";
+      std::cout << "deflation 4.2, set z(" << i << ") to zero because " << abs(col0(i)) << " < " << epsilon_strict << "  (diag(" << i << ")=" << diag(i) << ")\n";
 #endif
       col0(i) = 0;
     }
 
   //condition 4.3
   for (Index i=1;i<length; i++)
-    if (diag(i) < epsilon)
+    if (diag(i) < epsilon_coarse)
     {
 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
-      std::cout << "deflation 4.3, cancel z(" << i << ") because " << diag(i) << " < " << epsilon << "  (z[" << i << "]=" << col0(i) << ")\n";
+      std::cout << "deflation 4.3, cancel z(" << i << ")=" << col0(i) << " because diag(" << i << ")=" << diag(i) << " < " << epsilon_coarse << "\n";
 #endif
       deflation43(firstCol, shift, i, length);
     }
@@ -1010,14 +1048,22 @@
   assert(m_naiveV.allFinite());
   assert(m_computed.allFinite());
 #endif
-  
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  std::cout << "to be sorted: " << diag.transpose() << "\n\n";
+#endif
   {
+    // Check for total deflation
+    // If we have a total deflation, then we have to consider col0(0)==diag(0) as a singular value during sorting
+    bool total_deflation = (col0.tail(length-1).array()==RealScalar(0)).all();
+    
     // Sort the diagonal entries, since diag(1:k-1) and diag(k:length) are already sorted, let's do a sorted merge.
     // First, compute the respective permutation.
     Index *permutation = new Index[length]; // FIXME avoid repeated dynamic memory allocation
     {
       permutation[0] = 0;
       Index p = 1;
+      
+      // Move deflated diagonal entries at the end.
       for(Index i=1; i<length; ++i)
         if(diag(i)==0)
           permutation[p++] = i;
@@ -1032,6 +1078,22 @@
       }
     }
     
+    // If we have a total deflation, then we have to insert diag(0) at the right place
+    if(total_deflation)
+    {
+      for(Index i=1; i<length; ++i)
+      {
+        Index pi = permutation[i];
+        if(diag(pi)==0 || diag(0)<diag(pi))
+          permutation[i-1] = permutation[i];
+        else
+        {
+          permutation[i-1] = 0;
+          break;
+        }
+      }
+    }
+    
     // Current index of each col, and current column of each index
     Index *realInd = new Index[length];  // FIXME avoid repeated dynamic memory allocation
     Index *realCol = new Index[length];  // FIXME avoid repeated dynamic memory allocation
@@ -1042,15 +1104,15 @@
       realInd[pos] = pos;
     }
     
-    for(Index i = 1; i < length; i++)
+    for(Index i = total_deflation?0:1; i < length; i++)
     {
-      const Index pi = permutation[length - i];
+      const Index pi = permutation[length - (total_deflation ? i+1 : i)];
       const Index J = realCol[pi];
       
       using std::swap;
-      // swap diaognal and first column entries:
+      // swap diagonal and first column entries:
       swap(diag(i), diag(J));
-      swap(col0(i), col0(J));
+      if(i!=0 && J!=0) swap(col0(i), col0(J));
 
       // change columns
       if (m_compU) m_naiveU.col(firstCol+i).segment(firstCol, length + 1).swap(m_naiveU.col(firstCol+J).segment(firstCol, length + 1));
@@ -1068,23 +1130,31 @@
     delete[] realInd;
     delete[] realCol;
   }
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  std::cout << "sorted: " << diag.transpose().format(bdcsvdfmt) << "\n";
+  std::cout << "      : " << col0.transpose() << "\n\n";
+#endif
+    
+  //condition 4.4
+  {
+    Index i = length-1;
+    while(i>0 && (diag(i)==0 || col0(i)==0)) --i;
+    for(; i>1;--i)
+       if( (diag(i) - diag(i-1)) < NumTraits<RealScalar>::epsilon()*maxDiag )
+      {
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+        std::cout << "deflation 4.4 with i = " << i << " because " << (diag(i) - diag(i-1)) << " < " << NumTraits<RealScalar>::epsilon()*diag(i) << "\n";
+#endif
+        eigen_internal_assert(abs(diag(i) - diag(i-1))<epsilon_coarse && " diagonal entries are not properly sorted");
+        deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i-1, i, length);
+      }
+  }
   
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
-  for(int k=2;k<length;++k)
-    assert(diag(k-1)<=diag(k) || diag(k)==0);
+  for(Index j=2;j<length;++j)
+    assert(diag(j-1)<=diag(j) || diag(j)==0);
 #endif
   
-  //condition 4.4
-  for (Index i = 1; i+1<length && diag(i+1)!=0 && col0(i+1)!=0;i++)
-    if ((diag(i+1) - diag(i)) < NumTraits<RealScalar>::epsilon()*diag(i+1))
-//     if ((diag(i+1) - diag(i)) < epsilon)
-    {
-#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
-      std::cout << "deflation 4.4 with i = " << i << " because " << (diag(i+1) - diag(i)) << " < " << epsilon << "\n";
-#endif
-      deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i, i + 1, length);
-    }
-  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   assert(m_naiveU.allFinite());
   assert(m_naiveV.allFinite());