class Bidiagonal Divide and Conquer SVD More...

#include <BDCSVD.h>

Inheritance diagram for Eigen::BDCSVD< MatrixType_, Options_ >:

Public Types
enum	{ Options = Options_ , QRDecomposition = Options & internal::QRPreconditionerBits , ComputationOptions = Options & internal::ComputationOptionsBits , RowsAtCompileTime = Base::RowsAtCompileTime , ColsAtCompileTime = Base::ColsAtCompileTime , DiagSizeAtCompileTime = Base::DiagSizeAtCompileTime , MaxRowsAtCompileTime = Base::MaxRowsAtCompileTime , MaxColsAtCompileTime = Base::MaxColsAtCompileTime , MaxDiagSizeAtCompileTime = Base::MaxDiagSizeAtCompileTime , MatrixOptions = Base::MatrixOptions }

typedef MatrixType_	MatrixType

typedef Base::Scalar	Scalar

typedef Base::RealScalar	RealScalar

typedef NumTraits< RealScalar >::Literal	Literal

typedef Base::Index	Index

typedef Base::MatrixUType	MatrixUType

typedef Base::MatrixVType	MatrixVType

typedef Base::SingularValuesType	SingularValuesType

typedef Matrix< Scalar, Dynamic, Dynamic, ColMajor >	MatrixX

typedef Matrix< RealScalar, Dynamic, Dynamic, ColMajor >	MatrixXr

typedef Matrix< RealScalar, Dynamic, 1 >	VectorType

typedef Array< RealScalar, Dynamic, 1 >	ArrayXr

typedef Array< Index, 1, Dynamic >	ArrayXi

typedef Ref< ArrayXr >	ArrayRef

typedef Ref< ArrayXi >	IndicesRef

Public Types inherited from Eigen::SVDBase< BDCSVD< MatrixType_, Options_ > >
enum

typedef internal::traits< BDCSVD< MatrixType_, Options_ > >::MatrixType	MatrixType

typedef MatrixType::Scalar	Scalar

typedef NumTraits< typename MatrixType::Scalar >::Real	RealScalar

typedef Eigen::internal::traits< SVDBase >::StorageIndex	StorageIndex

typedef internal::make_proper_matrix_type< Scalar, RowsAtCompileTime, MatrixUColsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MatrixUMaxColsAtCompileTime >::type	MatrixUType

typedef internal::make_proper_matrix_type< Scalar, ColsAtCompileTime, MatrixVColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MatrixVMaxColsAtCompileTime >::type	MatrixVType

typedef internal::plain_diag_type< MatrixType, RealScalar >::type	SingularValuesType

Public Types inherited from Eigen::SolverBase< Derived >
enum	{ RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime , ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime , SizeAtCompileTime = (internal::size_of_xpr_at_compile_time<Derived>::ret) , MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime , MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime , MaxSizeAtCompileTime , IsVectorAtCompileTime , NumDimensions }

typedef EigenBase< Derived >	Base

typedef internal::traits< Derived >::Scalar	Scalar

typedef Scalar	CoeffReturnType

typedef Transpose< const Derived >	ConstTransposeReturnType

typedef std::conditional_t< NumTraits< Scalar >::IsComplex, CwiseUnaryOp< internal::scalar_conjugate_op< Scalar >, const ConstTransposeReturnType >, const ConstTransposeReturnType >	AdjointReturnType

Public Types inherited from Eigen::EigenBase< Derived >
typedef Eigen::Index	Index
	The interface type of indices. More...

typedef internal::traits< Derived >::StorageKind	StorageKind

Public Member Functions
	BDCSVD ()
	Default Constructor. More...

	BDCSVD (Index rows, Index cols)
	Default Constructor with memory preallocation. More...

EIGEN_DEPRECATED	BDCSVD (Index rows, Index cols, unsigned int computationOptions)
	Default Constructor with memory preallocation. More...

	BDCSVD (const MatrixType &matrix)
	Constructor performing the decomposition of given matrix, using the custom options specified with the Options template parameter. More...

EIGEN_DEPRECATED	BDCSVD (const MatrixType &matrix, unsigned int computationOptions)
	Constructor performing the decomposition of given matrix using specified options for computing unitaries. More...

	~BDCSVD ()

BDCSVD &	compute (const MatrixType &matrix)
	Method performing the decomposition of given matrix. Computes Thin/Full unitaries U/V if specified using the Options template parameter or the class constructor. More...

EIGEN_DEPRECATED BDCSVD &	compute (const MatrixType &matrix, unsigned int computationOptions)
	Method performing the decomposition of given matrix, as specified by the `computationOptions` parameter. More...

void	setSwitchSize (int s)

EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index	cols () const EIGEN_NOEXCEPT

EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index	rows () const EIGEN_NOEXCEPT

Public Member Functions inherited from Eigen::SVDBase< BDCSVD< MatrixType_, Options_ > >
BDCSVD< MatrixType_, Options_ > &	derived ()

const BDCSVD< MatrixType_, Options_ > &	derived () const

const MatrixUType &	matrixU () const

const MatrixVType &	matrixV () const

const SingularValuesType &	singularValues () const

Index	nonzeroSingularValues () const

Index	rank () const

BDCSVD< MatrixType_, Options_ > &	setThreshold (const RealScalar &threshold)

BDCSVD< MatrixType_, Options_ > &	setThreshold (Default_t)

RealScalar	threshold () const

bool	computeU () const

bool	computeV () const

Index	rows () const

Index	cols () const

Index	diagSize () const

EIGEN_DEVICE_FUNC ComputationInfo	info () const
	Reports whether previous computation was successful. More...

void	_solve_impl (const RhsType &rhs, DstType &dst) const

void	_solve_impl_transposed (const RhsType &rhs, DstType &dst) const

Public Member Functions inherited from Eigen::SolverBase< Derived >
	SolverBase ()

	~SolverBase ()

template<typename Rhs >
const Solve< Derived, Rhs >	solve (const MatrixBase< Rhs > &b) const

const ConstTransposeReturnType	transpose () const

const AdjointReturnType	adjoint () const

constexpr EIGEN_DEVICE_FUNC Derived &	derived ()

constexpr EIGEN_DEVICE_FUNC const Derived &	derived () const

Public Member Functions inherited from Eigen::EigenBase< Derived >
constexpr EIGEN_DEVICE_FUNC Derived &	derived ()

constexpr EIGEN_DEVICE_FUNC const Derived &	derived () const

EIGEN_DEVICE_FUNC Derived &	const_cast_derived () const

EIGEN_DEVICE_FUNC const Derived &	const_derived () const

EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index	rows () const EIGEN_NOEXCEPT

EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index	cols () const EIGEN_NOEXCEPT

EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index	size () const EIGEN_NOEXCEPT

template<typename Dest >
EIGEN_DEVICE_FUNC void	evalTo (Dest &dst) const

template<typename Dest >
EIGEN_DEVICE_FUNC void	addTo (Dest &dst) const

template<typename Dest >
EIGEN_DEVICE_FUNC void	subTo (Dest &dst) const

template<typename Dest >
EIGEN_DEVICE_FUNC void	applyThisOnTheRight (Dest &dst) const

template<typename Dest >
EIGEN_DEVICE_FUNC void	applyThisOnTheLeft (Dest &dst) const

template<typename Device >
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE DeviceWrapper< Derived, Device >	device (Device &device)

template<typename Device >
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE DeviceWrapper< const Derived, Device >	device (Device &device) const

Public Attributes
int	m_numIters

Protected Member Functions
void	allocate (Index rows, Index cols, unsigned int computationOptions)

Protected Member Functions inherited from Eigen::SVDBase< BDCSVD< MatrixType_, Options_ > >
void	_check_compute_assertions () const

void	_check_solve_assertion (const Rhs &b) const

bool	allocate (Index rows, Index cols, unsigned int computationOptions)

	SVDBase ()
	Default Constructor. More...

Protected Member Functions inherited from Eigen::SolverBase< Derived >
template<bool Transpose_, typename Rhs >
void	_check_solve_assertion (const Rhs &b) const

Protected Attributes
MatrixXr	m_naiveU

MatrixXr	m_naiveV

MatrixXr	m_computed

Index	m_nRec

ArrayXr	m_workspace

ArrayXi	m_workspaceI

int	m_algoswap

bool	m_isTranspose

bool	m_compU

bool	m_compV

bool	m_useQrDecomp

JacobiSVD< MatrixType, ComputationOptions >	smallSvd

HouseholderQR< MatrixX >	qrDecomp

internal::UpperBidiagonalization< MatrixX >	bid

MatrixX	copyWorkspace

MatrixX	reducedTriangle

Protected Attributes inherited from Eigen::SVDBase< BDCSVD< MatrixType_, Options_ > >
MatrixUType	m_matrixU

MatrixVType	m_matrixV

SingularValuesType	m_singularValues

ComputationInfo	m_info

bool	m_isInitialized

bool	m_isAllocated

bool	m_usePrescribedThreshold

bool	m_computeFullU

bool	m_computeThinU

