BDCSVD: restore perturbCol0 guard

libeigen/eigen!2440

Closes #1723

Eigen/src/SVD/BDCSVDImpl.h

diff --git a/Eigen/src/SVD/BDCSVDImpl.h b/Eigen/src/SVD/BDCSVDImpl.h
index 69e7156..7ba5bbf 100644
--- a/Eigen/src/SVD/BDCSVDImpl.h
+++ b/Eigen/src/SVD/BDCSVDImpl.h
@@ -580,7 +580,14 @@
       for (Index l = 0; l < m; ++l) {
         Index i = perm(l);
         if (i != k) {
-          Index j = i < k ? i : l > 0 ? perm(l - 1) : i;
+          // There is no valid predecessor when the first active index is already on the
+          // right of k. Treat this as a numerical issue and zero the product.
+          if (i >= k && l == 0) {
+            m_info = NumericalIssue;
+            prod = Literal(0);
+            break;
+          }
+          Index j = i < k ? i : perm(l - 1);
           prod *= ((singVals(j) + dk) / ((diag(i) + dk))) * ((mus(j) + (shifts(j) - dk)) / ((diag(i) - dk)));
         }
       }

test/bdcsvd.cpp

diff --git a/test/bdcsvd.cpp b/test/bdcsvd.cpp
index d6f639f..7ac28b0 100644
--- a/test/bdcsvd.cpp
+++ b/test/bdcsvd.cpp
@@ -16,7 +16,11 @@
 
 #include "main.h"
 #include "tridiag_test_matrices.h"
+#include <complex>
+#include <sstream>
+#define private public
 #include <Eigen/SVD>
+#undef private
 
 #define SVD_DEFAULT(M) BDCSVD<M>
 #define SVD_FOR_MIN_NORM(M) BDCSVD<M>
@@ -170,6 +174,28 @@
   }
 }
 
+void bdcsvd_perturbcol0_missing_predecessor() {
+  typedef Eigen::internal::bdcsvd_impl<double> Impl;
+
+  Impl impl;
+  Impl::ArrayXr col0(2), diag(2), shifts(2), mus(2), zhat(2);
+  Impl::ArrayXi perm(1);
+  Impl::VectorType singVals(2);
+
+  col0 << 1.0, 0.0;
+  diag << 1.0, 3.0;
+  singVals << 2.0, 4.0;
+  shifts << 0.0, 0.0;
+  mus << 1.0, 5.0;
+  perm << 1;
+
+  impl.perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat);
+
+  VERIFY(impl.info() == NumericalIssue);
+  VERIFY_IS_APPROX(zhat(0), 0.0);
+  VERIFY_IS_APPROX(zhat(1), 0.0);
+}
+
 EIGEN_DECLARE_TEST(bdcsvd) {
   CALL_SUBTEST_1((bdcsvd_verify_assert<Matrix3f>()));
   CALL_SUBTEST_2((bdcsvd_verify_assert<Matrix4d>()));
@@ -263,4 +289,5 @@
   // Bidiagonal SVD hard test cases
   CALL_SUBTEST_50((bdcsvd_bidiagonal_hard_cases<float>()));
   CALL_SUBTEST_51((bdcsvd_bidiagonal_hard_cases<double>()));
+  CALL_SUBTEST_52((bdcsvd_perturbcol0_missing_predecessor()));
 }