SVD: lower-level API

Overview

The stages and algorithm-specific entry points behind gesvd (gebrd, bdsqr, bdsdc, ormbr, blocked, CTA and Jacobi SVD) and their workspace queries.

Functions

gebrd_unblocked()

template<Backend B, typename T>

Event batchlas::gebrd_unblocked	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		const VectorView< typename base_type< T >::type > &	d,
		const VectorView< typename base_type< T >::type > &	e,
		const VectorView< T > &	tauq,
		const VectorView< T > &	taup )

Unblocked GEBRD-like reduction to real bidiagonal form: reduces each square A in place, returning (d,e) and the Householder scalars (tauq,taup).

Computes \( A = Q B P^H \) with B upper bidiagonal: d gets its n diagonal and e its n-1 superdiagonal entries (real), and the reflectors defining Q and P are left in A with their scalars in tauq and taup, as in LAPACK ?gebrd.

See also

GESVD (where the unblocked form stops paying)

gebrd_cta()

template<Backend B, typename T>

Event batchlas::gebrd_cta	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		const VectorView< typename base_type< T >::type > &	d,
		const VectorView< typename base_type< T >::type > &	e,
		const VectorView< T > &	tauq,
		const VectorView< T > &	taup,
		size_t	cta_wg_size_multiplier = 1 )

CTA-parallel small-matrix GEBRD for real square matrices, 1 <= n <= 32, where one CTA cooperatively reduces one matrix to upper bidiagonal form.

gebrd_blocked()

template<Backend B, typename T>

Event batchlas::gebrd_blocked	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		const VectorView< typename base_type< T >::type > &	d,
		const VectorView< typename base_type< T >::type > &	e,
		const VectorView< T > &	tauq,
		const VectorView< T > &	taup,
		const Span< std::byte > &	ws,
		int32_t	block_size = 16 )

Blocked GEBRD for real square dense matrices: DLABRD-style panel + GEMM.

Output layout as for gebrd_unblocked().

gebrd_blocked_buffer_size()

template<Backend B, typename T>

size_t batchlas::gebrd_blocked_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		const VectorView< typename base_type< T >::type > &	d,
		const VectorView< typename base_type< T >::type > &	e,
		const VectorView< T > &	tauq,
		const VectorView< T > &	taup,
		int32_t	block_size = 16 )

Required workspace, in bytes, for gebrd_blocked(); arguments as for the call.

bdsqr() [1/2]

template<Backend B, typename T>

Event batchlas::bdsqr	(	Queue &	ctx,
		const VectorView< T > &	d,
		const VectorView< T > &	e,
		Span< T >	singular_values_out,
		const Span< std::byte > &	ws,
		bool	sort_desc = true,
		Span< int32_t >	info = Span< int32_t >() )

Bidiagonal QR iteration for a real upper bidiagonal matrix.

The matrix overload also accumulates the alternating Givens rotations into u and vh, matching LAPACK BDSQR: it returns u * Q and P^T * vh.

Parameters

d,e	diagonal (n) and superdiagonal (n-1) per batch item
singular_values_out	n singular values per item, packed
ws	at least bdsqr_buffer_size() bytes
sort_desc	sort the singular values (and vectors) in descending order
info	per-item convergence status, or empty; see docs/cpp-api.md

bdsqr() [2/2]

template<Backend B, typename T>

Event batchlas::bdsqr	(	Queue &	ctx,
		const VectorView< T > &	d,
		const VectorView< T > &	e,
		Span< T >	singular_values_out,
		const Span< std::byte > &	ws,
		const MatrixView< T, MatrixFormat::Dense > &	u,
		const MatrixView< T, MatrixFormat::Dense > &	vh,
		bool	sort_desc = true,
		Span< int32_t >	info = Span< int32_t >() )

bdsqr() that also accumulates the rotations: u <- u * Q, vh <- P^T * vh.

Other parameters as for the values-only form.

bdsqr_buffer_size() [1/2]

template<typename T>

size_t batchlas::bdsqr_buffer_size	(	Queue &	ctx,
		const VectorView< T > &	d,
		const VectorView< T > &	e,
		Span< T >	singular_values_out )

Required workspace, in bytes, for bdsqr().

bdsqr_buffer_size() [2/2]

template<typename T>

size_t batchlas::bdsqr_buffer_size ( Queue & ctx, const VectorView< T > & d, const VectorView< T > & e, Span< T > singular_values_out, const MatrixView< T, MatrixFormat::Dense > & u, const MatrixView< T, MatrixFormat::Dense > & vh )