bool	m_computeFullV

bool	m_computeThinV

unsigned int	m_computationOptions

Index	m_nonzeroSingularValues

internal::variable_if_dynamic< Index, RowsAtCompileTime >	m_rows

internal::variable_if_dynamic< Index, ColsAtCompileTime >	m_cols

internal::variable_if_dynamic< Index, DiagSizeAtCompileTime >	m_diagSize

RealScalar	m_prescribedThreshold

Private Types
typedef SVDBase< BDCSVD >	Base

Private Member Functions
BDCSVD &	compute_impl (const MatrixType &matrix, 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 ArrayRef &col0, const ArrayRef &diag, const IndicesRef &perm, VectorType &singVals, ArrayRef shifts, ArrayRef mus)

void	perturbCol0 (const ArrayRef &col0, const ArrayRef &diag, const IndicesRef &perm, const VectorType &singVals, const ArrayRef &shifts, const ArrayRef &mus, ArrayRef zhat)

void	computeSingVecs (const ArrayRef &zhat, const ArrayRef &diag, const IndicesRef &perm, const VectorType &singVals, const ArrayRef &shifts, const ArrayRef &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)

void	structured_update (Block< MatrixXr, Dynamic, Dynamic > A, const MatrixXr &B, Index n1)

template<typename SVDType >
void	computeBaseCase (SVDType &svd, Index n, Index firstCol, Index firstRowW, Index firstColW, Index shift)

Static Private Member Functions
static RealScalar	secularEq (RealScalar x, const ArrayRef &col0, const ArrayRef &diag, const IndicesRef &perm, const ArrayRef &diagShifted, RealScalar shift)

Additional Inherited Members
Static Public Attributes inherited from Eigen::SVDBase< BDCSVD< MatrixType_, Options_ > >
static constexpr bool	ShouldComputeFullU

static constexpr bool	ShouldComputeThinU

static constexpr bool	ShouldComputeFullV

static constexpr bool	ShouldComputeThinV

Detailed Description

template<typename MatrixType_, int Options_>
class Eigen::BDCSVD< MatrixType_, Options_ >

class Bidiagonal Divide and Conquer SVD

Template Parameters

MatrixType_ the type of the matrix of which we are computing the SVD decomposition

Options_ this optional parameter allows one to specify options for computing unitaries U and V. Possible values are ComputeThinU, ComputeThinV, ComputeFullU, ComputeFullV, and DisableQRDecomposition. It is not possible to request both the thin and full version of U or V. By default, unitaries are not computed. BDCSVD uses R-Bidiagonalization to improve performance on tall and wide matrices. For backwards compatility, the option DisableQRDecomposition can be used to disable this optimization.

This class first reduces the input matrix to bi-diagonal form using class UpperBidiagonalization, and then performs a divide-and-conquer diagonalization. Small blocks are diagonalized using class JacobiSVD. You can control the switching size with the setSwitchSize() method, default is 16. For small matrice (<16), it is thus preferable to directly use JacobiSVD. For larger ones, BDCSVD is highly recommended and can several order of magnitude faster.

Warning: this algorithm is unlikely to provide accurate result when compiled with unsafe math optimizations. For instance, this concerns Intel's compiler (ICC), which performs such optimization by default unless you compile with the -fp-model precise option. Likewise, the -ffast-math option of GCC or clang will significantly degrade the accuracy.

See also: class JacobiSVD

Member Typedef Documentation

◆ ArrayRef

template<typename MatrixType_ , int Options_>

typedef Ref<ArrayXr> Eigen::BDCSVD< MatrixType_, Options_ >::ArrayRef

◆ ArrayXi

template<typename MatrixType_ , int Options_>

typedef Array<Index, 1, Dynamic> Eigen::BDCSVD< MatrixType_, Options_ >::ArrayXi

◆ ArrayXr

template<typename MatrixType_ , int Options_>

typedef Array<RealScalar, Dynamic, 1> Eigen::BDCSVD< MatrixType_, Options_ >::ArrayXr

◆ Base

template<typename MatrixType_ , int Options_>

typedef SVDBase<BDCSVD> Eigen::BDCSVD< MatrixType_, Options_ >::Base

private

◆ Index

template<typename MatrixType_ , int Options_>

typedef Base::Index Eigen::BDCSVD< MatrixType_, Options_ >::Index

◆ IndicesRef

template<typename MatrixType_ , int Options_>

typedef Ref<ArrayXi> Eigen::BDCSVD< MatrixType_, Options_ >::IndicesRef

◆ Literal

template<typename MatrixType_ , int Options_>

typedef NumTraits<RealScalar>::Literal Eigen::BDCSVD< MatrixType_, Options_ >::Literal

◆ MatrixType

template<typename MatrixType_ , int Options_>

typedef MatrixType_ Eigen::BDCSVD< MatrixType_, Options_ >::MatrixType

◆ MatrixUType

template<typename MatrixType_ , int Options_>

typedef Base::MatrixUType Eigen::BDCSVD< MatrixType_, Options_ >::MatrixUType

◆ MatrixVType

template<typename MatrixType_ , int Options_>

typedef Base::MatrixVType Eigen::BDCSVD< MatrixType_, Options_ >::MatrixVType

◆ MatrixX

template<typename MatrixType_ , int Options_>

typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> Eigen::BDCSVD< MatrixType_, Options_ >::MatrixX

◆ MatrixXr

template<typename MatrixType_ , int Options_>

typedef Matrix<RealScalar, Dynamic, Dynamic, ColMajor> Eigen::BDCSVD< MatrixType_, Options_ >::MatrixXr

◆ RealScalar

template<typename MatrixType_ , int Options_>

typedef Base::RealScalar Eigen::BDCSVD< MatrixType_, Options_ >::RealScalar

◆ Scalar

template<typename MatrixType_ , int Options_>

typedef Base::Scalar Eigen::BDCSVD< MatrixType_, Options_ >::Scalar

◆ SingularValuesType

template<typename MatrixType_ , int Options_>

typedef Base::SingularValuesType Eigen::BDCSVD< MatrixType_, Options_ >::SingularValuesType

◆ VectorType

template<typename MatrixType_ , int Options_>

typedef Matrix<RealScalar, Dynamic, 1> Eigen::BDCSVD< MatrixType_, Options_ >::VectorType

Member Enumeration Documentation

◆ anonymous enum

template<typename MatrixType_ , int Options_>

anonymous enum

Enumerator
Options
QRDecomposition
ComputationOptions
RowsAtCompileTime
ColsAtCompileTime
DiagSizeAtCompileTime
MaxRowsAtCompileTime
MaxColsAtCompileTime
MaxDiagSizeAtCompileTime
MatrixOptions

        {
     Options = Options_,
     QRDecomposition = Options & internal::QRPreconditionerBits,
     ComputationOptions = Options & internal::ComputationOptionsBits,
     RowsAtCompileTime = Base::RowsAtCompileTime,
     ColsAtCompileTime = Base::ColsAtCompileTime,
     DiagSizeAtCompileTime = Base::DiagSizeAtCompileTime,
     MaxRowsAtCompileTime = Base::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = Base::MaxColsAtCompileTime,
     MaxDiagSizeAtCompileTime = Base::MaxDiagSizeAtCompileTime,
     MatrixOptions = Base::MatrixOptions
   };

Constructor & Destructor Documentation

◆ BDCSVD() [1/5]

template<typename MatrixType_ , int Options_>

Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD ( )

inline

Default Constructor.

The default constructor is useful in cases in which the user intends to perform decompositions via BDCSVD::compute(const MatrixType&).

130 : m_algoswap(16), m_isTranspose(false), m_compU(false), m_compV(false), m_numIters(0) {}

Eigen::BDCSVD::m_isTranspose

bool m_isTranspose

Definition: BDCSVD.h:245

Eigen::BDCSVD::m_numIters

int m_numIters

Definition: BDCSVD.h:263

Eigen::BDCSVD::m_compU

bool m_compU

Definition: BDCSVD.h:245

Eigen::BDCSVD::m_algoswap

int m_algoswap

Definition: BDCSVD.h:244

Eigen::BDCSVD::m_compV

bool m_compV

Definition: BDCSVD.h:245

◆ BDCSVD() [2/5]

template<typename MatrixType_ , int Options_>

Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD	(	Index	rows,
		Index	cols
	)

inline

Default Constructor with memory preallocation.

Like the default constructor but with preallocation of the internal data according to the specified problem size and Options template parameter.

See also: BDCSVD()

                                  : m_algoswap(16), m_numIters(0) {
     allocate(rows, cols, internal::get_computation_options(Options));
   }

References Eigen::BDCSVD< MatrixType_, Options_ >::allocate(), Eigen::BDCSVD< MatrixType_, Options_ >::cols(), Eigen::internal::get_computation_options(), Eigen::BDCSVD< MatrixType_, Options_ >::Options, and Eigen::BDCSVD< MatrixType_, Options_ >::rows().

◆ BDCSVD() [3/5]

template<typename MatrixType_ , int Options_>

