naginterfaces.library.lapackeig.dtgsja

naginterfaces.library.lapackeig.dtgsja(jobu, jobv, jobq, k, l, a, b, tola, tolb, u, v, q)

dtgsja computes the generalized singular value decomposition (GSVD) of two real upper trapezoidal matrices $A$ and $B$, where $A$ is an $m\times n$ matrix and $B$ is a $p\times n$ matrix.

$A$ and $B$ are assumed to be in the form returned by dggsvp3().

For full information please refer to the NAG Library document for f08ye

https://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/f08/f08yef.html

Parameters

jobustr, length 1

If $jobu=\text{'U'}$, $u$ must contain an orthogonal matrix $U_1$ on entry, and the product $U_{1}U$ is returned.

If $jobu=\text{'I'}$, $u$ is initialized to the unit matrix, and the orthogonal matrix $U$ is returned.

If $jobu=\text{'N'}$, $U$ is not computed.

jobvstr, length 1

If $jobv=\text{'V'}$, $v$ must contain an orthogonal matrix $V_1$ on entry, and the product $V_{1}V$ is returned.

If $jobv=\text{'I'}$, $v$ is initialized to the unit matrix, and the orthogonal matrix $V$ is returned.

If $jobv=\text{'N'}$, $V$ is not computed.

jobqstr, length 1

If $jobq=\text{'Q'}$, $q$ must contain an orthogonal matrix $Q_1$ on entry, and the product $Q_{1}Q$ is returned.

If $jobq=\text{'I'}$, $q$ is initialized to the unit matrix, and the orthogonal matrix $Q$ is returned.

If $jobq=\text{'N'}$, $Q$ is not computed.

kint

$k$ and $l$ specify the sizes, $k$ and $l$, of the subblocks of $A$ and $B$, whose GSVD is to be computed by dtgsja.

lint

$k$ and $l$ specify the sizes, $k$ and $l$, of the subblocks of $A$ and $B$, whose GSVD is to be computed by dtgsja.

afloat, array-like, shape $(m,n)$

The $m\times n$ matrix $A$.

bfloat, array-like, shape $(p,n)$

The $p\times n$ matrix $B$.

tolafloat

$tola$ and $tolb$ are the convergence criteria for the Jacobi–Kogbetliantz iteration procedure. Generally, they should be the same as used in the preprocessing step performed by dggsvp3(), say

$$\begin{array}{c}\begin{array}{c}tola=max(m,n)\parallel A\parallel ϵ\text{,}\\ tolb=max(p,n)\parallel B\parallel ϵ\text{,}\end{array}\end{array}$$

where $ϵ$ is the machine precision.

tolbfloat

$tola$ and $tolb$ are the convergence criteria for the Jacobi–Kogbetliantz iteration procedure. Generally, they should be the same as used in the preprocessing step performed by dggsvp3(), say

$$\begin{array}{c}\begin{array}{c}tola=max(m,n)\parallel A\parallel ϵ\text{,}\\ tolb=max(p,n)\parallel B\parallel ϵ\text{,}\end{array}\end{array}$$

where $ϵ$ is the machine precision.

ufloat, array-like, shape $(:,:)$

Note: the required extent for this argument in dimension 1 is determined as follows: if $jobu\text{in}(\text{'U'},\text{'I'})$: $m$; otherwise: $1$.

Note: the required extent for this argument in dimension 2 is determined as follows: if $jobu\text{in}(\text{'U'},\text{'I'})$: $m$; otherwise: $1$.

If $jobu=\text{'U'}$, $u$ must contain an $m\times m$ matrix $U_1$ (usually the orthogonal matrix returned by dggsvp3()).

vfloat, array-like, shape $(:,:)$

Note: the required extent for this argument in dimension 1 is determined as follows: if $jobv\text{in}(\text{'V'},\text{'I'})$: $p$; otherwise: $1$.

Note: the required extent for this argument in dimension 2 is determined as follows: if $jobv\text{in}(\text{'V'},\text{'I'})$: $p$; otherwise: $1$.

If $jobv=\text{'V'}$, $v$ must contain an $p\times p$ matrix $V_1$ (usually the orthogonal matrix returned by dggsvp3()).

qfloat, array-like, shape $(:,:)$

Note: the required extent for this argument in dimension 1 is determined as follows: if $jobq\text{in}(\text{'Q'},\text{'I'})$: $n$; otherwise: $1$.

