# ?tgsja

Computes the generalized SVD of two upper triangular or trapezoidal matrices.

## Syntax

FORTRAN 77:

call stgsja(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)

call dtgsja(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)

call ctgsja(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)

call ztgsja(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)

Fortran 95:

call tgsja(a, b, tola, tolb, k, l [,u] [,v] [,q] [,jobu] [,jobv] [,jobq] [,alpha] [,beta] [,ncycle] [,info])

## Description

This routine computes the generalized singular value decomposition (GSVD) of two real/complex upper triangular (or trapezoidal) matrices A and B. On entry, it is assumed that matrices A and B have the following forms, which may be obtained by the preprocessing subroutine ?ggsvp from a general m-by-n matrix A and p-by-n matrix B:

where the k-by-k matrix A12 and l-by-l matrix B13 are nonsingular upper triangular; A23 is l-by-l upper triangular if `m-k-l ≥0`, otherwise A23 is (m-k)-by-l upper trapezoidal.

On exit,

`UH*A*Q = D1*(0 R)`, `VH*B*Q = D2*(0 R)`,

where U, V and Q are orthogonal/unitary matrices, R is a nonsingular upper triangular matrix, and D1 and D2 are “diagonal” matrices, which are of the following structures:

If `m-k-l ≥0`,

where

`C = diag(alpha(k+1),...,alpha(k+l))`

`S = diag(beta(k+1),...,beta(k+l))`

`C2 + S2 = I`

R is stored in a(1:k+l, n-k-l+1:n ) on exit.

If `m-k-l < 0`,

where

`C = diag(alpha(K+1),...,alpha(m))`,

`S = diag(beta(K+1),...,beta(m))`,

`C2 + S2 = I`

On exit, is stored in a(1:m, n-k-l+1:n ) and R33 is stored

in b(m-k+1:l, n+m-k-l+1:n ).

The computation of the orthogonal/unitary transformation matrices U, V or Q is optional. These matrices may either be formed explicitly, or they may be postmultiplied into input matrices U1, V1, or Q1.

## Input Parameters

jobu

  CHARACTER*1. Must be 'U', 'I', or 'N'.

  If `jobu = 'U'`, u must contain an orthogonal/unitary matrix U1 on entry.

  If `jobu = 'I'`, u is initialized to the unit matrix.

  If `jobu = 'N'`, u is not computed.

jobv

  CHARACTER*1. Must be 'V', 'I', or 'N'.

  If `jobv = 'V'`, v must contain an orthogonal/unitary matrix V1 on entry.

  If `jobv = 'I'`, v is initialized to the unit matrix.

  If `jobv = 'N'`, v is not computed.

jobq

  CHARACTER*1. Must be 'Q', 'I', or 'N'.

  If `jobq = 'Q'`, q must contain an orthogonal/unitary matrix Q1 on entry.

  If `jobq = 'I'`, q is initialized to the unit matrix.

  If `jobq = 'N'`, q is not computed.

m

  INTEGER. The number of rows of the matrix A (m ≥ 0).

p

  INTEGER. The number of rows of the matrix B (p ≥ 0).

n

  INTEGER. The number of columns of the matrices A and B (n ≥ 0).

k, l

  INTEGER. Specify the subblocks in the input matrices A and B, whose GSVD is computed.

a,b,u,v,q,work

  REAL for stgsja

  DOUBLE PRECISION for dtgsja

  COMPLEX for ctgsja

  DOUBLE COMPLEX for ztgsja.

  Arrays:

  a(lda,*) contains the m-by-n matrix A.

  The second dimension of a must be at least max(1, n).

  b(ldb,*) contains the p-by-n matrix B.

  The second dimension of b must be at least max(1, n).

  If `jobu = 'U'`, u(ldu,*) must contain a matrix U1 (usually the orthogonal/unitary matrix returned by ?ggsvp).

  The second dimension of u must be at least max(1, m).

  If `jobv = 'V'`, v(ldv,*) must contain a matrix V1 (usually the orthogonal/unitary matrix returned by ?ggsvp).

  The second dimension of v must be at least max(1, p).

  If `jobq = 'Q'`, q(ldq,*) must contain a matrix Q1 (usually the orthogonal/unitary matrix returned by ?ggsvp).

  The second dimension of q must be at least max(1, n).

  work(*) is a workspace array.

  The dimension of work must be at least max(1, 2n).

lda

  INTEGER. The first dimension of a; at least max(1, m).

