F08 Chapter Contents
F08 Chapter Introduction
NAG Library Manual

# NAG Library Routine DocumentF08YSF (ZTGSJA)

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

## 1  Purpose

F08YSF (ZTGSJA) computes the generalized singular value decomposition (GSVD) of two complex upper trapezoidal matrices $A$ and $B$, where $A$ is an $m$ by $n$ matrix and $B$ is a $p$ by $n$ matrix.
$A$ and $B$ are assumed to be in the form returned by F08VSF (ZGGSVP).

## 2  Specification

 SUBROUTINE F08YSF ( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, NCYCLE, INFO)
 INTEGER M, P, N, K, L, LDA, LDB, LDU, LDV, LDQ, NCYCLE, INFO REAL (KIND=nag_wp) TOLA, TOLB, ALPHA(N), BETA(N) COMPLEX (KIND=nag_wp) A(LDA,*), B(LDB,*), U(LDU,*), V(LDV,*), Q(LDQ,*), WORK(2*N) CHARACTER(1) JOBU, JOBV, JOBQ
The routine may be called by its LAPACK name ztgsja.

## 3  Description

F08YSF (ZTGSJA) computes the GSVD of the matrices $A$ and $B$ which are assumed to have the form as returned by F08VSF (ZGGSVP)
 $A= n-k-lklk(0A12A13) l 0 0 A23 m-k-l 0 0 0 , if ​ m-k-l ≥ 0; n-k-lklk(0A12A13) m-k 0 0 A23 , if ​ m-k-l < 0 ; B= n-k-lkll(00B13) p-l 0 0 0 ,$
where the $k$ by $k$ matrix ${A}_{12}$ and the $l$ by $l$ matrix ${B}_{13}$ are nonsingular upper triangular, ${A}_{23}$ is $l$ by $l$ upper triangular if $m-k-l\ge 0$ and is $\left(m-k\right)$ by $l$ upper trapezoidal otherwise.
F08YSF (ZTGSJA) computes unitary matrices $Q$, $U$ and $V$, diagonal matrices ${D}_{1}$ and ${D}_{2}$, and an upper triangular matrix $R$ such that
 $UHAQ = D1 0 R , VHBQ = D2 0 R .$
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 $\left(m-k-l\right)\ge 0$ then ${D}_{1}$, ${D}_{2}$ and $R$ have the form
 $D1= klk(I0) l 0 C m-k-l 0 0 ,$
 $D2= kll(0S) p-l 0 0 ,$
 $R = klk(R11R12) l 0 R22 ,$
where $C=\mathrm{diag}\left({\alpha }_{k+1},,,\dots ,,,{\alpha }_{k+l}\right)\text{, }S=\mathrm{diag}\left({\beta }_{k+1},,,\dots ,,,{\beta }_{k+l}\right)$.
If $\left(m-k-l\right)<0$ then ${D}_{1}$, ${D}_{2}$ and $R$ have the form
 $D1= km-kk+l-mk(I00) m-k 0 C 0 ,$
 $D2= km-kk+l-mm-k(0S0) k+l-m 0 0 I p-l 0 0 0 ,$
 $R = km-kk+l-mk(R11R12R13) m-k 0 R22 R23 k+l-m 0 0 R33 ,$
where $C=\mathrm{diag}\left({\alpha }_{k+1},,,\dots ,,,{\alpha }_{m}\right)\text{, }S=\mathrm{diag}\left({\beta }_{k+1},,,\dots ,,,{\beta }_{m}\right)$.
In both cases the diagonal matrix $C$ has real non-negative diagonal elements, the diagonal matrix $S$ has real 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.

## 4  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 http://www.netlib.org/lapack/lug
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

## 5  Parameters