Note: the required extent for this argument in dimension 2 is determined as follows: if $jobq\text{in}(\text{'Q'},\text{'I'})$: $n$; otherwise: $1$.

If $jobq=\text{'Q'}$, $q$ must contain an $n\times n$ matrix $Q_1$ (usually the orthogonal matrix returned by dggsvp3()).

Returns

afloat, ndarray, shape $(m,n)$

If $m-k-l\ge 0$, $a[0:k+l,n-k-l:n]$ contains the $(k+l)\times (k+l)$ upper triangular matrix $R$.

If $m-k-l<0$, $a[0:m,n-k-l:n]$ contains the first $m$ rows of the $(k+l)\times (k+l)$ upper triangular matrix $R$, and the submatrix $R_{33}$ is returned in $b[m-k:l,n+m-k-l:n]$.

bfloat, ndarray, shape $(p,n)$

If $m-k-l<0$, $b[m-k:l,n+m-k-l:n]$ contains the submatrix $R_{33}$ of $R$.

alphafloat, ndarray, shape $\left(n\right)$

See the description of $beta$.

betafloat, ndarray, shape $\left(n\right)$

$alpha$ and $beta$ contain the generalized singular value pairs of $A$ and $B$;

$alpha\left[\text{i}\right]=1$, $beta\left[\text{i}\right]=0$, for $\text{i}=0,\dots ,k-1$, and

if $m-k-l\ge 0$, $alpha\left[\text{i}\right]={\alpha }_{\text{i}}$, $beta\left[\text{i}\right]={\beta }_{\text{i}}$, for $\text{i}=k,\dots ,k+l-1$, or

if $m-k-l<0$, $alpha\left[\text{i}\right]={\alpha }_{\text{i}}$, $beta\left[\text{i}\right]={\beta }_{\text{i}}$, for $\text{i}=k,\dots ,m-1$ and $alpha\left[\text{i}\right]=0$, $beta\left[\text{i}\right]=1$, for $\text{i}=m,\dots ,k+l-1$.

Furthermore, if $k+l<n$, $alpha\left[\text{i}\right]=beta\left[\text{i}\right]=0$, for $\text{i}=k+l,\dots ,n-1$.

ufloat, ndarray, shape $(:,:)$

If $jobu=\text{'U'}$, $u$ contains the product $U_{1}U$.

If $jobu=\text{'I'}$, $u$ contains the orthogonal matrix $U$.

If $jobu=\text{'N'}$, $u$ is not referenced.

vfloat, ndarray, shape $(:,:)$

If $jobv=\text{'I'}$, $v$ contains the orthogonal matrix $V$.

If $jobv=\text{'V'}$, $v$ contains the product $V_{1}V$.

If $jobv=\text{'N'}$, $v$ is not referenced.

qfloat, ndarray, shape $(:,:)$

If $jobq=\text{'I'}$, $q$ contains the orthogonal matrix $Q$.

If $jobq=\text{'Q'}$, $q$ contains the product $Q_{1}Q$.

If $jobq=\text{'N'}$, $q$ is not referenced.

ncycleint

The number of cycles required for convergence.

Raises

NagValueError

(errno $-1$)

On entry, error in parameter $jobu$.

Constraint: $jobu=\text{'U'}$, $\text{'I'}$ or $\text{'N'}$.

(errno $-2$)

On entry, error in parameter $jobv$.

Constraint: $jobv=\text{'V'}$, $\text{'I'}$ or $\text{'N'}$.

(errno $-3$)

On entry, error in parameter $jobq$.

Constraint: $jobq=\text{'Q'}$, $\text{'I'}$ or $\text{'N'}$.

(errno $-4$)

On entry, error in parameter $\text{m}$.

Constraint: $m\ge 0$.

(errno $-5$)

On entry, error in parameter $\text{p}$.

Constraint: $p\ge 0$.

(errno $-6$)

On entry, error in parameter $\text{n}$.

Constraint: $n\ge 0$.

(errno $1$)

The procedure does not converge after $40$ cycles.

Notes

dtgsja computes the GSVD of the matrices $A$ and $B$ which are assumed to have the form as returned by dggsvp3()