inline

bdsqr_buffer_size() for the vector-accumulating form; u and vh do not change it.

Definition at line 1514 of file extensions.hh.

bdsdc() [1/2]

template<Backend B, typename T>

Event batchlas::bdsdc	(	Queue &	ctx,
		const VectorView< T > &	d,
		const VectorView< T > &	e,
		Span< T >	singular_values_out,
		const Span< std::byte > &	ws,
		bool	sort_desc = true,
		Span< int32_t >	info = Span< int32_t >() )

Bidiagonal divide-and-conquer SVD for a real upper bidiagonal matrix.

Same problem as bdsqr, via the symmetric tridiagonal eigenproblem of the interleaved Golub-Kahan form (order 2n) handed to stedc.

Unlike bdsqr, which accumulates into whatever it is handed (u <- u*Q), bdsdc WRITES the leading n x n block of u and of vh and leaves the rest untouched – seed them with the identity if the trailing columns matter. Workspace is dominated by a 2n x 2n eigenvector matrix per batch item. Parameters as for bdsqr().

See also

gesvd design, section "gesvd design: why bdsdc goes through Golub-Kahan"

bdsdc() [2/2]

template<Backend B, typename T>

Event batchlas::bdsdc	(	Queue &	ctx,
		const VectorView< T > &	d,
		const VectorView< T > &	e,
		Span< T >	singular_values_out,
		const Span< std::byte > &	ws,
		const MatrixView< T, MatrixFormat::Dense > &	u,
		const MatrixView< T, MatrixFormat::Dense > &	vh,
		bool	sort_desc = true,
		Span< int32_t >	info = Span< int32_t >() )

bdsdc() that also writes the leading n x n blocks of u and vh.

bdsdc_buffer_size()

template<Backend B, typename T>

size_t batchlas::bdsdc_buffer_size	(	Queue &	ctx,
		const VectorView< T > &	d,
		const VectorView< T > &	e,
		Span< T >	singular_values_out,
		bool	want_vectors )

Required workspace, in bytes, for bdsdc(); want_vectors selects the overload.

ormbr()

template<Backend B, typename T>

Event batchlas::ormbr	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		const VectorView< T > &	tau,
		const MatrixView< T, MatrixFormat::Dense > &	c,
		char	vect,
		Side	side,
		Transpose	trans,
		const Span< std::byte > &	ws,
		int32_t	block_size = 32 )

ORMBR/UNMBR-style application of bidiagonal reduction reflectors.

Supports vect='Q' (tauq, CTA/blocked ORMQR path) and vect='P' (taup, blocked compact-WY path over the right reflectors).

Overwrites c with op(Q) * c, c * op(Q), op(P) * c or c * op(P) according to vect, side and trans, where Q/P are the reflectors that gebrd_blocked() (or gebrd_unblocked()) left in a and tau, as LAPACK ?ormbr / ?unmbr: 'Q' is of order a.rows(), 'P' of order a.cols().

Precondition

tau is unit-stride and packed by batch

Exceptions

batchlas::invalid_argument	on mismatched batch sizes or orders, a vect other than 'Q'/'P', or a short or strided tau
batchlas::unsupported	for Transpose::Trans with complex T and 'P' (use ConjTrans)

ormbr_buffer_size()

template<Backend B, typename T>

size_t batchlas::ormbr_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		const VectorView< T > &	tau,
		const MatrixView< T, MatrixFormat::Dense > &	c,
		char	vect,
		Side	side,
		Transpose	trans,
		int32_t	block_size = 32 )

Required workspace, in bytes, for ormbr(); arguments as for the call.

gesvd_blocked() [1/2]

template<Backend B, typename T>

Event batchlas::gesvd_blocked	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a_in,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		const Span< std::byte > &	ws,
		Span< int32_t >	info = Span< int32_t >() )

Blocked native SVD for real dense matrices: GEBRD-style dense -> bidiagonal, then BDSQR (or BDSDC when explicitly selected), then ORMBR-style back-transforms for full U and V^H.

Precondition

in-order Queue

Exceptions

batchlas::invalid_argument	for an out-of-order Queue or mis-shaped outputs
batchlas::unsupported	for complex input (use the Hermitian overload)

gesvd_blocked() [2/2]

template<Backend B, typename T>

Event batchlas::gesvd_blocked	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a_in,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		Uplo	hermitian_uplo,
		const Span< std::byte > &	ws,
		Span< int32_t >	info = Span< int32_t >() )