EIGEN_DEPRECATED Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD	(	Index	rows,
		Index	cols,
		unsigned int	computationOptions
	)

inline

Default Constructor with memory preallocation.

Like the default constructor but with preallocation of the internal data according to the specified problem size and the computationOptions.

One cannot request unitaries using both the Options template parameter and the constructor. If possible, prefer using the Options template parameter.

Parameters

computationOptions specification for computing Thin/Full unitaries U/V

See also: BDCSVD()

Deprecated:: Will be removed in the next major Eigen version. Options should be specified in the Options template parameter.

                                                                                    : m_algoswap(16), m_numIters(0) {
     internal::check_svd_options_assertions<MatrixType, Options>(computationOptions, rows, cols);
     allocate(rows, cols, computationOptions);
   }

References Eigen::BDCSVD< MatrixType_, Options_ >::allocate(), Eigen::BDCSVD< MatrixType_, Options_ >::cols(), and Eigen::BDCSVD< MatrixType_, Options_ >::rows().

◆ BDCSVD() [4/5]

template<typename MatrixType_ , int Options_>

Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD ( const MatrixType & matrix )

inline

Constructor performing the decomposition of given matrix, using the custom options specified with the Options template parameter.

Parameters

matrix the matrix to decompose

                                    : m_algoswap(16), m_numIters(0) {
     compute_impl(matrix, internal::get_computation_options(Options));
   }

References Eigen::BDCSVD< MatrixType_, Options_ >::compute_impl(), Eigen::internal::get_computation_options(), matrix(), and Eigen::BDCSVD< MatrixType_, Options_ >::Options.

◆ BDCSVD() [5/5]

template<typename MatrixType_ , int Options_>

EIGEN_DEPRECATED Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD	(	const MatrixType &	matrix,
		unsigned int	computationOptions
	)

inline

Constructor performing the decomposition of given matrix using specified options for computing unitaries.

One cannot request unitaries using both the Options template parameter and the constructor. If possible, prefer using the Options template parameter.

Parameters

matrix	the matrix to decompose
computationOptions	specification for computing Thin/Full unitaries U/V

Deprecated:: Will be removed in the next major Eigen version. Options should be specified in the Options template parameter.

                                                                                      : m_algoswap(16), m_numIters(0) {
     internal::check_svd_options_assertions<MatrixType, Options>(computationOptions, matrix.rows(), matrix.cols());
     compute_impl(matrix, computationOptions);
   }

References Eigen::BDCSVD< MatrixType_, Options_ >::compute_impl(), and matrix().

◆ ~BDCSVD()

template<typename MatrixType_ , int Options_>

Eigen::BDCSVD< MatrixType_, Options_ >::~BDCSVD ( )

inline

187 {}

Member Function Documentation

◆ allocate()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::allocate	(	Index	rows,
		Index	cols,
		unsigned int	computationOptions
	)

protected

                                                                                                   {
   if (Base::allocate(rows, cols, computationOptions)) return;
  
   if (cols < m_algoswap)
     smallSvd.allocate(rows, cols, Options == 0 ? computationOptions : internal::get_computation_options(Options));
  
   m_computed = MatrixXr::Zero(diagSize() + 1, diagSize());
   m_compU = computeV();
   m_compV = computeU();
   m_isTranspose = (cols > rows);
   if (m_isTranspose) std::swap(m_compU, m_compV);
  
   // kMinAspectRatio is the crossover point that determines if we perform R-Bidiagonalization
   // or bidiagonalize the input matrix directly.
   // It is based off of LAPACK's dgesdd routine, which uses 11.0/6.0
   // we use a larger scalar to prevent a regression for relatively square matrices.
   constexpr Index kMinAspectRatio = 4;
   constexpr bool disableQrDecomp = static_cast<int>(QRDecomposition) == static_cast<int>(DisableQRDecomposition);
   m_useQrDecomp = !disableQrDecomp && ((rows / kMinAspectRatio > cols) || (cols / kMinAspectRatio > rows));
   if (m_useQrDecomp) {
     qrDecomp = HouseholderQR<MatrixX>((std::max)(rows, cols), (std::min)(rows, cols));
     reducedTriangle = MatrixX(diagSize(), diagSize());
   }
  
   copyWorkspace = MatrixX(m_isTranspose ? cols : rows, m_isTranspose ? rows : cols);
   bid = internal::UpperBidiagonalization<MatrixX>(m_useQrDecomp ? diagSize() : copyWorkspace.rows(),
                                                   m_useQrDecomp ? diagSize() : copyWorkspace.cols());
  
   if (m_compU)
     m_naiveU = MatrixXr::Zero(diagSize() + 1, diagSize() + 1);
   else
     m_naiveU = MatrixXr::Zero(2, diagSize() + 1);
  
   if (m_compV) m_naiveV = MatrixXr::Zero(diagSize(), diagSize());
  
   m_workspace.resize((diagSize() + 1) * (diagSize() + 1) * 3);
   m_workspaceI.resize(3 * diagSize());
 }  // end allocate

References cols, Eigen::DisableQRDecomposition, Eigen::internal::get_computation_options(), max, min, rows, swap(), and oomph::PseudoSolidHelper::Zero.

Referenced by Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD().

◆ cols()

template<typename MatrixType_ , int Options_>

EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index Eigen::EigenBase< Derived >::cols

inline

Returns: the number of columns.

See also: rows(), ColsAtCompileTime

61 { return derived().cols(); }

Eigen::SVDBase< BDCSVD< MatrixType_, Options_ > >::derived

BDCSVD< MatrixType_, Options_ > & derived()

Definition: SVDBase.h:160

Referenced by gdb.printers._MatrixEntryIterator::__next__(), Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD(), gdb.printers.EigenMatrixPrinter::children(), gdb.printers.EigenSparseMatrixPrinter::children(), gdb.printers.EigenMatrixPrinter::to_string(), and gdb.printers.EigenSparseMatrixPrinter::to_string().

◆ compute() [1/2]

template<typename MatrixType_ , int Options_>

BDCSVD& Eigen::BDCSVD< MatrixType_, Options_ >::compute ( const MatrixType & matrix )

inline

Method performing the decomposition of given matrix. Computes Thin/Full unitaries U/V if specified using the Options template parameter or the class constructor.

Parameters

matrix the matrix to decompose

194 { return compute_impl(matrix, m_computationOptions); }

Eigen::SVDBase< BDCSVD< MatrixType_, Options_ > >::m_computationOptions

unsigned int m_computationOptions

Definition: SVDBase.h:338

References Eigen::BDCSVD< MatrixType_, Options_ >::compute_impl(), Eigen::SVDBase< BDCSVD< MatrixType_, Options_ > >::m_computationOptions, and matrix().

Referenced by bench(), and compare_bdc_jacobi().

◆ compute() [2/2]

template<typename MatrixType_ , int Options_>

EIGEN_DEPRECATED BDCSVD& Eigen::BDCSVD< MatrixType_, Options_ >::compute	(	const MatrixType &	matrix,
		unsigned int	computationOptions
	)

inline

Method performing the decomposition of given matrix, as specified by the computationOptions parameter.

Parameters

matrix	the matrix to decompose
computationOptions	specify whether to compute Thin/Full unitaries U/V

Deprecated:: Will be removed in the next major Eigen version. Options should be specified in the Options template parameter.

                                                                                               {
     internal::check_svd_options_assertions<MatrixType, Options>(computationOptions, matrix.rows(), matrix.cols());
     return compute_impl(matrix, computationOptions);
   }

References Eigen::BDCSVD< MatrixType_, Options_ >::compute_impl(), and matrix().

◆ compute_impl()

template<typename MatrixType , int Options>

BDCSVD< MatrixType, Options > & Eigen::BDCSVD< MatrixType, Options >::compute_impl	(	const MatrixType &	matrix,
		unsigned int	computationOptions
	)