$$\begin{array}{c}\begin{array}{cc}A=& \{\begin{array}{c}\begin{array}{c}\\ k\\ l\\ m-k-l\end{array}\left(\begin{array}{ccc}n-k-l& k& l\\ 0& A_{12}& A_{13}\\ 0& 0& A_{23}\\ 0& 0& 0\end{array}\right)\text{, if}m-k-l\ge 0;\\ \begin{array}{c}\\ k\\ m-k\end{array}\left(\begin{array}{ccc}n-k-l& k& l\\ 0& A_{12}& A_{13}\\ 0& 0& A_{23}\end{array}\right)\text{, if}m-k-l<0;\end{array}\\ & \\ B=& \begin{array}{c}\\ l\\ p-l\end{array}\left(\begin{array}{ccc}n-k-l& k& l\\ 0& 0& B_{13}\\ 0& 0& 0\end{array}\right)\text{,}\end{array}\end{array}$$

where the $k\times k$ matrix $A_{12}$ and the $l\times l$ matrix $B_{13}$ are nonsingular upper triangular, $A_{23}$ is $l\times l$ upper triangular if $m-k-l\ge 0$ and is $(m-k)\times l$ upper trapezoidal otherwise.

dtgsja computes orthogonal matrices $Q$, $U$ and $V$, diagonal matrices $D_1$ and $D_2$, and an upper triangular matrix $R$ such that

$$U^{T}AQ=D_1\left(\begin{array}{cc}0& R\end{array}\right)\text{,}{V}^{T}BQ=D_2\left(\begin{array}{cc}0& R\end{array}\right)\text{.}$$

Optionally $Q$, $U$ and $V$ may or may not be computed, or they may be premultiplied by matrices $Q_1$, $U_1$ and $V_1$ respectively.

If $(m-k-l)\ge 0$ then $D_1$, $D_2$ and $R$ have the form

$$\begin{array}{c}{D}_1=\begin{array}{c}\\ k\\ l\\ m-k-l\end{array}\left(\begin{array}{cc}k& l\\ I& 0\\ 0& C\\ 0& 0\end{array}\right)\text{,}\end{array}$$

$$\begin{array}{c}{D}_2=\begin{array}{c}\\ l\\ p-l\end{array}\left(\begin{array}{cc}k& l\\ 0& S\\ 0& 0\end{array}\right)\text{,}\end{array}$$

$$\begin{array}{c}R=\begin{array}{c}\\ k\\ l\end{array}\left(\begin{array}{cc}k& l\\ R_{11}& R_{12}\\ 0& R_{22}\end{array}\right)\text{,}\end{array}$$

where $C=diag({\alpha }_{k+1},,,\dots ,,,{\alpha }_{k+l})\text{,}S=diag({\beta }_{k+1},,,\dots ,,,{\beta }_{k+l})$.

If $(m-k-l)<0$ then $D_1$, $D_2$ and $R$ have the form

$$\begin{array}{c}{D}_1=\begin{array}{c}\\ k\\ m-k\end{array}\left(\begin{array}{ccc}k& m-k& k+l-m\\ I& 0& 0\\ 0& C& 0\end{array}\right)\text{,}\end{array}$$

$$\begin{array}{c}{D}_2=\begin{array}{c}\\ m-k\\ k+l-m\\ p-l\end{array}\left(\begin{array}{ccc}k& m-k& k+l-m\\ 0& S& 0\\ 0& 0& I\\ 0& 0& 0\end{array}\right)\text{,}\end{array}$$

$$\begin{array}{c}R=\begin{array}{c}\\ k\\ m-k\\ k+l-m\end{array}\left(\begin{array}{ccc}k& m-k& k+l-m\\ R_{11}& R_{12}& R_{13}\\ 0& R_{22}& R_{23}\\ 0& 0& R_{33}\end{array}\right)\text{,}\end{array}$$

where $C=diag({\alpha }_{k+1},,,\dots ,,,{\alpha }_m)\text{,}S=diag({\beta }_{k+1},,,\dots ,,,{\beta }_m)$.

In both cases the diagonal matrix $C$ has non-negative diagonal elements, the diagonal matrix $S$ has positive diagonal elements, so that $S$ is nonsingular, and $C^2+S^2=1$. See Section 2.3.5.3 of Anderson et al. (1999) for further information.

References

Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, LAPACK Users’ Guide, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug

Golub, G H and Van Loan, C F, 1996, Matrix Computations, (3rd Edition), Johns Hopkins University Press, Baltimore