Solve symmetric eigenproblems and compute singular value decompositions. More...

Data Structures
struct	_tridiag
	Real tridiagonal matrix , represented by vectors containing the diagonal, sub- and superdiagonal. More...

Typedefs
typedef struct _tridiag	tridiag
	Tridiagonal matrix.

typedef tridiag *	ptridiag
	Pointer to a tridiag object.

typedef const tridiag *	pctridiag
	Pointer to a constant tridiag object.

Functions
ptridiag	init_tridiag (ptridiag T, uint size)
	Initialize a tridiag object. More...

ptridiag	init_sub_tridiag (ptridiag T, ptridiag src, uint size, uint off)
	Initialize a tridiag object to represent a submatrix. More...

ptridiag	init_vec_tridiag (ptridiag T, prealavector src, uint size)
	Initialize a tridiag object using pre-allocated storage. More...

void	uninit_tridiag (ptridiag T)
	Uninitialize a tridiag object. More...

ptridiag	new_tridiag (uint size)
	Create a new tridiag object. More...

void	del_tridiag (ptridiag T)
	Delete a tridiag object. More...

void	copy_tridiag (pctridiag T, ptridiag Tcopy)
	Copy a matrix into another matrix, $T_{\rm copy} \gets T$ . More...

real	check_tridiag (pctridiag T, pcamatrix Ts)
	Compute the Frobenius norm $\\|T - T_s\\|_F$ of the difference of a tridiagonal matrix and a general matrix . More...

real	check_lower_tridiag (pctridiag T, pcamatrix Ts)
	Compute the Frobenius norm $\\|T - T_s\\|_F$ of the difference of a lower bidiagonal matrix and a general matrix . More...

void	diageval_tridiag_amatrix (field alpha, bool atrans, pctridiag a, bool xtrans, pamatrix x)
	Multiply a matrix by the diagonal part of a tridiagonal matrix, $X \gets \alpha A X$ . More...

void	lowereval_tridiag_amatrix (field alpha, bool atrans, pctridiag a, bool xtrans, pamatrix x)
	Multiply a matrix by the lower bidiagonal part of a tridiagonal matrix, $X \gets \alpha A X$ . More...

void	qrstep_tridiag (ptridiag T, field shift, pamatrix Q)
	Perform one implicit QR step on a self-adjoint tridiagonal matrix. More...

uint	sb_muleig_tridiag (ptridiag T, pamatrix Q, uint maxiter)
	Solve a self-adjoint tridiagonal eigenproblem, self-made implementation without BLAS or LAPACK. More...

uint	muleig_tridiag (ptridiag T, pamatrix Q)
	Solve a self-adjoint tridiagonal eigenproblem. More...

uint	eig_tridiag (ptridiag T, pamatrix Q)
	Solve a self-adjoint tridiagonal eigenproblem. More...

void	sb_tridiagonalize_amatrix (pamatrix A, ptridiag T, pamatrix Q)
	Hessenberg tridiagonalization of a self-adjoint matrix, self-made implementation without BLAS or LAPACK. More...

void	tridiagonalize_amatrix (pamatrix A, ptridiag T, pamatrix Q)
	Hessenberg tridiagonalization of a self-adjoint matrix. More...

uint	sb_eig_amatrix (pamatrix A, prealavector lambda, pamatrix Q, uint maxiter)
	Solve a self-adjoint eigenproblem, self-made implementation without BLAS or LAPACK. More...

uint	eig_amatrix (pamatrix A, prealavector lambda, pamatrix Q)
	Solve a self-adjoint eigenproblem, $A e = \lambda e$ . More...

uint	geig_amatrix (pamatrix A, pamatrix M, prealavector lambda, pamatrix Q)
	Solve self-adjoint generalized eigenproblem, $A e = \lambda M e$ . More...

void	svdstep_tridiag (ptridiag T, field shift, pamatrix U, pamatrix Vt)
	Perform one step of the Golub-Kahan iteration on a lower bidiagonal matrix. More...

uint	sb_mulsvd_tridiag (ptridiag T, pamatrix U, pamatrix Vt, uint maxiter)
	Compute the SVD of a bidiagonal matrix, self-made implementation without BLAS or LAPACK. More...

uint	mulsvd_tridiag (ptridiag T, pamatrix U, pamatrix Vt)
	Compute the SVD of a bidiagonal matrix. More...

uint	svd_tridiag (ptridiag T, pamatrix U, pamatrix Vt)
	Compute the SVD of a bidiagonal matrix. More...