private

                                                                                                         {
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
   std::cout << "\n\n\n================================================================================================="
                "=====================\n\n\n";
 #endif
   using std::abs;
  
   allocate(matrix.rows(), matrix.cols(), computationOptions);
  
   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
  
   //**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return
   if (matrix.cols() < m_algoswap) {
     smallSvd.compute(matrix);
     m_isInitialized = true;
     m_info = smallSvd.info();
     if (m_info == Success || m_info == NoConvergence) {
       if (computeU()) m_matrixU = smallSvd.matrixU();
       if (computeV()) m_matrixV = smallSvd.matrixV();
       m_singularValues = smallSvd.singularValues();
       m_nonzeroSingularValues = smallSvd.nonzeroSingularValues();
     }
     return *this;
   }
  
   //**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows
   RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
   if (!(numext::isfinite)(scale)) {
     m_isInitialized = true;
     m_info = InvalidInput;
     return *this;
   }
  
   if (numext::is_exactly_zero(scale)) scale = Literal(1);
  
   if (m_isTranspose)
     copyWorkspace = matrix.adjoint() / scale;
   else
     copyWorkspace = matrix / scale;
  
   //**** step 1 - Bidiagonalization.
   // If the problem is sufficiently rectangular, we perform R-Bidiagonalization: compute A = Q(R/0)
   // and then bidiagonalize R. Otherwise, if the problem is relatively square, we
   // bidiagonalize the input matrix directly.
   if (m_useQrDecomp) {
     qrDecomp.compute(copyWorkspace);
     reducedTriangle = qrDecomp.matrixQR().topRows(diagSize());
     reducedTriangle.template triangularView<StrictlyLower>().setZero();
     bid.compute(reducedTriangle);
   } else {
     bid.compute(copyWorkspace);
   }
  
   //**** step 2 - Divide & Conquer
   m_naiveU.setZero();
   m_naiveV.setZero();
   // FIXME this line involves a temporary matrix
   m_computed.topRows(diagSize()) = bid.bidiagonal().toDenseMatrix().transpose();
   m_computed.template bottomRows<1>().setZero();
   divide(0, diagSize() - 1, 0, 0, 0);
   if (m_info != Success && m_info != NoConvergence) {
     m_isInitialized = true;
     return *this;
   }
  
   //**** step 3 - Copy singular values and vectors
   for (int i = 0; i < diagSize(); i++) {
     RealScalar a = abs(m_computed.coeff(i, i));
     m_singularValues.coeffRef(i) = a * scale;
     if (a < considerZero) {
       m_nonzeroSingularValues = i;
       m_singularValues.tail(diagSize() - i - 1).setZero();
       break;
     } else if (i == diagSize() - 1) {
       m_nonzeroSingularValues = i + 1;
       break;
     }
   }
  
   //**** step 4 - Finalize unitaries U and V
   if (m_isTranspose)
     copyUV(bid.householderV(), bid.householderU(), m_naiveV, m_naiveU);
   else
     copyUV(bid.householderU(), bid.householderV(), m_naiveU, m_naiveV);
  
   if (m_useQrDecomp) {
     if (m_isTranspose && computeV())
       m_matrixV.applyOnTheLeft(qrDecomp.householderQ());
     else if (!m_isTranspose && computeU())
       m_matrixU.applyOnTheLeft(qrDecomp.householderQ());
   }
  
   m_isInitialized = true;
   return *this;
 }  // end compute

References a, abs(), i, Eigen::InvalidInput, Eigen::numext::is_exactly_zero(), Eigen::numext::isfinite(), matrix(), min, Eigen::NoConvergence, and Eigen::Success.

Referenced by Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD(), and Eigen::BDCSVD< MatrixType_, Options_ >::compute().

◆ computeBaseCase()

template<typename MatrixType , int Options>

template<typename SVDType >

void Eigen::BDCSVD< MatrixType, Options >::computeBaseCase	(	SVDType &	svd,
		Index	n,
		Index	firstCol,
		Index	firstRowW,
		Index	firstColW,
		Index	shift
	)

private

                                                                                 {
   svd.compute(m_computed.block(firstCol, firstCol, n + 1, n));
   m_info = svd.info();
   if (m_info != Success && m_info != NoConvergence) return;
   if (m_compU)
     m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = svd.matrixU();
   else {
     m_naiveU.row(0).segment(firstCol, n + 1).real() = svd.matrixU().row(0);
     m_naiveU.row(1).segment(firstCol, n + 1).real() = svd.matrixU().row(n);
   }
   if (m_compV) m_naiveV.block(firstRowW, firstColW, n, n).real() = svd.matrixV();
   m_computed.block(firstCol + shift, firstCol + shift, n + 1, n).setZero();
   m_computed.diagonal().segment(firstCol + shift, n) = svd.singularValues().head(n);
 }

References n, Eigen::NoConvergence, Eigen::Success, and svd().

◆ computeSingVals()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::computeSingVals	(	const ArrayRef &	col0,
		const ArrayRef &	diag,
		const IndicesRef &	perm,
		VectorType &	singVals,
		ArrayRef	shifts,
		ArrayRef	mus
	)

private

                                                                                                        {
   using std::abs;
   using std::sqrt;
   using std::swap;
  
   Index n = col0.size();
   Index actual_n = n;
   // Note that here actual_n is computed based on col0(i)==0 instead of diag(i)==0 as above
   // because 1) we have diag(i)==0 => col0(i)==0 and 2) if col0(i)==0, then diag(i) is already a singular value.
   while (actual_n > 1 && numext::is_exactly_zero(col0(actual_n - 1))) --actual_n;
  
   for (Index k = 0; k < n; ++k) {
     if (numext::is_exactly_zero(col0(k)) || 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) = k == 0 ? col0(0) : diag(k);
       mus(k) = Literal(0);
       shifts(k) = k == 0 ? col0(0) : diag(k);
       continue;
     }
  
     // otherwise, use secular equation to find singular value
     RealScalar left = diag(k);
     RealScalar right;  // was: = (k != actual_n-1) ? diag(k+1) : (diag(actual_n-1) + col0.matrix().norm());
     if (k == actual_n - 1)
       right = (diag(actual_n - 1) + col0.matrix().norm());
     else {
       // Skip deflated singular values,
       // recall that at this stage we assume that z[j]!=0 and all entries for which z[j]==0 have been put aside.
       // This should be equivalent to using perm[]
       Index l = k + 1;
       while (numext::is_exactly_zero(col0(l))) {
         ++l;
         eigen_internal_assert(l < actual_n);
       }
       right = diag(l);
     }
  
     // first decide whether it's closer to the left end or the right end
     RealScalar mid = left + (right - left) / Literal(2);
     RealScalar fMid = secularEq(mid, col0, diag, perm, diag, Literal(0));
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
     std::cout << "right-left = " << right - left << "\n";
     //     std::cout << "fMid = " << fMid << " " << secularEq(mid-left, col0, diag, perm, ArrayXr(diag-left), left)
     //                            << " " << secularEq(mid-right, col0, diag, perm, ArrayXr(diag-right), right)   <<
     //                            "\n";
     std::cout << "     = " << secularEq(left + RealScalar(0.000001) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.1) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.2) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.3) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.4) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.49) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.5) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.51) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.6) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.7) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.8) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.9) * (right - left), col0, diag, perm, diag, 0) << " "
               << secularEq(left + RealScalar(0.999999) * (right - left), col0, diag, perm, diag, 0) << "\n";
 #endif
     RealScalar shift = (k == actual_n - 1 || fMid > Literal(0)) ? left : right;
  
     // measure everything relative to shift
     Map<ArrayXr> diagShifted(m_workspace.data() + 4 * n, n);
     diagShifted = diag - shift;
  
     if (k != actual_n - 1) {
       // check that after the shift, f(mid) is still negative:
       RealScalar midShifted = (right - left) / RealScalar(2);
       // we can test exact equality here, because shift comes from `... ? left : right`
       if (numext::equal_strict(shift, right)) midShifted = -midShifted;
       RealScalar fMidShifted = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
       if (fMidShifted > 0) {
         // fMid was erroneous, fix it:
         shift = fMidShifted > Literal(0) ? left : right;
         diagShifted = diag - shift;
       }
     }
  
     // initial guess
     RealScalar muPrev, muCur;
     // we can test exact equality here, because shift comes from `... ? left : right`
     if (numext::equal_strict(shift, left)) {
       muPrev = (right - left) * RealScalar(0.1);
       if (k == actual_n - 1)
         muCur = right - left;
       else
         muCur = (right - left) * RealScalar(0.5);
     } else {
       muPrev = -(right - left) * RealScalar(0.1);
       muCur = -(right - left) * RealScalar(0.5);
     }
  
     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);
       swap(muPrev, muCur);
     }
  
     // 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 > Literal(0);
     while (!numext::is_exactly_zero(fCur) &&
            abs(muCur - muPrev) >
                Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(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) / (Literal(1) / muCur - Literal(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);
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
       eigen_internal_assert((numext::isfinite)(fZero));
 #endif
  
       muPrev = muCur;
       fPrev = fCur;
       muCur = muZero;
       fCur = fZero;
  
       // we can test exact equality here, because shift comes from `... ? left : right`
       if (numext::equal_strict(shift, left) && (muCur < Literal(0) || muCur > right - left)) useBisection = true;
       if (numext::equal_strict(shift, right) && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true;
       if (abs(fCur) > abs(fPrev)) useBisection = true;
     }
  
     // fall back on bisection method if rational interpolation did not work
     if (useBisection) {
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
       std::cout << "useBisection for k = " << k << ", actual_n = " << actual_n << "\n";
 #endif
       RealScalar leftShifted, rightShifted;
       // we can test exact equality here, because shift comes from `... ? left : right`
       if (numext::equal_strict(shift, left)) {
         // to avoid overflow, we must have mu > max(real_min, |z(k)|/sqrt(real_max)),
         // the factor 2 is to be more conservative
         leftShifted =
             numext::maxi<RealScalar>((std::numeric_limits<RealScalar>::min)(),
                                      Literal(2) * abs(col0(k)) / sqrt((std::numeric_limits<RealScalar>::max)()));
  
         // check that we did it right:
         eigen_internal_assert(
             (numext::isfinite)((col0(k) / leftShifted) * (col0(k) / (diag(k) + shift + leftShifted))));
         // I don't understand why the case k==0 would be special there:
         // if (k == 0) rightShifted = right - left; else
         rightShifted = (k == actual_n - 1)
                            ? right
                            : ((right - left) * RealScalar(0.51));  // theoretically we can take 0.5, but let's be safe
       } else {
         leftShifted = -(right - left) * RealScalar(0.51);
         if (k + 1 < n)
           rightShifted = -numext::maxi<RealScalar>((std::numeric_limits<RealScalar>::min)(),
                                                    abs(col0(k + 1)) / sqrt((std::numeric_limits<RealScalar>::max)()));
         else
           rightShifted = -(std::numeric_limits<RealScalar>::min)();
       }
  
       RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
       eigen_internal_assert(fLeft < Literal(0));
  
 #if defined EIGEN_BDCSVD_DEBUG_VERBOSE || defined EIGEN_BDCSVD_SANITY_CHECKS || defined EIGEN_INTERNAL_DEBUGGING
       RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift);
 #endif
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
       if (!(numext::isfinite)(fLeft))
         std::cout << "f(" << leftShifted << ") =" << fLeft << " ; " << left << " " << shift << " " << right << "\n";
       eigen_internal_assert((numext::isfinite)(fLeft));
  
       if (!(numext::isfinite)(fRight))
         std::cout << "f(" << rightShifted << ") =" << fRight << " ; " << left << " " << shift << " " << right << "\n";
         // eigen_internal_assert((numext::isfinite)(fRight));
 #endif
  
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
       if (!(fLeft * fRight < 0)) {
         std::cout << "f(leftShifted) using  leftShifted=" << leftShifted
                   << " ;  diagShifted(1:10):" << diagShifted.head(10).transpose() << "\n ; "
                   << "left==shift=" << bool(left == shift) << " ; left-shift = " << (left - shift) << "\n";
         std::cout << "k=" << k << ", " << fLeft << " * " << fRight << " == " << fLeft * fRight << "  ;  "
                   << "[" << left << " .. " << right << "] -> [" << leftShifted << " " << rightShifted
                   << "], shift=" << shift << " ,  f(right)=" << secularEq(0, col0, diag, perm, diagShifted, shift)
                   << " == " << secularEq(right, col0, diag, perm, diag, 0) << " == " << fRight << "\n";
       }
 #endif
       eigen_internal_assert(fLeft * fRight < Literal(0));
  
       if (fLeft < Literal(0)) {
         while (rightShifted - leftShifted > Literal(2) * NumTraits<RealScalar>::epsilon() *
                                                 numext::maxi<RealScalar>(abs(leftShifted), abs(rightShifted))) {
           RealScalar midShifted = (leftShifted + rightShifted) / Literal(2);
           fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
           eigen_internal_assert((numext::isfinite)(fMid));
  
           if (fLeft * fMid < Literal(0)) {
             rightShifted = midShifted;
           } else {
             leftShifted = midShifted;
             fLeft = fMid;
           }
         }
         muCur = (leftShifted + rightShifted) / Literal(2);
       } else {
         // We have a problem as shifting on the left or right give either a positive or negative value
         // at the middle of [left,right]...
         // Instead of abbording or entering an infinite loop,
         // let's just use the middle as the estimated zero-crossing:
         muCur = (right - left) * RealScalar(0.5);
         // we can test exact equality here, because shift comes from `... ? left : right`
         if (numext::equal_strict(shift, right)) muCur = -muCur;
       }
     }
  
     singVals[k] = shift + muCur;
     shifts[k] = shift;
     mus[k] = muCur;
  
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
     if (k + 1 < n)
       std::cout << "found " << singVals[k] << " == " << shift << " + " << muCur << " from " << diag(k) << " .. "
                 << diag(k + 1) << "\n";
 #endif
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
     eigen_internal_assert(k == 0 || singVals[k] >= singVals[k - 1]);
     eigen_internal_assert(singVals[k] >= diag(k));
 #endif
  
     // perturb singular value slightly if it equals diagonal entry to avoid division by zero later
     // (deflation is supposed to avoid this from happening)
     // - this does no seem to be necessary anymore -
     // if (singVals[k] == left) singVals[k] *= 1 + NumTraits<RealScalar>::epsilon();
     // if (singVals[k] == right) singVals[k] *= 1 - NumTraits<RealScalar>::epsilon();
   }
 }

