373 for(
Index j=0; j<n; ++j)
375 if( (A.col(j).head(n1).array()!=Literal(0)).any() )
377 A1.col(k1) = A.col(j).head(n1);
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);
394 Map<MatrixXr,Aligned> tmp(m_workspace.data(),n,n);
410template<
typename MatrixType>
417 const Index n = lastCol - firstCol + 1;
419 const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
423 RealScalar lambda, phi, c0, s0;
430 JacobiSVD<MatrixXr> b(m_computed.block(firstCol, firstCol, n + 1, n), ComputeFullU | (m_compV ? ComputeFullV : 0));
432 if (m_info != Success && m_info != NoConvergence)
return;
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);
452 if (m_info != Success && m_info != NoConvergence)
return;
453 divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1);
454 if (m_info != Success && m_info != NoConvergence)
return;
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);
466 r0 = sqrt((abs(alphaK * lambda) * abs(alphaK * lambda)) + abs(betaK * phi) * abs(betaK * phi));
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;
512 RealScalar q1 = m_naiveU(0, firstCol + k);
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()>1e-12*tmp2.abs()).transpose() <<
"\n";
548 static int count = 0;
549 std::cout <<
"# " << ++count <<
"\n\n";
550 assert((tmp1-tmp2).matrix().norm() < 1e-14*tmp2.matrix().norm());
556 MatrixXr UofSVD, VofSVD;
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);
569 Map<Matrix<RealScalar,2,Dynamic>,
Aligned> tmp(m_workspace.data(),2,n+1);
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;
594template <
typename MatrixType>
595void BDCSVD<MatrixType>::computeSVDofM(
Eigen::Index firstCol,
Eigen::Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V)
597 const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
599 ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n);
600 m_workspace.head(n) = m_computed.block(firstCol, firstCol, n, n).diagonal();
601 ArrayRef diag = m_workspace.head(n);
602 diag(0) = Literal(0);
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";
618 while(actual_n>1 && diag(actual_n-1)==Literal(0)) {--actual_n; eigen_internal_assert(col0(actual_n)==Literal(0)); }
620 for(Index k=0;k<actual_n;++k)
621 if(abs(col0(k))>considerZero)
622 m_workspaceI(m++) = k;
623 Map<ArrayXi> perm(m_workspaceI.data(),m);
625 Map<ArrayXr> shifts(m_workspace.data()+1*n, n);
626 Map<ArrayXr> mus(m_workspace.data()+2*n, n);
627 Map<ArrayXr> zhat(m_workspace.data()+3*n, n);
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());
688 for(Index i=0; i<actual_n-1; ++i)
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
715 JacobiSVD<MatrixXr> jsvd(m_computed.block(firstCol, firstCol, n, n) );
716 std::cout <<
" * j: " << jsvd.singularValues().transpose() <<
"\n\n";
717 std::cout <<
" * sing-val: " << singVals.transpose() <<
"\n";
722template <
typename MatrixType>
723typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar mu,
const ArrayRef& col0,
const ArrayRef& diag,
const IndicesRef &perm,
const ArrayRef& diagShifted, RealScalar shift)
725 Index m = perm.size();
726 RealScalar res = Literal(1);
727 for(Index i=0; i<m; ++i)
732 res += (col0(j) / (diagShifted(j) - mu)) * (col0(j) / (diag(j) + shift + mu));
738template <
typename MatrixType>
739void BDCSVD<MatrixType>::computeSingVals(
const ArrayRef& col0,
const ArrayRef& diag,
const IndicesRef &perm,
740 VectorType& singVals, ArrayRef shifts, ArrayRef mus)
746 Index n = col0.size();
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);
765 RealScalar left = diag(k);
768 right = (diag(actual_n-1) + col0.matrix().norm());
775 while(col0(l)==Literal(0)) { ++l; eigen_internal_assert(l<actual_n); }
780 RealScalar mid = left + (right-left) / Literal(2);
781 RealScalar fMid = secularEq(mid, col0, diag, perm, diag, Literal(0));
782#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
783 std::cout <<
"right-left = " << right-left <<
"\n";
786 std::cout <<
" = " << secularEq(left+RealScalar(0.000001)*(right-left), col0, diag, perm, diag, 0)
787 <<
" " << secularEq(left+RealScalar(0.1) *(right-left), col0, diag, perm, diag, 0)
788 <<
" " << secularEq(left+RealScalar(0.2) *(right-left), col0, diag, perm, diag, 0)
789 <<
" " << secularEq(left+RealScalar(0.3) *(right-left), col0, diag, perm, diag, 0)
790 <<
" " << secularEq(left+RealScalar(0.4) *(right-left), col0, diag, perm, diag, 0)
791 <<
" " << secularEq(left+RealScalar(0.49) *(right-left), col0, diag, perm, diag, 0)
792 <<
" " << secularEq(left+RealScalar(0.5) *(right-left), col0, diag, perm, diag, 0)
793 <<
" " << secularEq(left+RealScalar(0.51) *(right-left), col0, diag, perm, diag, 0)
794 <<
" " << secularEq(left+RealScalar(0.6) *(right-left), col0, diag, perm, diag, 0)
795 <<
" " << secularEq(left+RealScalar(0.7) *(right-left), col0, diag, perm, diag, 0)
796 <<
" " << secularEq(left+RealScalar(0.8) *(right-left), col0, diag, perm, diag, 0)
797 <<
" " << secularEq(left+RealScalar(0.9) *(right-left), col0, diag, perm, diag, 0)
798 <<
" " << secularEq(left+RealScalar(0.999999)*(right-left), col0, diag, perm, diag, 0) <<
"\n";
800 RealScalar shift = (k == actual_n-1 || fMid > Literal(0)) ? left : right;
803 Map<ArrayXr> diagShifted(m_workspace.data()+4*n, n);
804 diagShifted = diag - shift;
809 RealScalar midShifted = (right - left) / RealScalar(2);
811 midShifted = -midShifted;
812 RealScalar fMidShifted = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
816 shift = fMidShifted > Literal(0) ? left : right;
817 diagShifted = diag - shift;
822 RealScalar muPrev, muCur;
825 muPrev = (right - left) * RealScalar(0.1);
826 if (k == actual_n-1) muCur = right - left;
827 else muCur = (right - left) * RealScalar(0.5);
831 muPrev = -(right - left) * RealScalar(0.1);
832 muCur = -(right - left) * RealScalar(0.5);
835 RealScalar fPrev = secularEq(muPrev, col0, diag, perm, diagShifted, shift);
836 RealScalar fCur = secularEq(muCur, col0, diag, perm, diagShifted, shift);
837 if (abs(fPrev) < abs(fCur))
845 bool useBisection = fPrev*fCur>Literal(0);
846 while (fCur!=Literal(0) && abs(muCur - muPrev) > Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
851 RealScalar a = (fCur - fPrev) / (Literal(1)/muCur - Literal(1)/muPrev);
852 RealScalar b = fCur - a / muCur;
854 RealScalar muZero = -a/b;
855 RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift);
857#ifdef EIGEN_BDCSVD_SANITY_CHECKS
858 assert((numext::isfinite)(fZero));
866 if (shift == left && (muCur < Literal(0) || muCur > right - left)) useBisection =
true;
867 if (shift == right && (muCur < -(right - left) || muCur > Literal(0))) useBisection =
true;
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";
877 RealScalar leftShifted, rightShifted;
882 leftShifted = numext::maxi<RealScalar>( (std::numeric_limits<RealScalar>::min)(), Literal(2) * abs(col0(k)) / sqrt((std::numeric_limits<RealScalar>::max)()) );
885 eigen_internal_assert( (numext::isfinite)( (col0(k)/leftShifted)*(col0(k)/(diag(k)+shift+leftShifted)) ) );
888 rightShifted = (k==actual_n-1) ? right : ((right - left) * RealScalar(0.51));
892 leftShifted = -(right - left) * RealScalar(0.51);
894 rightShifted = -numext::maxi<RealScalar>( (std::numeric_limits<RealScalar>::min)(), abs(col0(k+1)) / sqrt((std::numeric_limits<RealScalar>::max)()) );
896 rightShifted = -(std::numeric_limits<RealScalar>::min)();
899 RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
900 eigen_internal_assert(fLeft<Literal(0));
902#if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_SANITY_CHECKS
903 RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift);
906#ifdef EIGEN_BDCSVD_SANITY_CHECKS
907 if(!(numext::isfinite)(fLeft))
908 std::cout <<
"f(" << leftShifted <<
") =" << fLeft <<
" ; " << left <<
" " << shift <<
" " << right <<
"\n";
909 assert((numext::isfinite)(fLeft));
911 if(!(numext::isfinite)(fRight))
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";
927 eigen_internal_assert(fLeft * fRight < Literal(0));
931 while (rightShifted - leftShifted > Literal(2) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(leftShifted), abs(rightShifted)))
933 RealScalar midShifted = (leftShifted + rightShifted) / Literal(2);
934 fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
935 eigen_internal_assert((numext::isfinite)(fMid));
937 if (fLeft * fMid < Literal(0))
939 rightShifted = midShifted;
943 leftShifted = midShifted;
947 muCur = (leftShifted + rightShifted) / Literal(2);
955 muCur = (right - left) * RealScalar(0.5);
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));
984template <
typename MatrixType>
985void BDCSVD<MatrixType>::perturbCol0
986 (
const ArrayRef& col0,
const ArrayRef& diag,
const IndicesRef &perm,
const VectorType& singVals,
987 const ArrayRef& shifts,
const ArrayRef& mus, ArrayRef zhat)
990 Index n = col0.size();
991 Index m = perm.size();
997 Index lastIdx = perm(m-1);
999 for (Index k = 0; k < n; ++k)
1001 if (col0(k) == Literal(0))
1002 zhat(k) = Literal(0);
1006 RealScalar dk = diag(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";
1017 for(Index l = 0; l<m; ++l)
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";
1028 Index j = (i<k ) ? i : perm(l-1);
1029 std::cout <<
" " <<
"j=" << j <<
"\n";
1032 Index j = i<k ? i : perm(l-1);
1033#ifdef EIGEN_BDCSVD_SANITY_CHECKS
1034 if(!(dk!=Literal(0) || diag(i)!=Literal(0)))
1036 std::cout <<
"k=" << k <<
", i=" << i <<
", l=" << l <<
", perm.size()=" << perm.size() <<
"\n";
1038 assert(dk!=Literal(0) || diag(i)!=Literal(0));
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
1045 if(i!=k && numext::abs(((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) - 1) > 0.9 )
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";
1054 RealScalar tmp = sqrt(prod);
1055#ifdef EIGEN_BDCSVD_SANITY_CHECKS
1056 assert((numext::isfinite)(tmp));
1058 zhat(k) = col0(k) > Literal(0) ? RealScalar(tmp) : RealScalar(-tmp);
1064template <
typename MatrixType>
1065void BDCSVD<MatrixType>::computeSingVecs
1066 (
const ArrayRef& zhat,
const ArrayRef& diag,
const IndicesRef &perm,
const VectorType& singVals,
1067 const ArrayRef& shifts,
const ArrayRef& mus, MatrixXr& U, MatrixXr& V)
1069 Index n = zhat.size();
1070 Index m = perm.size();
1072 for (Index k = 0; k < n; ++k)
1074 if (zhat(k) == Literal(0))
1076 U.col(k) = VectorType::Unit(n+1, k);
1077 if (m_compV) V.col(k) = VectorType::Unit(n, k);
1082 for(Index l=0;l<m;++l)
1085 U(i,k) = zhat(i)/(((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
1087 U(n,k) = Literal(0);
1088 U.col(k).normalize();
1093 for(Index l=1;l<m;++l)
1096 V(i,k) = diag(i) * zhat(i) / (((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
1098 V(0,k) = Literal(-1);
1099 V.col(k).normalize();
1103 U.col(n) = VectorType::Unit(n+1, n);
1110template <
typename MatrixType>
1116 Index start = firstCol + shift;
1117 RealScalar c = m_computed(start, start);
1118 RealScalar s = m_computed(start+i, start);
1119 RealScalar r = numext::hypot(c,s);
1120 if (r == Literal(0))
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);
1129 JacobiRotation<RealScalar> J(c/r,-s/r);
1130 if (m_compU) m_naiveU.middleRows(firstCol, size+1).applyOnTheRight(firstCol, firstCol+i, J);
1131 else m_naiveU.applyOnTheRight(firstCol, firstCol+i, J);
1139template <
typename MatrixType>
1146 RealScalar c = m_computed(firstColm+i, firstColm);
1147 RealScalar s = m_computed(firstColm+j, firstColm);
1148 RealScalar r = sqrt(numext::abs2(c) + numext::abs2(s));
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);
1171 JacobiRotation<RealScalar> J(c,-s);
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);
1179template <
typename MatrixType>
1184 const Index length = lastCol + 1 - firstCol;
1186 Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1);
1187 Diagonal<MatrixXr> fulldiag(m_computed);
1188 VectorBlock<Diagonal<MatrixXr>,
Dynamic> diag(fulldiag, firstCol+shift, length);
1190 const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
1191 RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff();
1192 RealScalar epsilon_strict = numext::maxi<RealScalar>(considerZero,NumTraits<RealScalar>::epsilon() * maxDiag);
1193 RealScalar epsilon_coarse = Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag);
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;
1215 for (Index i=1;i<length;++i)
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";
1221 col0(i) = Literal(0);
1225 for (Index i=1;i<length; i++)
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();
1256 for(Index i=1; i<length; ++i)
1257 if(abs(diag(i))<considerZero)
1258 permutation[p++] = i;
1261 for( ; p < length; ++p)
1263 if (i > k) permutation[p] = j++;
1264 else if (j >= length) permutation[p] = i++;
1265 else if (diag(i) < diag(j)) permutation[p] = j++;
1266 else permutation[p] = i++;
1273 for(Index i=1; i<length; ++i)
1275 Index pi = permutation[i];
1276 if(abs(diag(pi))<considerZero || diag(0)<diag(pi))
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;
1290 for(
int pos = 0; pos< length; pos++)
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];
1303 swap(diag(i), diag(J));
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;
1329 if( (diag(i) - diag(i-1)) < NumTraits<RealScalar>::epsilon()*maxDiag )
1331#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
1332 std::cout <<
"deflation 4.4 with i = " << i <<
" because " << diag(i) <<
" - " << diag(i-1) <<
" == " << (diag(i) - diag(i-1)) <<
" < " << NumTraits<RealScalar>::epsilon()*maxDiag <<
"\n";
1334 eigen_internal_assert(abs(diag(i) - diag(i-1))<epsilon_coarse &&
" diagonal entries are not properly sorted");
1335 deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i-1, i, length);
1339#ifdef EIGEN_BDCSVD_SANITY_CHECKS
1340 for(Index j=2;j<length;++j)
1341 assert(diag(j-1)<=diag(j) || abs(diag(j))<considerZero);
1344#ifdef EIGEN_BDCSVD_SANITY_CHECKS
1345 assert(m_naiveU.allFinite());
1346 assert(m_naiveV.allFinite());
1347 assert(m_computed.allFinite());
1357template<
typename Derived>
1358BDCSVD<typename MatrixBase<Derived>::PlainObject>