ldb

  INTEGER. The first dimension of b; at least max(1, p).

ldu

  INTEGER. The first dimension of the array u .

  `ldu ≥ max(1, m)` if `jobu = 'U'`; `ldu ≥ 1` otherwise.

ldv

  INTEGER. The first dimension of the array v .

  `ldv ≥ max(1, p)` if `jobv = 'V'`; `ldv ≥ 1` otherwise.

ldq

  INTEGER. The first dimension of the array q .

  `ldq ≥ max(1, n)` if `jobq = 'Q'`; `ldq ≥ 1` otherwise.

tola, tolb

  REAL for single-precision flavors

  DOUBLE PRECISION for double-precision flavors.

  tola and tolb are the convergence criteria for the Jacobi-Kogbetliantz iteration procedure. Generally, they are the same as used in ?ggsvp:

  `tola = max(m, n)*|A|*MACHEPS`,

  `tolb = max(p, n)*|B|*MACHEPS`.

## Output Parameters

a

  On exit, a(n-k+1:n, 1:min(k+l, m)) contains the triangular matrix R or part of R.

b

  On exit, if necessary, b(m-k+1: l, n+m-k-l+1: n)) contains a part of R.

alpha, beta

  REAL for single-precision flavors

  DOUBLE PRECISION for double-precision flavors.

  Arrays, DIMENSION at least max(1, n). Contain the generalized singular value pairs of A and B:

  `alpha(1:k) = 1`,

  `beta(1:k) = 0`,

  and if `m-k-l ≥ 0`,

  `alpha(k+1:k+l) = diag(C)`,

  `beta(k+1:k+l) = diag(S)`,

  or if `m-k-l < 0`,

  `alpha(k+1:m)= C`, `alpha(m+1:k+l)=0`

  `beta(K+1:M) = S`,

  `beta(m+1:k+l) = 1`.

  Furthermore, if `k+l < n`,

  `alpha(k+l+1:n)= 0` and

  `beta(k+l+1:n) = 0`.

u

  If `jobu = 'I'`, u contains the orthogonal/unitary matrix U.

  If `jobu = 'U'`, u contains the product U1*U.

  If `jobu = 'N'`, u is not referenced.

v

  If `jobv = 'I'`, v contains the orthogonal/unitary matrix U.

  If `jobv = 'V'`, v contains the product V1*V.

  If `jobv = 'N'`, v is not referenced.

q

  If `jobq = 'I'`, q contains the orthogonal/unitary matrix U.

  If `jobq = 'Q'`, q contains the product Q1*Q.

  If `jobq = 'N'`, q is not referenced.

ncycle

  INTEGER. The number of cycles required for convergence.

info

  INTEGER.

  If `info = 0`, the execution is successful.

  If `info = -i`, the i-th parameter had an illegal value.

  If `info = 1`, the procedure does not converge after MAXIT cycles.

## Fortran 95 Interface Notes

Routines in Fortran 95 interface have fewer arguments in the calling sequence than their FORTRAN 77 counterparts. For general conventions applied to skip redundant or restorable arguments, see Fortran 95 Interface Conventions.

Specific details for the routine tgsja interface are the following:

a

  Holds the matrix A of size (m,n).

b

  Holds the matrix B of size (p,n).

u

  Holds the matrix U of size (m,m).

v

  Holds the matrix V of size (p,p).

q

  Holds the matrix Q of size (n,n).

alpha

  Holds the vector of length (n).

beta

  Holds the vector of length (n).

jobu

  If omitted, this argument is restored based on the presence of argument u as follows:

  `jobu = 'U'`, if u is present,

  `jobu = 'N'`, if u is omitted.

  If present, jobu must be equal to 'I' or 'U' and the argument u must also be present.

  Note that there will be an error condition if jobu is present and u omitted.

jobv

  If omitted, this argument is restored based on the presence of argument v as follows:

  `jobv = 'V'`, if v is present,

  `jobv = 'N'`, if v is omitted.

  If present, jobv must be equal to 'I' or 'V' and the argument v must also be present.

  Note that there will be an error condition if jobv is present and v omitted.

jobq

  If omitted, this argument is restored based on the presence of argument q as follows:

  `jobq = 'Q'`, if q is present,

  `jobq = 'N'`, if q is omitted.

  If present, jobq must be equal to 'I' or 'Q' and the argument q must also be present.

  Note that there will be an error condition if jobq is present and q omitted.

Parent topic: Generalized Singular Value Decomposition