References a, abs(), actual_n, b, diag, eigen_internal_assert, Eigen::numext::equal_strict(), Eigen::numext::is_exactly_zero(), Eigen::numext::isfinite(), k, max, min, plotDoE::mus, n, sqrt(), Eigen::swap(), and swap().

◆ computeSingVecs()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::computeSingVecs	(	const ArrayRef &	zhat,
		const ArrayRef &	diag,
		const IndicesRef &	perm,
		const VectorType &	singVals,
		const ArrayRef &	shifts,
		const ArrayRef &	mus,
		MatrixXr &	U,
		MatrixXr &	V
	)

private

                                                                                                  {
   Index n = zhat.size();
   Index m = perm.size();
  
   for (Index k = 0; k < n; ++k) {
     if (numext::is_exactly_zero(zhat(k))) {
       U.col(k) = VectorType::Unit(n + 1, k);
       if (m_compV) V.col(k) = VectorType::Unit(n, k);
     } else {
       U.col(k).setZero();
       for (Index l = 0; l < m; ++l) {
         Index i = perm(l);
         U(i, k) = zhat(i) / (((diag(i) - shifts(k)) - mus(k))) / ((diag(i) + singVals[k]));
       }
       U(n, k) = Literal(0);
       U.col(k).normalize();
  
       if (m_compV) {
         V.col(k).setZero();
         for (Index l = 1; l < m; ++l) {
           Index i = perm(l);
           V(i, k) = diag(i) * zhat(i) / (((diag(i) - shifts(k)) - mus(k))) / ((diag(i) + singVals[k]));
         }
         V(0, k) = Literal(-1);
         V.col(k).normalize();
       }
     }
   }
   U.col(n) = VectorType::Unit(n + 1, n);
 }

References diag, i, Eigen::numext::is_exactly_zero(), k, m, plotDoE::mus, n, RachelsAdvectionDiffusion::U, and V.

◆ computeSVDofM()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::computeSVDofM	(	Index	firstCol,
		Index	n,
		MatrixXr &	U,
		VectorType &	singVals,
		MatrixXr &	V
	)

private

                                                              {
   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
   using std::abs;
   ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n);
   m_workspace.head(n) = m_computed.block(firstCol, firstCol, n, n).diagonal();
   ArrayRef diag = m_workspace.head(n);
   diag(0) = Literal(0);
  
   // Allocate space for singular values and vectors
   singVals.resize(n);
   U.resize(n + 1, n + 1);
   if (m_compV) V.resize(n, n);
  
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
   if (col0.hasNaN() || diag.hasNaN()) std::cout << "\n\nHAS NAN\n\n";
 #endif
  
   // 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 && numext::is_exactly_zero(diag(actual_n - 1))) {
     --actual_n;
     eigen_internal_assert(numext::is_exactly_zero(col0(actual_n)));
   }
   Index m = 0;  // size of the deflated problem
   for (Index k = 0; k < actual_n; ++k)
     if (abs(col0(k)) > considerZero) m_workspaceI(m++) = k;
   Map<ArrayXi> perm(m_workspaceI.data(), m);
  
   Map<ArrayXr> shifts(m_workspace.data() + 1 * n, n);
   Map<ArrayXr> mus(m_workspace.data() + 2 * n, n);
   Map<ArrayXr> zhat(m_workspace.data() + 3 * n, n);
  
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
   std::cout << "computeSVDofM using:\n";
   std::cout << "  z: " << col0.transpose() << "\n";
   std::cout << "  d: " << diag.transpose() << "\n";
 #endif
  
   // Compute singVals, shifts, and 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";
  
   {
     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";
     eigen_internal_assert((((singVals.array() - (shifts + mus)) / singVals.array()).head(actual_n) >= 0).all());
     std::cout << "    check2 (>0)      : " << ((singVals.array() - diag) / singVals.array()).head(actual_n).transpose()
               << "\n\n";
     eigen_internal_assert((((singVals.array() - diag) / singVals.array()).head(actual_n) >= 0).all());
   }
 #endif
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(singVals.allFinite());
   eigen_internal_assert(mus.allFinite());
   eigen_internal_assert(shifts.allFinite());
 #endif
  
   // Compute zhat
   perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat);
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
   std::cout << "  zhat: " << zhat.transpose() << "\n";
 #endif
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(zhat.allFinite());
 #endif
  
   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";
   std::cout << "V^T V: " << (V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(), V.cols()))).norm() << "\n";
 #endif
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(m_naiveU.allFinite());
   eigen_internal_assert(m_naiveV.allFinite());
   eigen_internal_assert(m_computed.allFinite());
   eigen_internal_assert(U.allFinite());
   eigen_internal_assert(V.allFinite());
 //   eigen_internal_assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() <
 //   100*NumTraits<RealScalar>::epsilon() * n); eigen_internal_assert((V.transpose() * V -
 //   MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 100*NumTraits<RealScalar>::epsilon() * n);
 #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));
     }
   }
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   {
     bool singular_values_sorted =
         (((singVals.segment(1, actual_n - 1) - singVals.head(actual_n - 1))).array() >= 0).all();
     if (!singular_values_sorted)
       std::cout << "Singular values are not sorted: " << singVals.segment(1, actual_n).transpose() << "\n";
     eigen_internal_assert(singular_values_sorted);
   }
 #endif
  
   // Reverse order so that singular values in increased order
   // Because of deflation, the zeros singular-values are already at the end
   singVals.head(actual_n).reverseInPlace();
   U.leftCols(actual_n).rowwise().reverseInPlace();
   if (m_compV) V.leftCols(actual_n).rowwise().reverseInPlace();
  
 #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
 }