1:     JOBU – CHARACTER(1)Input
On entry: if ${\mathbf{JOBU}}=\text{'U'}$, U must contain a unitary matrix ${U}_{1}$ on entry, and the product ${U}_{1}U$ is returned.
If ${\mathbf{JOBU}}=\text{'I'}$, U is initialized to the unit matrix, and the unitary matrix $U$ is returned.
If ${\mathbf{JOBU}}=\text{'N'}$, $U$ is not computed.
Constraint: ${\mathbf{JOBU}}=\text{'I'}$, $\text{'U'}$ or $\text{'N'}$.
2:     JOBV – CHARACTER(1)Input
On entry: if ${\mathbf{JOBV}}=\text{'V'}$, V must contain a unitary matrix ${V}_{1}$ on entry, and the product ${V}_{1}V$ is returned.
If ${\mathbf{JOBV}}=\text{'I'}$, V is initialized to the unit matrix, and the unitary matrix $V$ is returned.
If ${\mathbf{JOBV}}=\text{'N'}$, $V$ is not computed.
Constraint: ${\mathbf{JOBV}}=\text{'V'}$, $\text{'I'}$ or $\text{'N'}$.
3:     JOBQ – CHARACTER(1)Input
On entry: if ${\mathbf{JOBQ}}=\text{'Q'}$, Q must contain a unitary matrix ${Q}_{1}$ on entry, and the product ${Q}_{1}Q$ is returned.
If ${\mathbf{JOBQ}}=\text{'I'}$, Q is initialized to the unit matrix, and the unitary matrix $Q$ is returned.
If ${\mathbf{JOBQ}}=\text{'N'}$, $Q$ is not computed.
Constraint: ${\mathbf{JOBQ}}=\text{'Q'}$, $\text{'I'}$ or $\text{'N'}$.
4:     M – INTEGERInput
On entry: $m$, the number of rows of the matrix $A$.
Constraint: ${\mathbf{M}}\ge 0$.
5:     P – INTEGERInput
On entry: $p$, the number of rows of the matrix $B$.
Constraint: ${\mathbf{P}}\ge 0$.
6:     N – INTEGERInput
On entry: $n$, the number of columns of the matrices $A$ and $B$.
Constraint: ${\mathbf{N}}\ge 0$.
7:     K – INTEGERInput
8:     L – INTEGERInput
On entry: K and L specify the sizes, $k$ and $l$, of the subblocks of $A$ and $B$, whose GSVD is to be computed by F08YSF (ZTGSJA).
9:     A(LDA,$*$) – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array A must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: the $m$ by $n$ matrix $A$.
On exit: if $m-k-l\ge 0$, ${\mathbf{A}}\left(1:k+l,n-k-l+1:n\right)$ contains the $\left(k+l\right)$ by $\left(k+l\right)$ upper triangular matrix $R$.
If $m-k-l<0$, ${\mathbf{A}}\left(1:m,n-k-l+1:n\right)$ contains the first $m$ rows of the $\left(k+l\right)$ by $\left(k+l\right)$ upper triangular matrix $R$, and the submatrix ${R}_{33}$ is returned in ${\mathbf{B}}\left(m-k+1:l,n+m-k-l+1:n\right)$.
10:   LDA – INTEGERInput
On entry: the first dimension of the array A as declared in the (sub)program from which F08YSF (ZTGSJA) is called.
Constraint: ${\mathbf{LDA}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{M}}\right)$.
11:   B(LDB,$*$) – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array B must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: the $p$ by $n$ matrix $B$.
On exit: if $m-k-l<0$, ${\mathbf{B}}\left(m-k+1:l,n+m-k-l+1:n\right)$ contains the submatrix ${R}_{33}$ of $R$.
12:   LDB – INTEGERInput
On entry: the first dimension of the array B as declared in the (sub)program from which F08YSF (ZTGSJA) is called.
Constraint: ${\mathbf{LDB}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{P}}\right)$.
13:   TOLA – REAL (KIND=nag_wp)Input
14:   TOLB – REAL (KIND=nag_wp)Input
On entry: 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 F08VSF (ZGGSVP), say
 $TOLA=maxM,NAε, TOLB=maxP,NBε,$