void	sb_bidiagonalize_amatrix (pamatrix A, ptridiag T, pamatrix U, pamatrix Vt)
	Golub-Kahan bidiagonalization of a matrix, self-made implementation without BLAS or LAPACK. More...

void	bidiagonalize_amatrix (pamatrix A, ptridiag T, pamatrix U, pamatrix Vt)
	Golub-Kahan bidiagonalization of a matrix. More...

uint	sb_svd_amatrix (pamatrix A, prealavector sigma, pamatrix U, pamatrix Vt, uint maxiter)
	Compute the SVD of a matrix, self-made implementation without BLAS or LAPACK. More...

uint	svd_amatrix (pamatrix A, prealavector sigma, pamatrix U, pamatrix Vt)
	Compute the SVD of a matrix. More...

Detailed Description

Solve symmetric eigenproblems and compute singular value decompositions.

Function Documentation

void bidiagonalize_amatrix	(	pamatrix	A,
		ptridiag	T,
		pamatrix	U,
		pamatrix	Vt
	)

Golub-Kahan bidiagonalization of a matrix.

Compute unitary matrices $U$ and $V$ and a lower bidiagonal matrix $T$ such that $A = U T V^*$ .

Remarks: This function is only intended as an intermediate step of svd_amatrix on systems that do not offer a LAPACK library.

Parameters

A	Matrix , will be overwritten by the function.
T	Bidiagonal matrix .
U	If `U!=0`, this matrix will be filled with the unitary transformation .
Vt	If `Vt!=0`, this matrix will be filled with the adjoint unitary transformation .

real check_lower_tridiag	(	pctridiag	T,
		pcamatrix	Ts
	)

Compute the Frobenius norm $\|T - T_s\|_F$ of the difference of a lower bidiagonal matrix $T$ and a general matrix $T_s$ .

Parameters

T	First matrix, represented as tridiag object, where the superdiagonal `T->u` is ignored.
Ts	Second matrix, represented as amatrix object.

Returns: Frobenius norm $\|T - T_s\|_F$ .

real check_tridiag	(	pctridiag	T,
		pcamatrix	Ts
	)

Compute the Frobenius norm $\|T - T_s\|_F$ of the difference of a tridiagonal matrix $T$ and a general matrix $T_s$ .

Parameters

T	First matrix, represented as tridiag object.
Ts	Second matrix, represented as amatrix object.

Returns: Frobenius norm $\|T - T_s\|_F$ .

void copy_tridiag	(	pctridiag	T,
		ptridiag	Tcopy
	)

Copy a matrix into another matrix, $T_{\rm copy} \gets T$ .

Parameters

T	Source matrix .
Tcopy	Target matrix $T_{\rm copy}$ .

void del_tridiag ( ptridiag T )

Delete a tridiag object.

Releases the storage corresponding to the object.

Attention: Make sure that there are no submatrix objects referring to this object, since they will otherwise keep using pointers to invalid storage.

Parameters

T	Object to be deleted.

void diageval_tridiag_amatrix	(	field	alpha,
		bool	atrans,
		pctridiag	a,
		bool	xtrans,
		pamatrix	x
	)

Multiply a matrix $X$ by the diagonal part $A$ of a tridiagonal matrix, $X \gets \alpha A X$ .

Parameters

alpha	Scaling factor $\alpha$ .
atrans	Set if is to be used instead of .
a	Matrix , only its diagonal will be used.
xtrans	Set if is to be used instead of .
x	Matrix , will be overwritten by the result.

uint eig_amatrix	(	pamatrix	A,
		prealavector	lambda,
		pamatrix	Q
	)

Solve a self-adjoint eigenproblem, $A e = \lambda e$ .

Parameters

A	Matrix , will be overwritten by the function.
lambda	Eigenvalues of , in ascending order.
Q	If `Q!=0`, the columns of this matrix will contain an orthonormal basis of eigenvectors corresponding to the eigenvalues in `lambda`.

Returns: Zero on successful completion, non-zero if the QR iteration did not converge.

uint eig_tridiag	(	ptridiag	T,
		pamatrix	Q
	)

Solve a self-adjoint tridiagonal eigenproblem.

Eigenvalues will be on the diagonal of T, in ascending order.

Parameters

T	Matrix , will be overwritten by a diagonal matrix such that with a unitary matrix .
Q	If `Q!=0`, this matrix will be filled with an orthonormal basis of eigenvectors.

Returns: Zero on successful completion, non-zero if the QR iteration did not converge.