References abs(), actual_n, Eigen::placeholders::all, diag, eigen_internal_assert, i, Eigen::numext::is_exactly_zero(), k, m, min, plotDoE::mus, n, Eigen::PlainObjectBase< Derived >::resize(), Eigen::SVDBase< Derived >::singularValues(), Eigen::swap(), swap(), RachelsAdvectionDiffusion::U, and V.

◆ copyUV()

template<typename MatrixType , int Options>

template<typename HouseholderU , typename HouseholderV , typename NaiveU , typename NaiveV >

void Eigen::BDCSVD< MatrixType, Options >::copyUV	(	const HouseholderU &	householderU,
		const HouseholderV &	householderV,
		const NaiveU &	naiveU,
		const NaiveV &	naivev
	)

private

                                                                                      {
   // Note exchange of U and V: m_matrixU is set from m_naiveV and vice versa
   if (computeU()) {
     Index Ucols = m_computeThinU ? diagSize() : rows();
     m_matrixU = MatrixX::Identity(rows(), Ucols);
     m_matrixU.topLeftCorner(diagSize(), diagSize()) =
         naiveV.template cast<Scalar>().topLeftCorner(diagSize(), diagSize());
     // FIXME the following conditionals involve temporary buffers
     if (m_useQrDecomp)
       m_matrixU.topLeftCorner(householderU.cols(), diagSize()).applyOnTheLeft(householderU);
     else
       m_matrixU.applyOnTheLeft(householderU);
   }
   if (computeV()) {
     Index Vcols = m_computeThinV ? diagSize() : cols();
     m_matrixV = MatrixX::Identity(cols(), Vcols);
     m_matrixV.topLeftCorner(diagSize(), diagSize()) =
         naiveU.template cast<Scalar>().topLeftCorner(diagSize(), diagSize());
     // FIXME the following conditionals involve temporary buffers
     if (m_useQrDecomp)
       m_matrixV.topLeftCorner(householderV.cols(), diagSize()).applyOnTheLeft(householderV);
     else
       m_matrixV.applyOnTheLeft(householderV);
   }
 }

References cols, and rows.

◆ deflation()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::deflation	(	Index	firstCol,
		Index	lastCol,
		Index	k,
		Index	firstRowW,
		Index	firstColW,
		Index	shift
	)

private

                                                          {
   using std::abs;
   using std::sqrt;
   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);
  
   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
   RealScalar maxDiag = diag.tail((std::max)(Index(1), length - 1)).cwiseAbs().maxCoeff();
   RealScalar epsilon_strict = numext::maxi<RealScalar>(considerZero, NumTraits<RealScalar>::epsilon() * maxDiag);
   RealScalar epsilon_coarse =
       Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag);
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(m_naiveU.allFinite());
   eigen_internal_assert(m_naiveV.allFinite());
   eigen_internal_assert(m_computed.allFinite());
 #endif
  
 #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_coarse) {
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
     std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon_coarse << "\n";
 #endif
     diag(0) = epsilon_coarse;
   }
  
   // condition 4.2
   for (Index i = 1; i < length; ++i)
     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_strict
                 << "  (diag(" << i << ")=" << diag(i) << ")\n";
 #endif
       col0(i) = Literal(0);
     }
  
   // condition 4.3
   for (Index i = 1; i < length; i++)
     if (diag(i) < epsilon_coarse) {
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
       std::cout << "deflation 4.3, cancel z(" << i << ")=" << col0(i) << " because diag(" << i << ")=" << diag(i)
                 << " < " << epsilon_coarse << "\n";
 #endif
       deflation43(firstCol, shift, i, length);
     }
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(m_naiveU.allFinite());
   eigen_internal_assert(m_naiveV.allFinite());
   eigen_internal_assert(m_computed.allFinite());
 #endif
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
   std::cout << "to be sorted: " << diag.transpose() << "\n\n";
   std::cout << "            : " << col0.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.
     const bool total_deflation = (col0.tail(length - 1).array().abs() < considerZero).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 = m_workspaceI.data();
     {
       permutation[0] = 0;
       Index p = 1;
  
       // Move deflated diagonal entries at the end.
       for (Index i = 1; i < length; ++i)
         if (diag(i) < considerZero) permutation[p++] = i;
  
       Index i = 1, j = k + 1;
       for (; p < length; ++p) {
         if (i > k)
           permutation[p] = j++;
         else if (j >= length)
           permutation[p] = i++;
         else if (diag(i) < diag(j))
           permutation[p] = j++;
         else
           permutation[p] = i++;
       }
     }
  
     // 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) < considerZero || 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 = m_workspaceI.data() + length;
     Index* realCol = m_workspaceI.data() + 2 * length;
  
     for (int pos = 0; pos < length; pos++) {
       realCol[pos] = pos;
       realInd[pos] = pos;
     }
  
     for (Index i = total_deflation ? 0 : 1; i < length; i++) {
       const Index pi = permutation[length - (total_deflation ? i + 1 : i)];
       const Index J = realCol[pi];
  
       using std::swap;
       // swap diagonal and first column entries:
       swap(diag(i), diag(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));
       else
         m_naiveU.col(firstCol + i).segment(0, 2).swap(m_naiveU.col(firstCol + J).segment(0, 2));
       if (m_compV)
         m_naiveV.col(firstColW + i)
             .segment(firstRowW, length)
             .swap(m_naiveV.col(firstColW + J).segment(firstRowW, length));
  
       // update real pos
       const Index realI = realInd[i];
       realCol[realI] = J;
       realCol[pi] = i;
       realInd[J] = realI;
       realInd[i] = pi;
     }
   }
 #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;
     // Find last non-deflated entry.
     while (i > 0 && (diag(i) < considerZero || abs(col0(i)) < considerZero)) --i;
  
     for (; i > 1; --i)
       if ((diag(i) - diag(i - 1)) < epsilon_strict) {
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
         std::cout << "deflation 4.4 with i = " << i << " because " << diag(i) << " - " << diag(i - 1)
                   << " == " << (diag(i) - diag(i - 1)) << " < " << epsilon_strict << "\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, i - 1, length);
       }
   }
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   for (Index j = 2; j < length; ++j) eigen_internal_assert(diag(j - 1) <= diag(j) || abs(diag(j)) < considerZero);
 #endif
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(m_naiveU.allFinite());
   eigen_internal_assert(m_naiveV.allFinite());
   eigen_internal_assert(m_computed.allFinite());
 #endif
 }  // end deflation

References abs(), Eigen::placeholders::all, diag, Eigen::Dynamic, eigen_internal_assert, oomph::SarahBL::epsilon, i, J, j, k, max, min, p, constants::pi, sqrt(), Eigen::swap(), and swap().

◆ deflation43()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::deflation43	(	Index	firstCol,
		Index	shift,
		Index	i,
		Index	size
	)

private

                                                                                               {
   using std::abs;
   using std::pow;
   using std::sqrt;
   Index start = firstCol + shift;
   RealScalar c = m_computed(start, start);
   RealScalar s = m_computed(start + i, start);
   RealScalar r = numext::hypot(c, s);
   if (numext::is_exactly_zero(r)) {
     m_computed(start + i, start + i) = Literal(0);
     return;
   }
   m_computed(start, start) = r;
   m_computed(start + i, start) = Literal(0);
   m_computed(start + i, start + i) = Literal(0);
  
   JacobiRotation<RealScalar> J(c / r, -s / r);
   if (m_compU)
     m_naiveU.middleRows(firstCol, size + 1).applyOnTheRight(firstCol, firstCol + i, J);
   else
     m_naiveU.applyOnTheRight(firstCol, firstCol + i, J);
 }  // end deflation 43

References abs(), calibrate::c, i, Eigen::numext::is_exactly_zero(), J, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, s, size, sqrt(), and oomph::CumulativeTimings::start().

◆ deflation44()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::deflation44	(	Index	firstColu,
		Index	firstColm,
		Index	firstRowW,
		Index	firstColW,
		Index	i,
		Index	j,
		Index	size
	)

