20 #ifndef EIGEN_BDCSVD_H
21 #define EIGEN_BDCSVD_H
25 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
26 #undef eigen_internal_assert
27 #define eigen_internal_assert(X) assert(X);
32 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
33 IOFormat bdcsvdfmt(8, 0,
", ",
"\n",
" [",
"]");
36 template<
typename _MatrixType>
class BDCSVD;
40 template<
typename _MatrixType>
72 template<
typename _MatrixType>
175 eigen_assert(
s>3 &&
"BDCSVD the size of the algo switch has to be greater than 3");
189 template<
typename HouseholderU,
typename HouseholderV,
typename NaiveU,
typename NaiveV>
190 void copyUV(
const HouseholderU &householderU,
const HouseholderV &householderV,
const NaiveU &naiveU,
const NaiveV &naivev);
221 template<
typename MatrixType>
226 if (Base::allocate(
rows,
cols, computationOptions))
229 m_computed = MatrixXr::Zero(m_diagSize + 1, m_diagSize );
230 m_compU = computeV();
231 m_compV = computeU();
235 if (m_compU) m_naiveU = MatrixXr::Zero(m_diagSize + 1, m_diagSize + 1 );
236 else m_naiveU = MatrixXr::Zero(2, m_diagSize + 1 );
238 if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize);
240 m_workspace.resize((m_diagSize+1)*(m_diagSize+1)*3);
241 m_workspaceI.resize(3*m_diagSize);
244 template<
typename MatrixType>
247 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
248 std::cout <<
"\n\n\n======================================================================================================================\n\n\n";
250 allocate(
matrix.rows(),
matrix.cols(), computationOptions);
256 if(
matrix.cols() < m_algoswap)
260 m_isInitialized =
true;
261 m_info = jsvd.info();
263 if(computeU()) m_matrixU = jsvd.matrixU();
264 if(computeV()) m_matrixV = jsvd.matrixV();
265 m_singularValues = jsvd.singularValues();
266 m_nonzeroSingularValues = jsvd.nonzeroSingularValues();
274 m_isInitialized =
true;
292 m_computed.topRows(m_diagSize) = bid.
bidiagonal().toDenseMatrix().transpose();
293 m_computed.template bottomRows<1>().setZero();
294 divide(0, m_diagSize - 1, 0, 0, 0);
296 m_isInitialized =
true;
301 for (
int i=0;
i<m_diagSize;
i++)
304 m_singularValues.coeffRef(
i) =
a *
scale;
307 m_nonzeroSingularValues =
i;
308 m_singularValues.tail(m_diagSize -
i - 1).setZero();
311 else if (
i == m_diagSize - 1)
313 m_nonzeroSingularValues =
i + 1;
318 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
325 m_isInitialized =
true;
330 template<
typename MatrixType>
331 template<
typename HouseholderU,
typename HouseholderV,
typename NaiveU,
typename NaiveV>
332 void BDCSVD<MatrixType>::copyUV(
const HouseholderU &householderU,
const HouseholderV &householderV,
const NaiveU &naiveU,
const NaiveV &naiveV)
337 Index Ucols = m_computeThinU ? m_diagSize : householderU.cols();
338 m_matrixU = MatrixX::Identity(householderU.cols(), Ucols);
339 m_matrixU.topLeftCorner(m_diagSize, m_diagSize) = naiveV.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
340 householderU.applyThisOnTheLeft(m_matrixU);
344 Index Vcols = m_computeThinV ? m_diagSize : householderV.cols();
345 m_matrixV = MatrixX::Identity(householderV.cols(), Vcols);
346 m_matrixV.topLeftCorner(m_diagSize, m_diagSize) = naiveU.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
347 householderV.applyThisOnTheLeft(m_matrixV);
359 template<
typename MatrixType>
375 if( (
A.col(
j).head(
n1).array()!=
Literal(0)).any() )
378 B1.row(
k1) =
B.row(
j);
381 if( (
A.col(
j).tail(
n2).array()!=
Literal(0)).any() )
383 A2.col(k2) =
A.col(
j).tail(
n2);
384 B2.row(k2) =
B.row(
j);
389 A.topRows(
n1).noalias() =
A1.leftCols(
k1) * B1.topRows(
k1);
390 A.bottomRows(
n2).noalias() =
A2.leftCols(k2) * B2.topRows(k2);
410 template<
typename MatrixType>
417 const Index n = lastCol - firstCol + 1;
434 m_naiveU.block(firstCol, firstCol,
n + 1,
n + 1).real() =
b.matrixU();
437 m_naiveU.row(0).segment(firstCol,
n + 1).real() =
b.matrixU().row(0);
438 m_naiveU.row(1).segment(firstCol,
n + 1).real() =
b.matrixU().row(
n);
440 if (m_compV) m_naiveV.block(firstRowW, firstColW,
n,
n).real() =
b.matrixV();
441 m_computed.block(firstCol + shift, firstCol + shift,
n + 1,
n).setZero();
442 m_computed.diagonal().segment(firstCol + shift,
n) =
b.singularValues().head(
n);
446 alphaK = m_computed(firstCol + k, firstCol + k);
447 betaK = m_computed(firstCol + k + 1, firstCol + k);
451 divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift);
453 divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1);
458 lambda = m_naiveU(firstCol + k, firstCol + k);
459 phi = m_naiveU(firstCol + k + 1, lastCol + 1);
463 lambda = m_naiveU(1, firstCol + k);
464 phi = m_naiveU(0, lastCol + 1);
469 l = m_naiveU.row(firstCol + k).segment(firstCol, k);
470 f = m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1,
n - k - 1);
474 l = m_naiveU.row(1).segment(firstCol, k);
475 f = m_naiveU.row(0).segment(firstCol + k + 1,
n - k - 1);
477 if (m_compV) m_naiveV(firstRowW+k, firstColW) =
Literal(1);
485 c0 = alphaK *
lambda / r0;
486 s0 = betaK * phi / r0;
489 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
490 assert(m_naiveU.allFinite());
491 assert(m_naiveV.allFinite());
492 assert(m_computed.allFinite());
497 MatrixXr q1 (m_naiveU.col(firstCol + k).segment(firstCol, k + 1));
499 for (
Index i = firstCol + k - 1;
i >= firstCol;
i--)
500 m_naiveU.col(
i + 1).segment(firstCol, k + 1) = m_naiveU.col(
i).segment(firstCol, k + 1);
502 m_naiveU.col(firstCol).segment( firstCol, k + 1) = (q1 * c0);
504 m_naiveU.col(lastCol + 1).segment(firstCol, k + 1) = (q1 * ( - s0));
506 m_naiveU.col(firstCol).segment(firstCol + k + 1,
n - k) = m_naiveU.col(lastCol + 1).segment(firstCol + k + 1,
n - k) * s0;
508 m_naiveU.col(lastCol + 1).segment(firstCol + k + 1,
n - k) *= c0;
514 for (
Index i = firstCol + k - 1;
i >= firstCol;
i--)
515 m_naiveU(0,
i + 1) = m_naiveU(0,
i);
517 m_naiveU(0, firstCol) = (q1 * c0);
519 m_naiveU(0, lastCol + 1) = (q1 * ( - s0));
521 m_naiveU(1, firstCol) = m_naiveU(1, lastCol + 1) *s0;
523 m_naiveU(1, lastCol + 1) *= c0;
524 m_naiveU.row(1).segment(firstCol + 1, k).setZero();
525 m_naiveU.row(0).segment(firstCol + k + 1,
n - k - 1).setZero();
528 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
529 assert(m_naiveU.allFinite());
530 assert(m_naiveV.allFinite());
531 assert(m_computed.allFinite());
534 m_computed(firstCol + shift, firstCol + shift) = r0;
535 m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) = alphaK *
l.transpose().real();
536 m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1,
n - k - 1) = betaK *
f.transpose().real();
538 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
539 ArrayXr tmp1 = (m_computed.block(firstCol+shift, firstCol+shift,
n,
n)).jacobiSvd().singularValues();
542 deflation(firstCol, lastCol, k, firstRowW, firstColW, shift);
543 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
544 ArrayXr tmp2 = (m_computed.block(firstCol+shift, firstCol+shift,
n,
n)).jacobiSvd().singularValues();
545 std::cout <<
"\n\nj1 = " << tmp1.transpose().format(bdcsvdfmt) <<
"\n";
546 std::cout <<
"j2 = " << tmp2.transpose().format(bdcsvdfmt) <<
"\n\n";
547 std::cout <<
"err: " << ((tmp1-tmp2).
abs()>1
e-12*tmp2.abs()).transpose() <<
"\n";
548 static int count = 0;
549 std::cout <<
"# " << ++count <<
"\n\n";
550 assert((tmp1-tmp2).
matrix().norm() < 1
e-14*tmp2.matrix().norm());
558 computeSVDofM(firstCol + shift,
n, UofSVD, singVals, VofSVD);
560 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
561 assert(UofSVD.allFinite());
562 assert(VofSVD.allFinite());
566 structured_update(m_naiveU.block(firstCol, firstCol,
n + 1,
n + 1), UofSVD, (
n+2)/2);
570 tmp.noalias() = m_naiveU.middleCols(firstCol,
n+1) * UofSVD;
571 m_naiveU.middleCols(firstCol,
n + 1) = tmp;
574 if (m_compV) structured_update(m_naiveV.block(firstRowW, firstColW,
n,
n), VofSVD, (
n+1)/2);
576 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
577 assert(m_naiveU.allFinite());
578 assert(m_naiveV.allFinite());
579 assert(m_computed.allFinite());
582 m_computed.block(firstCol + shift, firstCol + shift,
n,
n).setZero();
583 m_computed.block(firstCol + shift, firstCol + shift,
n,
n).diagonal() = singVals;
594 template <
typename MatrixType>
599 ArrayRef col0 = m_computed.col(firstCol).segment(firstCol,
n);
600 m_workspace.head(
n) = m_computed.block(firstCol, firstCol,
n,
n).diagonal();
607 if (m_compV)
V.resize(
n,
n);
609 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
610 if (col0.hasNaN() ||
diag.hasNaN())
611 std::cout <<
"\n\nHAS NAN\n\n";
620 for(
Index k=0;k<actual_n;++k)
621 if(
abs(col0(k))>considerZero)
622 m_workspaceI(
m++) = k;
629 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
630 std::cout <<
"computeSVDofM using:\n";
631 std::cout <<
" z: " << col0.transpose() <<
"\n";
632 std::cout <<
" d: " <<
diag.transpose() <<
"\n";
636 computeSingVals(col0,
diag,
perm, singVals, shifts, mus);
638 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
639 std::cout <<
" j: " << (m_computed.block(firstCol, firstCol,
n,
n)).jacobiSvd().singularValues().transpose().reverse() <<
"\n\n";
640 std::cout <<
" sing-val: " << singVals.transpose() <<
"\n";
641 std::cout <<
" mu: " << mus.transpose() <<
"\n";
642 std::cout <<
" shift: " << shifts.transpose() <<
"\n";
645 std::cout <<
"\n\n mus: " << mus.head(actual_n).transpose() <<
"\n\n";
646 std::cout <<
" check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).
head(actual_n).transpose() <<
"\n\n";
647 assert((((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n) >= 0).all());
648 std::cout <<
" check2 (>0) : " << ((singVals.array()-
diag) / singVals.array()).
head(actual_n).transpose() <<
"\n\n";
649 assert((((singVals.array()-
diag) / singVals.array()).head(actual_n) >= 0).all());
653 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
654 assert(singVals.allFinite());
655 assert(mus.allFinite());
656 assert(shifts.allFinite());
660 perturbCol0(col0,
diag,
perm, singVals, shifts, mus, zhat);
661 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
662 std::cout <<
" zhat: " << zhat.transpose() <<
"\n";
665 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
666 assert(zhat.allFinite());
669 computeSingVecs(zhat,
diag,
perm, singVals, shifts, mus,
U,
V);
671 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
672 std::cout <<
"U^T U: " << (
U.transpose() *
U -
MatrixXr(MatrixXr::Identity(
U.cols(),
U.cols()))).norm() <<
"\n";
673 std::cout <<
"V^T V: " << (
V.transpose() *
V -
MatrixXr(MatrixXr::Identity(
V.cols(),
V.cols()))).norm() <<
"\n";
676 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
677 assert(m_naiveU.allFinite());
678 assert(m_naiveV.allFinite());
679 assert(m_computed.allFinite());
680 assert(
U.allFinite());
681 assert(
V.allFinite());
690 if(singVals(
i)>singVals(
i+1))
693 swap(singVals(
i),singVals(
i+1));
694 U.col(
i).swap(
U.col(
i+1));
695 if(m_compV)
V.col(
i).swap(
V.col(
i+1));
699 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
701 bool singular_values_sorted = (((singVals.segment(1,actual_n-1)-singVals.head(actual_n-1))).array() >= 0).
all();
702 if(!singular_values_sorted)
703 std::cout <<
"Singular values are not sorted: " << singVals.segment(1,actual_n).transpose() <<
"\n";
704 assert(singular_values_sorted);
710 singVals.head(actual_n).reverseInPlace();
711 U.leftCols(actual_n).rowwise().reverseInPlace();
712 if (m_compV)
V.leftCols(actual_n).rowwise().reverseInPlace();
714 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
716 std::cout <<
" * j: " << jsvd.singularValues().transpose() <<
"\n\n";
717 std::cout <<
" * sing-val: " << singVals.transpose() <<
"\n";
722 template <
typename MatrixType>
732 res += (col0(
j) / (diagShifted(
j) -
mu)) * (col0(
j) / (
diag(
j) + shift +
mu));
738 template <
typename MatrixType>
750 while(actual_n>1 && col0(actual_n-1)==
Literal(0)) --actual_n;
752 for (
Index k = 0; k <
n; ++k)
754 if (col0(k) ==
Literal(0) || actual_n==1)
758 singVals(k) = k==0 ? col0(0) :
diag(k);
760 shifts(k) = k==0 ? col0(0) :
diag(k);
768 right = (
diag(actual_n-1) + col0.matrix().norm());
782 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
783 std::cout <<
"right-left = " <<
right-
left <<
"\n";
804 diagShifted =
diag - shift;
811 midShifted = -midShifted;
817 diagShifted =
diag - shift;
826 if (k == actual_n-1) muCur =
right -
left;
837 if (
abs(fPrev) <
abs(fCur))
845 bool useBisection = fPrev*fCur>
Literal(0);
857 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
868 if (
abs(fCur)>
abs(fPrev)) useBisection =
true;
874 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
875 std::cout <<
"useBisection for k = " << k <<
", actual_n = " << actual_n <<
"\n";
902 #if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_SANITY_CHECKS
906 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
908 std::cout <<
"f(" << leftShifted <<
") =" << fLeft <<
" ; " <<
left <<
" " << shift <<
" " <<
right <<
"\n";
912 std::cout <<
"f(" << rightShifted <<
") =" << fRight <<
" ; " <<
left <<
" " << shift <<
" " <<
right <<
"\n";
916 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
917 if(!(fLeft * fRight<0))
919 std::cout <<
"f(leftShifted) using leftShifted=" << leftShifted <<
" ; diagShifted(1:10):" << diagShifted.head(10).transpose() <<
"\n ; "
920 <<
"left==shift=" << bool(
left==shift) <<
" ; left-shift = " << (
left-shift) <<
"\n";
921 std::cout <<
"k=" << k <<
", " << fLeft <<
" * " << fRight <<
" == " << fLeft * fRight <<
" ; "
922 <<
"[" <<
left <<
" .. " <<
right <<
"] -> [" << leftShifted <<
" " << rightShifted <<
"], shift=" << shift
923 <<
" , f(right)=" << secularEq(0, col0,
diag,
perm, diagShifted, shift)
924 <<
" == " << secularEq(
right, col0,
diag,
perm,
diag, 0) <<
" == " << fRight <<
"\n";
934 fMid = secularEq(midShifted, col0,
diag,
perm, diagShifted, shift);
939 rightShifted = midShifted;
943 leftShifted = midShifted;
947 muCur = (leftShifted + rightShifted) /
Literal(2);
961 singVals[k] = shift + muCur;
965 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
967 std::cout <<
"found " << singVals[k] <<
" == " << shift <<
" + " << muCur <<
" from " <<
diag(k) <<
" .. " <<
diag(k+1) <<
"\n";
969 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
970 assert(k==0 || singVals[k]>=singVals[k-1]);
971 assert(singVals[k]>=
diag(k));
984 template <
typename MatrixType>
999 for (
Index k = 0; k <
n; ++k)
1007 RealScalar prod = (singVals(lastIdx) + dk) * (mus(lastIdx) + (shifts(lastIdx) - dk));
1008 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1010 std::cout <<
"k = " << k <<
" ; z(k)=" << col0(k) <<
", diag(k)=" << dk <<
"\n";
1011 std::cout <<
"prod = " <<
"(" << singVals(lastIdx) <<
" + " << dk <<
") * (" << mus(lastIdx) <<
" + (" << shifts(lastIdx) <<
" - " << dk <<
"))" <<
"\n";
1012 std::cout <<
" = " << singVals(lastIdx) + dk <<
" * " << mus(lastIdx) + (shifts(lastIdx) - dk) <<
"\n";
1022 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1023 if(
i>=k && (
l==0 ||
l-1>=
m))
1025 std::cout <<
"Error in perturbCol0\n";
1026 std::cout <<
" " << k <<
"/" <<
n <<
" " <<
l <<
"/" <<
m <<
" " <<
i <<
"/" <<
n <<
" ; " << col0(k) <<
" " <<
diag(k) <<
" " <<
"\n";
1027 std::cout <<
" " <<
diag(
i) <<
"\n";
1029 std::cout <<
" " <<
"j=" <<
j <<
"\n";
1033 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1036 std::cout <<
"k=" << k <<
", i=" <<
i <<
", l=" <<
l <<
", perm.size()=" <<
perm.size() <<
"\n";
1040 prod *= ((singVals(
j)+dk) / ((
diag(
i)+dk))) * ((mus(
j)+(shifts(
j)-dk)) / ((
diag(
i)-dk)));
1041 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1044 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1046 std::cout <<
" " << ((singVals(
j)+dk)*(mus(
j)+(shifts(
j)-dk)))/((
diag(
i)+dk)*(
diag(
i)-dk)) <<
" == (" << (singVals(
j)+dk) <<
" * " << (mus(
j)+(shifts(
j)-dk))
1047 <<
") / (" << (
diag(
i)+dk) <<
" * " << (
diag(
i)-dk) <<
")\n";
1051 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1052 std::cout <<
"zhat(" << k <<
") = sqrt( " <<
prod <<
") ; " << (singVals(lastIdx) + dk) <<
" * " << mus(lastIdx) + shifts(lastIdx) <<
" - " << dk <<
"\n";
1055 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1064 template <
typename MatrixType>
1072 for (
Index k = 0; k <
n; ++k)
1076 U.col(k) = VectorType::Unit(
n+1, k);
1077 if (m_compV)
V.col(k) = VectorType::Unit(
n, k);
1085 U(
i,k) = zhat(
i)/(((
diag(
i) - shifts(k)) - mus(k)) )/( (
diag(
i) + singVals[k]));
1088 U.col(k).normalize();
1096 V(
i,k) =
diag(
i) * zhat(
i) / (((
diag(
i) - shifts(k)) - mus(k)) )/( (
diag(
i) + singVals[k]));
1099 V.col(k).normalize();
1103 U.col(
n) = VectorType::Unit(
n+1,
n);
1110 template <
typename MatrixType>
1116 Index start = firstCol + shift;
1122 m_computed(start+
i, start+
i) =
Literal(0);
1125 m_computed(start,start) = r;
1126 m_computed(start+
i, start) =
Literal(0);
1127 m_computed(start+
i, start+
i) =
Literal(0);
1130 if (m_compU) m_naiveU.middleRows(firstCol,
size+1).applyOnTheRight(firstCol, firstCol+
i,
J);
1131 else m_naiveU.applyOnTheRight(firstCol, firstCol+
i,
J);
1139 template <
typename MatrixType>
1149 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1150 std::cout <<
"deflation 4.4: " <<
i <<
"," <<
j <<
" -> " <<
c <<
" " <<
s <<
" " << r <<
" ; "
1151 << m_computed(firstColm +
i-1, firstColm) <<
" "
1152 << m_computed(firstColm +
i, firstColm) <<
" "
1153 << m_computed(firstColm +
i+1, firstColm) <<
" "
1154 << m_computed(firstColm +
i+2, firstColm) <<
"\n";
1155 std::cout << m_computed(firstColm +
i-1, firstColm +
i-1) <<
" "
1156 << m_computed(firstColm +
i, firstColm+
i) <<
" "
1157 << m_computed(firstColm +
i+1, firstColm+
i+1) <<
" "
1158 << m_computed(firstColm +
i+2, firstColm+
i+2) <<
"\n";
1162 m_computed(firstColm +
i, firstColm +
i) = m_computed(firstColm +
j, firstColm +
j);
1167 m_computed(firstColm +
i, firstColm) = r;
1168 m_computed(firstColm +
j, firstColm +
j) = m_computed(firstColm +
i, firstColm +
i);
1169 m_computed(firstColm +
j, firstColm) =
Literal(0);
1172 if (m_compU) m_naiveU.middleRows(firstColu,
size+1).applyOnTheRight(firstColu +
i, firstColu +
j,
J);
1173 else m_naiveU.applyOnTheRight(firstColu+
i, firstColu+
j,
J);
1174 if (m_compV) m_naiveV.middleRows(firstRowW,
size).applyOnTheRight(firstColW +
i, firstColW +
j,
J);
1179 template <
typename MatrixType>
1184 const Index length = lastCol + 1 - firstCol;
1195 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1196 assert(m_naiveU.allFinite());
1197 assert(m_naiveV.allFinite());
1198 assert(m_computed.allFinite());
1201 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1202 std::cout <<
"\ndeflate:" <<
diag.head(k+1).transpose() <<
" | " <<
diag.segment(k+1,length-k-1).transpose() <<
"\n";
1206 if (
diag(0) < epsilon_coarse)
1208 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1209 std::cout <<
"deflation 4.1, because " <<
diag(0) <<
" < " << epsilon_coarse <<
"\n";
1211 diag(0) = epsilon_coarse;
1216 if (
abs(col0(
i)) < epsilon_strict)
1218 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1219 std::cout <<
"deflation 4.2, set z(" <<
i <<
") to zero because " <<
abs(col0(
i)) <<
" < " << epsilon_strict <<
" (diag(" <<
i <<
")=" <<
diag(
i) <<
")\n";
1226 if (
diag(
i) < epsilon_coarse)
1228 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1229 std::cout <<
"deflation 4.3, cancel z(" <<
i <<
")=" << col0(
i) <<
" because diag(" <<
i <<
")=" <<
diag(
i) <<
" < " << epsilon_coarse <<
"\n";
1231 deflation43(firstCol, shift,
i, length);
1234 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1235 assert(m_naiveU.allFinite());
1236 assert(m_naiveV.allFinite());
1237 assert(m_computed.allFinite());
1239 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1240 std::cout <<
"to be sorted: " <<
diag.transpose() <<
"\n\n";
1241 std::cout <<
" : " << col0.transpose() <<
"\n\n";
1246 bool total_deflation = (col0.tail(length-1).array()<considerZero).
all();
1250 Index *permutation = m_workspaceI.data();
1258 permutation[
p++] =
i;
1261 for( ;
p < length; ++
p)
1263 if (
i > k) permutation[
p] =
j++;
1264 else if (
j >= length) permutation[
p] =
i++;
1266 else permutation[
p] =
i++;
1275 Index pi = permutation[
i];
1277 permutation[
i-1] = permutation[
i];
1280 permutation[
i-1] = 0;
1287 Index *realInd = m_workspaceI.data()+length;
1288 Index *realCol = m_workspaceI.data()+2*length;
1296 for(
Index i = total_deflation?0:1;
i < length;
i++)
1298 const Index pi = permutation[length - (total_deflation ?
i+1 :
i)];
1299 const Index J = realCol[pi];
1304 if(
i!=0 &&
J!=0)
swap(col0(
i), col0(
J));
1307 if (m_compU) m_naiveU.col(firstCol+
i).segment(firstCol, length + 1).swap(m_naiveU.col(firstCol+
J).segment(firstCol, length + 1));
1308 else m_naiveU.col(firstCol+
i).segment(0, 2) .swap(m_naiveU.col(firstCol+
J).segment(0, 2));
1309 if (m_compV) m_naiveV.col(firstColW +
i).segment(firstRowW, length).swap(m_naiveV.col(firstColW +
J).segment(firstRowW, length));
1312 const Index realI = realInd[
i];
1319 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1320 std::cout <<
"sorted: " <<
diag.transpose().format(bdcsvdfmt) <<
"\n";
1321 std::cout <<
" : " << col0.transpose() <<
"\n\n";
1327 while(
i>0 && (
abs(
diag(
i))<considerZero ||
abs(col0(
i))<considerZero)) --
i;
1331 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1335 deflation44(firstCol, firstCol + shift, firstRowW, firstColW,
i-1,
i, length);
1339 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1344 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
1345 assert(m_naiveU.allFinite());
1346 assert(m_naiveV.allFinite());
1347 assert(m_computed.allFinite());
1357 template<
typename Derived>