uint geig_amatrix	(	pamatrix	A,
		pamatrix	M,
		prealavector	lambda,
		pamatrix	Q
	)

Solve self-adjoint generalized eigenproblem, $A e = \lambda M e$ .

Parameters

A	Matrix , will be overwritten by the function.
M	Matrix , will be overwritten by the Cholesky factorization of .
lambda	Eigenvalues of , in ascending order.
Q	If `Q!=0`, the columns of this matrix will contain an -orthonormal basis of eigenvectors corresponding to the eigenvalues in `lambda`.

Returns: Zero on successful completion, non-zero if the QR iteration did not converge.

ptridiag init_sub_tridiag	(	ptridiag	T,
		ptridiag	src,
		uint	size,
		uint	off
	)

Initialize a tridiag object to represent a submatrix.

Sets up the components of the object and uses part of the storage of another amatrix for the coefficient vectors, leading to a new matrix representing a submatrix of the source.

Remarks: Should always be matched by a call to uninit_tridiag that will not release the coefficient storage.

Parameters

T	Object to be initialized.
src	Source matrix .
size	Dimension of submatrix.
off	Row and column offset.

Returns: Initialized tridiag object.

ptridiag init_tridiag	(	ptridiag	T,
		uint	size
	)

Initialize a tridiag object.

Sets up the components of the object and allocates storage for the coefficient vectors.

Remarks: Should always be matched by a call to uninit_tridiag.

Parameters

T	Object to be initialized.
size	Dimension of the new matrix.

Returns: Initialized tridiag object.

ptridiag init_vec_tridiag	(	ptridiag	T,
		prealavector	src,
		uint	size
	)

Initialize a tridiag object using pre-allocated storage.

Sets up the components of the object and uses the storage of a vector of at least dimension 3*dim-2 for the coefficient vectors.

Remarks: Should always be matched by a call to uninit_tridiag that will not release the coefficient storage.

Parameters

T	Object to be initialized.
src	Source vector.
size	Dimension of the new matrix.

Returns: Initialized tridiag object.

void lowereval_tridiag_amatrix	(	field	alpha,
		bool	atrans,
		pctridiag	a,
		bool	xtrans,
		pamatrix	x
	)

Multiply a matrix $X$ by the lower bidiagonal part $A$ of a tridiagonal matrix, $X \gets \alpha A X$ .

Parameters

alpha	Scaling factor $\alpha$ .
atrans	Set if is to be used instead of .
a	Matrix , only its diagonal and lower subdiagonal will be used.
xtrans	Set if is to be used instead of .
x	Matrix , will be overwritten by the result.

uint muleig_tridiag	(	ptridiag	T,
		pamatrix	Q
	)

Solve a self-adjoint tridiagonal eigenproblem.

Eigenvalues will be on the diagonal of T, in ascending order.

Parameters

T	Matrix , will be overwritten by a diagonal matrix such that $T = \widehat{Q} D \widehat{Q}^*$ with a unitary matrix $\widehat{Q}$ .
Q	Unitary matrix used to accumulate transformations, will be overwritten by $Q \widehat{Q}$ . If `Q==0`, the transformations will not be accumulated.

Returns: Zero on successful completion, non-zero if the QR iteration did not converge.

uint mulsvd_tridiag	(	ptridiag	T,
		pamatrix	U,
		pamatrix	Vt
	)

Compute the SVD of a bidiagonal matrix.

Compute the factorization $T = U \Sigma V^*$ with a real diagonal matrix $\Sigma$ and unitary matrices $U$ and $V$ .

Parameters

T	Bidiagonal matrix , will be overwritten by a diagonal matrix such that $T = \widehat{U} D \widehat{V}^*$ with unitary matrices $\widehat{U}$ and $\widehat{V}$ .
U	Unitary matrix used to accumulate transformations, will be overwritten by $U \widehat{U}$ . If `U==0`, the transformations will not be accumulated.
Vt	Adjoint of unitary matrix used to accumulate transformations, will be overwritten by $\widehat{V}^* V_t = (V \widehat{V})^*$ . If `Vt==0`, the transformations will not be accumulated.

Returns: Zero on successful completion, non-zero if the QR iteration did not converge.

ptridiag new_tridiag ( uint size )

Create a new tridiag object.

Allocates storage for the object and sets up its components.

Remarks: Should always be matched by a call to del_tridiag.

Parameters

size	Dimension of the new matrix.

Returns: New tridiag object.

void qrstep_tridiag	(	ptridiag	T,
		field	shift,
		pamatrix	Q
	)