private

                                                                             {
   using std::abs;
   using std::conj;
   using std::pow;
   using std::sqrt;
  
   RealScalar s = m_computed(firstColm + i, firstColm);
   RealScalar c = m_computed(firstColm + j, firstColm);
   RealScalar r = numext::hypot(c, s);
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
   std::cout << "deflation 4.4: " << i << "," << j << " -> " << c << " " << s << " " << r << " ; "
             << m_computed(firstColm + i - 1, firstColm) << " " << m_computed(firstColm + i, firstColm) << " "
             << m_computed(firstColm + i + 1, firstColm) << " " << m_computed(firstColm + i + 2, firstColm) << "\n";
   std::cout << m_computed(firstColm + i - 1, firstColm + i - 1) << " " << m_computed(firstColm + i, firstColm + i)
             << " " << m_computed(firstColm + i + 1, firstColm + i + 1) << " "
             << m_computed(firstColm + i + 2, firstColm + i + 2) << "\n";
 #endif
   if (numext::is_exactly_zero(r)) {
     m_computed(firstColm + j, firstColm + j) = m_computed(firstColm + i, firstColm + i);
     return;
   }
   c /= r;
   s /= r;
   m_computed(firstColm + j, firstColm) = r;
   m_computed(firstColm + j, firstColm + j) = m_computed(firstColm + i, firstColm + i);
   m_computed(firstColm + i, firstColm) = Literal(0);
  
   JacobiRotation<RealScalar> J(c, -s);
   if (m_compU)
     m_naiveU.middleRows(firstColu, size + 1).applyOnTheRight(firstColu + j, firstColu + i, J);
   else
     m_naiveU.applyOnTheRight(firstColu + j, firstColu + i, J);
   if (m_compV) m_naiveV.middleRows(firstRowW, size).applyOnTheRight(firstColW + j, firstColW + i, J);
 }  // end deflation 44

References abs(), calibrate::c, conj(), i, Eigen::numext::is_exactly_zero(), J, j, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, s, size, and sqrt().

◆ divide()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::divide	(	Index	firstCol,
		Index	lastCol,
		Index	firstRowW,
		Index	firstColW,
		Index	shift
	)