gesvd_blocked() for Hermitian input: only the hermitian_uplo triangle of a_in is read.

gesvd_blocked_buffer_size() [1/2]

template<Backend B, typename T>

size_t batchlas::gesvd_blocked_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh )

Required workspace, in bytes, for gesvd_blocked().

gesvd_blocked_buffer_size() [2/2]

template<Backend B, typename T>

size_t batchlas::gesvd_blocked_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		Uplo	hermitian_uplo )

Required workspace, in bytes, for the Hermitian gesvd_blocked().

gesvd_cta() [1/2]

template<Backend B, typename T>

Event batchlas::gesvd_cta	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a_in,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		const Span< std::byte > &	ws,
		Span< int32_t >	info = Span< int32_t >() )

CTA-oriented native SVD for small matrices, max(m, n) <= 32: bidiagonal reduction, then the eigenproblem of \( B^T B \) (the normal equations).

Real input, rectangular either way.

The normal equations square the condition number; use gesvdj_cta() where small singular values matter.

Precondition

in-order Queue

Exceptions

batchlas::invalid_argument	for max(m, n) > 32, a genuinely Thin U or V^H, or an out-of-order Queue
batchlas::unsupported	for complex input (use the Hermitian overload)

See also

GESVD

gesvd_cta() [2/2]

template<Backend B, typename T>

Event batchlas::gesvd_cta	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a_in,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		Uplo	hermitian_uplo,
		const Span< std::byte > &	ws,
		Span< int32_t >	info = Span< int32_t >() )

gesvd_cta() for Hermitian input: only the hermitian_uplo triangle is read.

gesvd_cta_buffer_size() [1/2]

template<Backend B, typename T>

size_t batchlas::gesvd_cta_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh )

Required workspace, in bytes, for gesvd_cta().

gesvd_cta_buffer_size() [2/2]

template<Backend B, typename T>

size_t batchlas::gesvd_cta_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		Uplo	hermitian_uplo )

Required workspace, in bytes, for the Hermitian gesvd_cta().

gesvdj_cta()

template<Backend B, typename T>

Event batchlas::gesvdj_cta	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a_in,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		const Span< std::byte > &	ws = Span< std::byte >(),
		GesvdjParams< T >	params = GesvdjParams< T >(),
		Span< int32_t >	info = Span< int32_t >() )

One-sided (Hestenes) Jacobi SVD for batches of small matrices.

Computes \( A = U \operatorname{diag}(s) V^H \) with high RELATIVE accuracy: the singular-value error is governed by the condition number of the column-equilibrated matrix rather than of A, so graded input keeps its small values; gesvd_cta() / gesvd_blocked() do not.

max(m, n) <= 64 (32 for complex<double> with vectors), real and complex, rectangular either way; needs sub-group size 32. A is DESTROYED, singular values are returned descending, and no workspace is required.

Exceptions

batchlas::invalid_argument	for a size above the cap or undersized outputs
batchlas::unsupported	if the device has no sub-group size 32

See also

gesvd design, GESVD

gesvdj_cta_buffer_size()

template<Backend B, typename T>

size_t batchlas::gesvdj_cta_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	a,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	u_out,
		const MatrixView< T, MatrixFormat::Dense > &	vh_out,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		GesvdjParams< T >	params = GesvdjParams< T >() )

Required workspace, in bytes, for gesvdj_cta() (currently 0: all state is in local memory).

gesvd_buffer_size() [1/2]

template<Backend B, typename T>

size_t batchlas::gesvd_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	A,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	U,
		const MatrixView< T, MatrixFormat::Dense > &	Vh,
		SvdVectors	jobu,
		SvdVectors	jobvh )

Workspace, in bytes, that gesvd() needs for the same arguments.

Canonicalises jobu/jobvh and makes the same kernel choice as the call does, so the size is for the kernel that will run. info does not affect it.

gesvd_buffer_size() [2/2]

template<Backend B, typename T>

size_t batchlas::gesvd_buffer_size	(	Queue &	ctx,
		const MatrixView< T, MatrixFormat::Dense > &	A,
		Span< typename base_type< T >::type >	singular_values,
		const MatrixView< T, MatrixFormat::Dense > &	U,
		const MatrixView< T, MatrixFormat::Dense > &	Vh,
		SvdVectors	jobu,
		SvdVectors	jobvh,
		Uplo	hermitian_uplo )

Workspace, in bytes, for the Hermitian gesvd().