Perform one implicit QR step on a self-adjoint tridiagonal matrix.

Parameters

T	Matrix , will be overwritten by $\widehat{Q}^* T \widehat{Q}$ , where $\widehat{Q}$ is the unitary transformation corresponding to the QR step.
shift	Shift parameter $\mu$ , the matrix is shifted to $T-\mu I$ before computing the QR factorization.
Q	Can be used to accumulate transformations, will be overwritten by $Q \widehat{Q}$ . If `Q==0`, the transformations will not be accumulated.

void sb_bidiagonalize_amatrix	(	pamatrix	A,
		ptridiag	T,
		pamatrix	U,
		pamatrix	Vt
	)

Golub-Kahan bidiagonalization of a matrix, self-made implementation without BLAS or LAPACK.

Compute unitary matrices $U$ and $V$ and a lower bidiagonal matrix $T$ such that $A = U T V^*$ .

Remarks: This function is only intended as an intermediate step of svd_amatrix on systems that do not offer a LAPACK library.

Parameters

A	Matrix , will be overwritten by the function.
T	Bidiagonal matrix .
U	If `U!=0`, this matrix will be filled with the unitary transformation .
Vt	If `Vt!=0`, this matrix will be filled with the adjoint unitary transformation .

uint sb_eig_amatrix	(	pamatrix	A,
		prealavector	lambda,
		pamatrix	Q,
		uint	maxiter
	)

Solve a self-adjoint eigenproblem, self-made implementation without BLAS or LAPACK.

Remarks: This is a very simple implementation using Wilkinson shifts. It is only intended as a stand-in on systems that do not offer a LAPACK library. Whenever possible, eig_amatrix should be used instead.

Parameters

A	Matrix , will be overwritten by the function.
lambda	Eigenvalues of , in ascending order.
Q	If `Q!=0`, the columns of this matrix will contain an orthonormal basis of eigenvectors corresponding to the eigenvalues in `lambda`.
maxiter	Upper bound for the number of QR steps.

Returns: Number of QR steps.

uint sb_muleig_tridiag	(	ptridiag	T,
		pamatrix	Q,
		uint	maxiter
	)

Solve a self-adjoint tridiagonal eigenproblem, self-made implementation without BLAS or LAPACK.

Eigenvalues will be on the diagonal of T, in ascending order.

Remarks: This is a very simple implementation using Wilkinson shifts. It is only intended as a stand-in on systems that do not offer a LAPACK library. Whenever possible, muleig_tridiag should be used instead.

Parameters

T	Matrix , will be overwritten by a diagonal matrix such that $T = \widehat{Q} D \widehat{Q}^*$ with a unitary matrix $\widehat{Q}$ .
Q	Unitary matrix used to accumulate transformations, will be overwritten by $Q \widehat{Q}$ . If `Q==0`, the transformations will not be accumulated.
maxiter	Upper bound for the number of QR steps.

Returns: Number of QR steps.

uint sb_mulsvd_tridiag	(	ptridiag	T,
		pamatrix	U,
		pamatrix	Vt,
		uint	maxiter
	)

Compute the SVD of a bidiagonal matrix, self-made implementation without BLAS or LAPACK.

Singular values will be on the diagonal of T, in descending order.

Remarks: This is a very simple implementation using Wilkinson shifts. It is only intended as a stand-in on systems that do not offer a LAPACK library. Whenever possible, mulsvd_tridiag should be used instead.

Parameters

T	Bidiagonal matrix , will be overwritten by a diagonal matrix such that $T = \widehat{U} D \widehat{V}^*$ with unitary matrices $\widehat{U}$ and $\widehat{V}$ .
U	Unitary matrix used to accumulate transformations, will be overwritten by $U \widehat{U}$ . If `U==0`, the transformations will not be accumulated.
Vt	Adjoint of unitary matrix used to accumulate transformations, will be overwritten by $\widehat{V}^* V_t = (V \widehat{V})^*$ . If `Vt==0`, the transformations will not be accumulated.
maxiter	Upper bound for the number of Golub-Kahan steps.

Returns: Number of Golub-Kahan steps.

uint sb_svd_amatrix	(	pamatrix	A,
		prealavector	sigma,
		pamatrix	U,
		pamatrix	Vt,
		uint	maxiter
	)

Compute the SVD of a matrix, self-made implementation without BLAS or LAPACK.