where $\epsilon$ is the machine precision.
15:   ALPHA(N) – REAL (KIND=nag_wp) arrayOutput
On exit: see the description of BETA.
16:   BETA(N) – REAL (KIND=nag_wp) arrayOutput
On exit: ALPHA and BETA contain the generalized singular value pairs of $A$ and $B$;
• ${\mathbf{ALPHA}}\left(\mathit{i}\right)=1$, ${\mathbf{BETA}}\left(\mathit{i}\right)=0$, for $\mathit{i}=1,2,\dots ,k$, and
• if $m-k-l\ge 0$, ${\mathbf{ALPHA}}\left(\mathit{i}\right)={\alpha }_{\mathit{i}}$, ${\mathbf{BETA}}\left(\mathit{i}\right)={\beta }_{\mathit{i}}$, for $\mathit{i}=k+1,\dots ,k+l$, or
• if $m-k-l<0$, ${\mathbf{ALPHA}}\left(\mathit{i}\right)={\alpha }_{\mathit{i}}$, ${\mathbf{BETA}}\left(\mathit{i}\right)={\beta }_{\mathit{i}}$, for $\mathit{i}=k+1,\dots ,m$ and ${\mathbf{ALPHA}}\left(\mathit{i}\right)=0$, ${\mathbf{BETA}}\left(\mathit{i}\right)=1$, for $\mathit{i}=m+1,\dots ,k+l$.
Furthermore, if $k+l, ${\mathbf{ALPHA}}\left(\mathit{i}\right)={\mathbf{BETA}}\left(\mathit{i}\right)=0$, for $\mathit{i}=k+l+1,\dots ,n$.
17:   U(LDU,$*$) – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array U must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{M}}\right)$ if ${\mathbf{JOBU}}=\text{'U'}$ or $\text{'I'}$, and at least $1$ otherwise.
On entry: if ${\mathbf{JOBU}}=\text{'U'}$, U must contain an $m$ by $m$ matrix ${U}_{1}$ (usually the unitary matrix returned by F08VSF (ZGGSVP)).
On exit: if ${\mathbf{JOBU}}=\text{'I'}$, U contains the unitary matrix $U$.
If ${\mathbf{JOBU}}=\text{'U'}$, U contains the product ${U}_{1}U$.
If ${\mathbf{JOBU}}=\text{'N'}$, U is not referenced.
18:   LDU – INTEGERInput
On entry: the first dimension of the array U as declared in the (sub)program from which F08YSF (ZTGSJA) is called.
Constraints:
• if ${\mathbf{JOBU}}\ne \text{'N'}$, ${\mathbf{LDU}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{M}}\right)$;
• otherwise ${\mathbf{LDU}}\ge 1$.
19:   V(LDV,$*$) – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array V must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{P}}\right)$ if ${\mathbf{JOBV}}=\text{'V'}$ or $\text{'I'}$, and at least $1$ otherwise.
On entry: if ${\mathbf{JOBV}}=\text{'V'}$, V must contain an $p$ by $p$ matrix ${V}_{1}$ (usually the unitary matrix returned by F08VSF (ZGGSVP)).
On exit: if ${\mathbf{JOBV}}=\text{'I'}$, V contains the unitary matrix $V$.
If ${\mathbf{JOBV}}=\text{'V'}$, V contains the product ${V}_{1}V$.
If ${\mathbf{JOBV}}=\text{'N'}$, V is not referenced.
20:   LDV – INTEGERInput
On entry: the first dimension of the array V as declared in the (sub)program from which F08YSF (ZTGSJA) is called.
Constraints:
• if ${\mathbf{JOBV}}\ne \text{'N'}$, ${\mathbf{LDV}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{P}}\right)$;
• otherwise ${\mathbf{LDV}}\ge 1$.
21:   Q(LDQ,$*$) – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array Q must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$ if ${\mathbf{JOBQ}}=\text{'Q'}$ or $\text{'I'}$, and at least $1$ otherwise.
On entry: if ${\mathbf{JOBQ}}=\text{'Q'}$, Q must contain an $n$ by $n$ matrix ${Q}_{1}$ (usually the unitary matrix returned by F08VSF (ZGGSVP)).
On exit: if ${\mathbf{JOBQ}}=\text{'I'}$, Q contains the unitary matrix $Q$.
If ${\mathbf{JOBQ}}=\text{'Q'}$, Q contains the product ${Q}_{1}Q$.
If ${\mathbf{JOBQ}}=\text{'N'}$, Q is not referenced.
22:   LDQ – INTEGERInput
On entry: the first dimension of the array Q as declared in the (sub)program from which F08YSF (ZTGSJA) is called.
Constraints:
• if ${\mathbf{JOBQ}}\ne \text{'N'}$, ${\mathbf{LDQ}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$;
• otherwise ${\mathbf{LDQ}}\ge 1$.
23:   WORK($2×{\mathbf{N}}$) – COMPLEX (KIND=nag_wp) arrayWorkspace
24:   NCYCLE – INTEGEROutput
On exit: the number of cycles required for convergence.
25:   INFO – INTEGEROutput
On exit: ${\mathbf{INFO}}=0$ unless the routine detects an error (see Section 6).

## 6  Error Indicators and Warnings

Errors or warnings detected by the routine:
${\mathbf{INFO}}<0$
If ${\mathbf{INFO}}=-i$, argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.
${\mathbf{INFO}}=1$
The procedure does not converge after $40$ cycles.

## 7  Accuracy

The computed generalized singular value decomposition is nearly the exact generalized singular value decomposition for nearby matrices $\left(A+E\right)$ and $\left(B+F\right)$, where
 $E2 = O⁡ε A2 and F2= O⁡ε B2 ,$
and $\epsilon$ is the machine precision. See Section 4.12 of Anderson et al. (1999) for further details.

The real analogue of this routine is F08YEF (DTGSJA).

## 9  Example

This example finds the generalized singular value decomposition
 $A = UΣ1 0 R QH , B= VΣ2 0 R QH ,$
of the matrix pair $\left(A,B\right)$, where
 $A = 0.96-0.81i -0.03+0.96i -0.91+2.06i -0.05+0.41i -0.98+1.98i -1.20+0.19i -0.66+0.42i -0.81+0.56i 0.62-0.46i 1.01+0.02i 0.63-0.17i -1.11+0.60i 0.37+0.38i 0.19-0.54i -0.98-0.36i 0.22-0.20i 0.83+0.51i 0.20+0.01i -0.17-0.46i 1.47+1.59i 1.08-0.28i 0.20-0.12i -0.07+1.23i 0.26+0.26i$
and
 $B = 1 0 -1 0 0 1 0 -1 .$

### 9.1  Program Text

Program Text (f08ysfe.f90)

### 9.2  Program Data

Program Data (f08ysfe.d)

### 9.3  Program Results

Program Results (f08ysfe.r)