private

                                                                                                                      {
   // requires rows = cols + 1;
   using std::abs;
   using std::pow;
   using std::sqrt;
   const Index n = lastCol - firstCol + 1;
   const Index k = n / 2;
   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
   RealScalar alphaK;
   RealScalar betaK;
   RealScalar r0;
   RealScalar lambda, phi, c0, s0;
   VectorType l, f;
   // We use the other algorithm which is more efficient for small
   // matrices.
   if (n < m_algoswap) {
     // FIXME this block involves temporaries
     if (m_compV) {
       JacobiSVD<MatrixXr, ComputeFullU | ComputeFullV> baseSvd;
       computeBaseCase(baseSvd, n, firstCol, firstRowW, firstColW, shift);
     } else {
       JacobiSVD<MatrixXr, ComputeFullU> baseSvd;
       computeBaseCase(baseSvd, n, firstCol, firstRowW, firstColW, shift);
     }
     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);
   if (m_info != Success && m_info != NoConvergence) return;
   divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1);
   if (m_info != Success && m_info != NoConvergence) return;
  
   if (m_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 (m_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 (m_compV) m_naiveV(firstRowW + k, firstColW) = Literal(1);
   if (r0 < considerZero) {
     c0 = Literal(1);
     s0 = Literal(0);
   } else {
     c0 = alphaK * lambda / r0;
     s0 = betaK * phi / r0;
   }
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(m_naiveU.allFinite());
   eigen_internal_assert(m_naiveV.allFinite());
   eigen_internal_assert(m_computed.allFinite());
 #endif
  
   if (m_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();
   }
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(m_naiveU.allFinite());
   eigen_internal_assert(m_naiveV.allFinite());
   eigen_internal_assert(m_computed.allFinite());
 #endif
  
   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();
  
 #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";
   eigen_internal_assert((tmp1 - tmp2).matrix().norm() < 1e-14 * tmp2.matrix().norm());
 //   eigen_internal_assert(count<681);
 //   eigen_internal_assert(((tmp1-tmp2).abs()<1e-13*tmp2.abs()).all());
 #endif
  
   // Third part: compute SVD of combined matrix
   MatrixXr UofSVD, VofSVD;
   VectorType singVals;
   computeSVDofM(firstCol + shift, n, UofSVD, singVals, VofSVD);
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(UofSVD.allFinite());
   eigen_internal_assert(VofSVD.allFinite());
 #endif
  
   if (m_compU)
     structured_update(m_naiveU.block(firstCol, firstCol, n + 1, n + 1), UofSVD, (n + 2) / 2);
   else {
     Map<Matrix<RealScalar, 2, Dynamic>, Aligned> tmp(m_workspace.data(), 2, n + 1);
     tmp.noalias() = m_naiveU.middleCols(firstCol, n + 1) * UofSVD;
     m_naiveU.middleCols(firstCol, n + 1) = tmp;
   }
  
   if (m_compV) structured_update(m_naiveV.block(firstRowW, firstColW, n, n), VofSVD, (n + 1) / 2);
  
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   eigen_internal_assert(m_naiveU.allFinite());
   eigen_internal_assert(m_naiveV.allFinite());
   eigen_internal_assert(m_computed.allFinite());
 #endif
  
   m_computed.block(firstCol + shift, firstCol + shift, n, n).setZero();
   m_computed.block(firstCol + shift, firstCol + shift, n, n).diagonal() = singVals;
 }  // end divide

References abs(), Eigen::Aligned, e(), eigen_internal_assert, f(), i, k, lambda, matrix(), min, n, Eigen::NoConvergence, Eigen::bfloat16_impl::pow(), sqrt(), Eigen::Success, tmp, and anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose().

◆ perturbCol0()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::perturbCol0	(	const ArrayRef &	col0,
		const ArrayRef &	diag,
		const IndicesRef &	perm,
		const VectorType &	singVals,
		const ArrayRef &	shifts,
		const ArrayRef &	mus,
		ArrayRef	zhat
	)

private

                                                              {
   using std::sqrt;
   Index n = col0.size();
   Index m = perm.size();
   if (m == 0) {
     zhat.setZero();
     return;
   }
   Index lastIdx = perm(m - 1);
   // The offset permits to skip deflated entries while computing zhat
   for (Index k = 0; k < n; ++k) {
     if (numext::is_exactly_zero(col0(k)))  // deflated
       zhat(k) = Literal(0);
     else {
       // see equation (3.6)
       RealScalar dk = diag(k);
       RealScalar prod = (singVals(lastIdx) + dk) * (mus(lastIdx) + (shifts(lastIdx) - dk));
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
       if (prod < 0) {
         std::cout << "k = " << k << " ;  z(k)=" << col0(k) << ", diag(k)=" << dk << "\n";
         std::cout << "prod = "
                   << "(" << singVals(lastIdx) << " + " << dk << ") * (" << mus(lastIdx) << " + (" << shifts(lastIdx)
                   << " - " << dk << "))"
                   << "\n";
         std::cout << "     = " << singVals(lastIdx) + dk << " * " << mus(lastIdx) + (shifts(lastIdx) - dk) << "\n";
       }
       eigen_internal_assert(prod >= 0);
 #endif
  
       for (Index l = 0; l < m; ++l) {
         Index i = perm(l);
         if (i != k) {
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
           if (i >= k && (l == 0 || l - 1 >= m)) {
             std::cout << "Error in perturbCol0\n";
             std::cout << "  " << k << "/" << n << " " << l << "/" << m << " " << i << "/" << n << " ; " << col0(k)
                       << " " << diag(k) << " "
                       << "\n";
             std::cout << "  " << diag(i) << "\n";
             Index j = (i < k /*|| l==0*/) ? i : perm(l - 1);
             std::cout << "  "
                       << "j=" << j << "\n";
           }
 #endif
           Index j = i < k ? i : l > 0 ? perm(l - 1) : i;
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
           if (!(dk != Literal(0) || diag(i) != Literal(0))) {
             std::cout << "k=" << k << ", i=" << i << ", l=" << l << ", perm.size()=" << perm.size() << "\n";
           }
           eigen_internal_assert(dk != Literal(0) || diag(i) != Literal(0));
 #endif
           prod *= ((singVals(j) + dk) / ((diag(i) + dk))) * ((mus(j) + (shifts(j) - dk)) / ((diag(i) - dk)));
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
           eigen_internal_assert(prod >= 0);
 #endif
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
           if (i != k &&
               numext::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";
 #endif
         }
       }
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
       std::cout << "zhat(" << k << ") =  sqrt( " << prod << ")  ;  " << (singVals(lastIdx) + dk) << " * "
                 << mus(lastIdx) + shifts(lastIdx) << " - " << dk << "\n";
 #endif
       RealScalar tmp = sqrt(prod);
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
       eigen_internal_assert((numext::isfinite)(tmp));
 #endif
       zhat(k) = col0(k) > Literal(0) ? RealScalar(tmp) : RealScalar(-tmp);
     }
   }
 }

References Eigen::numext::abs(), diag, eigen_internal_assert, i, Eigen::numext::is_exactly_zero(), Eigen::numext::isfinite(), j, k, m, plotDoE::mus, n, Eigen::prod(), sqrt(), and tmp.

◆ rows()

template<typename MatrixType_ , int Options_>

EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index Eigen::EigenBase< Derived >::rows

inline

Returns: the number of rows.

See also: cols(), RowsAtCompileTime

59 { return derived().rows(); }

Referenced by gdb.printers._MatrixEntryIterator::__next__(), Eigen::BDCSVD< MatrixType_, Options_ >::BDCSVD(), gdb.printers.EigenMatrixPrinter::children(), gdb.printers.EigenSparseMatrixPrinter::children(), gdb.printers.EigenMatrixPrinter::to_string(), and gdb.printers.EigenSparseMatrixPrinter::to_string().

◆ secularEq()

template<typename MatrixType , int Options>

BDCSVD< MatrixType, Options >::RealScalar Eigen::BDCSVD< MatrixType, Options >::secularEq	(	RealScalar	x,
		const ArrayRef &	col0,
		const ArrayRef &	diag,
		const IndicesRef &	perm,
		const ArrayRef &	diagShifted,
		RealScalar	shift
	)

staticprivate

                       {
   Index m = perm.size();
   RealScalar res = Literal(1);
   for (Index i = 0; i < m; ++i) {
     Index j = perm(i);
     // The following expression could be rewritten to involve only a single division,
     // but this would make the expression more sensitive to overflow.
     res += (col0(j) / (diagShifted(j) - mu)) * (col0(j) / (diag(j) + shift + mu));
   }
   return res;
 }

References diag, i, j, m, Global_Parameters::mu, and res.

◆ setSwitchSize()

template<typename MatrixType_ , int Options_>

void Eigen::BDCSVD< MatrixType_, Options_ >::setSwitchSize ( int s )

inline

                             {
     eigen_assert(s >= 3 && "BDCSVD the size of the algo switch has to be at least 3.");
     m_algoswap = s;
   }

References eigen_assert, Eigen::BDCSVD< MatrixType_, Options_ >::m_algoswap, and s.

Referenced by compare_bdc_jacobi().

◆ structured_update()

template<typename MatrixType , int Options>

void Eigen::BDCSVD< MatrixType, Options >::structured_update	(	Block< MatrixXr, Dynamic, Dynamic >	A,
		const MatrixXr &	B,
		Index	n1
	)

private

Performs A = A * B exploiting the special structure of the matrix A. Splitting A as: A = [A1] [A2] such that A1.rows()==n1, then we assume that at least half of the columns of A1 and A2 are zeros. We can thus pack them prior to the the matrix product. However, this is only worth the effort if the matrix is large enough.

                                                                                                                     {
   Index n = A.rows();
   if (n > 100) {
     // If the matrices are large enough, let's exploit the sparse structure of A by
     // splitting it in half (wrt n1), and packing the non-zero columns.
     Index n2 = n - n1;
     Map<MatrixXr> A1(m_workspace.data(), n1, n);
     Map<MatrixXr> A2(m_workspace.data() + n1 * n, n2, n);
     Map<MatrixXr> B1(m_workspace.data() + n * n, n, n);
     Map<MatrixXr> B2(m_workspace.data() + 2 * n * n, n, n);
     Index k1 = 0, k2 = 0;
     for (Index j = 0; j < n; ++j) {
       if ((A.col(j).head(n1).array() != Literal(0)).any()) {
         A1.col(k1) = A.col(j).head(n1);
         B1.row(k1) = B.row(j);
         ++k1;
       }
       if ((A.col(j).tail(n2).array() != Literal(0)).any()) {
         A2.col(k2) = A.col(j).tail(n2);
         B2.row(k2) = B.row(j);
         ++k2;
       }
     }
  
     A.topRows(n1).noalias() = A1.leftCols(k1) * B1.topRows(k1);
     A.bottomRows(n2).noalias() = A2.leftCols(k2) * B2.topRows(k2);
   } else {
     Map<MatrixXr, Aligned> tmp(m_workspace.data(), n, n);
     tmp.noalias() = A * B;
     A = tmp;
   }
 }

References j, n, Eigen::PlainObjectBase< Derived >::rows(), and tmp.

Member Data Documentation

◆ bid

template<typename MatrixType_ , int Options_>

internal::UpperBidiagonalization<MatrixX> Eigen::BDCSVD< MatrixType_, Options_ >::bid

protected

◆ copyWorkspace

template<typename MatrixType_ , int Options_>

MatrixX Eigen::BDCSVD< MatrixType_, Options_ >::copyWorkspace

protected

◆ m_algoswap

template<typename MatrixType_ , int Options_>

int Eigen::BDCSVD< MatrixType_, Options_ >::m_algoswap

protected

Referenced by Eigen::BDCSVD< MatrixType_, Options_ >::setSwitchSize().

◆ m_compU

template<typename MatrixType_ , int Options_>

bool Eigen::BDCSVD< MatrixType_, Options_ >::m_compU

protected

◆ m_computed

template<typename MatrixType_ , int Options_>

MatrixXr Eigen::BDCSVD< MatrixType_, Options_ >::m_computed

protected

◆ m_compV

template<typename MatrixType_ , int Options_>

bool Eigen::BDCSVD< MatrixType_, Options_ >::m_compV

protected

◆ m_isTranspose

template<typename MatrixType_ , int Options_>

bool Eigen::BDCSVD< MatrixType_, Options_ >::m_isTranspose

protected

◆ m_naiveU

template<typename MatrixType_ , int Options_>

MatrixXr Eigen::BDCSVD< MatrixType_, Options_ >::m_naiveU

protected

◆ m_naiveV

template<typename MatrixType_ , int Options_>

MatrixXr Eigen::BDCSVD< MatrixType_, Options_ >::m_naiveV

protected

◆ m_nRec

template<typename MatrixType_ , int Options_>

Index Eigen::BDCSVD< MatrixType_, Options_ >::m_nRec

protected

◆ m_numIters

template<typename MatrixType_ , int Options_>

int Eigen::BDCSVD< MatrixType_, Options_ >::m_numIters

◆ m_useQrDecomp

template<typename MatrixType_ , int Options_>

bool Eigen::BDCSVD< MatrixType_, Options_ >::m_useQrDecomp

protected

◆ m_workspace

template<typename MatrixType_ , int Options_>

ArrayXr Eigen::BDCSVD< MatrixType_, Options_ >::m_workspace

protected

◆ m_workspaceI

template<typename MatrixType_ , int Options_>

ArrayXi Eigen::BDCSVD< MatrixType_, Options_ >::m_workspaceI

protected

◆ qrDecomp

template<typename MatrixType_ , int Options_>

HouseholderQR<MatrixX> Eigen::BDCSVD< MatrixType_, Options_ >::qrDecomp

protected

◆ reducedTriangle

template<typename MatrixType_ , int Options_>

MatrixX Eigen::BDCSVD< MatrixType_, Options_ >::reducedTriangle

protected

◆ smallSvd

template<typename MatrixType_ , int Options_>

JacobiSVD<MatrixType, ComputationOptions> Eigen::BDCSVD< MatrixType_, Options_ >::smallSvd

protected

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Private Types

Private Member Functions

Static Private Member Functions

Additional Inherited Members

Detailed Description

template<typename MatrixType_, int Options_> class Eigen::BDCSVD< MatrixType_, Options_ >

Member Typedef Documentation

◆ ArrayRef

◆ ArrayXi

◆ ArrayXr

◆ Base

◆ Index

◆ IndicesRef

◆ Literal

◆ MatrixType

◆ MatrixUType

◆ MatrixVType

◆ MatrixX

◆ MatrixXr

◆ RealScalar

◆ Scalar

◆ SingularValuesType

◆ VectorType

Member Enumeration Documentation

◆ anonymous enum

Constructor & Destructor Documentation

◆ BDCSVD() [1/5]

◆ BDCSVD() [2/5]

◆ BDCSVD() [3/5]

◆ BDCSVD() [4/5]

◆ BDCSVD() [5/5]

◆ ~BDCSVD()

Member Function Documentation

◆ allocate()

◆ cols()

◆ compute() [1/2]

◆ compute() [2/2]

◆ compute_impl()

◆ computeBaseCase()

◆ computeSingVals()

◆ computeSingVecs()

◆ computeSVDofM()

◆ copyUV()

◆ deflation()

◆ deflation43()

◆ deflation44()

◆ divide()

◆ perturbCol0()

◆ rows()

◆ secularEq()

◆ setSwitchSize()

◆ structured_update()

Member Data Documentation

◆ bid

◆ copyWorkspace

◆ m_algoswap

◆ m_compU

◆ m_computed

◆ m_compV

◆ m_isTranspose

◆ m_naiveU

◆ m_naiveV

◆ m_nRec

◆ m_numIters

◆ m_useQrDecomp

◆ m_workspace

◆ m_workspaceI

◆ qrDecomp

◆ reducedTriangle

◆ smallSvd

template<typename MatrixType_, int Options_>
class Eigen::BDCSVD< MatrixType_, Options_ >