Compute the factorization $A = U \Sigma V^*$ with a real diagonal matrix $\Sigma=\mathop{\operatorname{diag}}(\sigma_1,\ldots,\sigma_k)$ , $\sigma_1\geq\sigma_2\geq\ldots$ and unitary matrices $U$ and $V$ .

Remarks: This function is only intended as an intermediate step of svd_amatrix on systems that do not offer a LAPACK library.

Parameters

A	Matrix , will be overwritten by the function.
sigma	Singular values $\sigma_1,\sigma_2,\ldots$ .
U	If `U!=0`, this matrix will be filled with the unitary transformation .
Vt	If `Vt!=0`, this matrix will be filled with the adjoint unitary transformation .
maxiter	Upper bound for the number of Golub-Kahan steps.

Returns: Number of Golub-Kahan steps.

void sb_tridiagonalize_amatrix	(	pamatrix	A,
		ptridiag	T,
		pamatrix	Q
	)

Hessenberg tridiagonalization of a self-adjoint matrix, self-made implementation without BLAS or LAPACK.

Compute a unitary matrix $Q$ and a tridiagonal self-adjoint matrix $T$ such that $A = Q T Q^*$ .

Remarks: This function is only intended as an intermediate step of eig_amatrix on systems that do not offer a LAPACK library.

Parameters

A	Matrix , will be overwritten by the function.
T	Tridiagonal matrix .
Q	If `Q!=0`, this matrix will be filled with the unitary transformation .

uint svd_amatrix	(	pamatrix	A,
		prealavector	sigma,
		pamatrix	U,
		pamatrix	Vt
	)

Compute the SVD of a matrix.

Compute the factorization $A = U \Sigma V^*$ with a real diagonal matrix $\Sigma=\mathop{\operatorname{diag}}(\sigma_1,\ldots,\sigma_k)$ , $\sigma_1\geq\sigma_2\geq\ldots$ and unitary matrices $U$ and $V$ .

Parameters

A	Matrix , will be overwritten by the function.
sigma	Singular values $\sigma_1,\sigma_2,\ldots$ .
U	If `U!=0`, this matrix will be filled with the unitary transformation .
Vt	If `Vt!=0`, this matrix will be filled with the adjoint unitary transformation .

Returns: Zero on successful completion, non-zero if the Golub-Kahan iteration did not converge.

uint svd_tridiag	(	ptridiag	T,
		pamatrix	U,
		pamatrix	Vt
	)

Compute the SVD of a bidiagonal matrix.

Singular values will be on the diagonal of T, in descending order.

Parameters

T	Bidiagonal matrix , will be overwritten by a diagonal matrix such that $T = \widehat{U} D \widehat{V}^*$ with unitary matrices $\widehat{U}$ and $\widehat{V}$ .
U	If `U!=0`, the columns will filled by the leading left singular vectors.
Vt	If `Vt!=0`, the rows will be filled by the conjugated leading right singular vectors.

Returns: Zero on successful completion, non-zero if the QR iteration did not converge.

void svdstep_tridiag	(	ptridiag	T,
		field	shift,
		pamatrix	U,
		pamatrix	Vt
	)

Perform one step of the Golub-Kahan iteration on a lower bidiagonal matrix.

Parameters

T	Bidiagonal matrix , will be overwritten by $\widehat{U}^* T \widehat{V}$ .
shift	Shift parameter $\mu$ , the matrix is shifted to $T T^* - \mu I$ before computing the QR factorization.
U	Can be used to accumulate transformations, will be overwritten by $U \widehat{U}$ . If `U==0`, the transformations will not be accumulated.
Vt	Can be used to accumulate transformations, will be overwritten by $\widehat{V}^* V_t$ . If `Vt==0`, the transformations will not be accumulated.

void tridiagonalize_amatrix	(	pamatrix	A,
		ptridiag	T,
		pamatrix	Q
	)

Hessenberg tridiagonalization of a self-adjoint matrix.

Compute a unitary matrix $Q$ and a tridiagonal self-adjoint matrix $T$ such that $A = Q T Q^*$ .

Remarks: This function is only intended as an intermediate step of eig_amatrix on systems that do not offer a LAPACK library.

Parameters

A	Matrix , will be overwritten by the function.
T	Tridiagonal matrix .
Q	If `Q!=0`, this matrix will be filled with the unitary transformation .

void uninit_tridiag ( ptridiag T )

Uninitialize a tridiag object.

Invalidates pointers, freeing corresponding storage if owner==0, and prepares the object for deletion.

Parameters

T	Object to be uninitialized.

Data Structures

Typedefs

Functions

Detailed Description

Function Documentation