Contents


NAME

     zggsvp - compute unitary matrices U, V and Q such that    N-
     K-L K L  U'*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0

SYNOPSIS

     SUBROUTINE ZGGSVP(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
           TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, TAU, WORK,
           INFO)

     CHARACTER * 1 JOBU, JOBV, JOBQ
     DOUBLE  COMPLEX  A(LDA,*),  B(LDB,*),  U(LDU,*),   V(LDV,*),
     Q(LDQ,*), TAU(*), WORK(*)
     INTEGER M, P, N, LDA, LDB, K, L, LDU, LDV, LDQ, INFO
     INTEGER IWORK(*)
     DOUBLE PRECISION TOLA, TOLB
     DOUBLE PRECISION RWORK(*)

     SUBROUTINE ZGGSVP_64(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
           TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, TAU, WORK,
           INFO)

     CHARACTER * 1 JOBU, JOBV, JOBQ
     DOUBLE  COMPLEX  A(LDA,*),  B(LDB,*),  U(LDU,*),   V(LDV,*),
     Q(LDQ,*), TAU(*), WORK(*)
     INTEGER*8 M, P, N, LDA, LDB, K, L, LDU, LDV, LDQ, INFO
     INTEGER*8 IWORK(*)
     DOUBLE PRECISION TOLA, TOLB
     DOUBLE PRECISION RWORK(*)

  F95 INTERFACE
     SUBROUTINE GGSVP(JOBU, JOBV, JOBQ, [M], [P], [N], A, [LDA], B, [LDB],
            TOLA, TOLB, K, L, U, [LDU], V, [LDV], Q, [LDQ], [IWORK], [RWORK],
            [TAU], [WORK], [INFO])

     CHARACTER(LEN=1) :: JOBU, JOBV, JOBQ
     COMPLEX(8), DIMENSION(:) :: TAU, WORK
     COMPLEX(8), DIMENSION(:,:) :: A, B, U, V, Q
     INTEGER :: M, P, N, LDA, LDB, K, L, LDU, LDV, LDQ, INFO
     INTEGER, DIMENSION(:) :: IWORK
     REAL(8) :: TOLA, TOLB
     REAL(8), DIMENSION(:) :: RWORK

     SUBROUTINE GGSVP_64(JOBU, JOBV, JOBQ, [M], [P], [N], A, [LDA], B,
            [LDB], TOLA, TOLB, K, L, U, [LDU], V, [LDV], Q, [LDQ], [IWORK],
            [RWORK], [TAU], [WORK], [INFO])

     CHARACTER(LEN=1) :: JOBU, JOBV, JOBQ
     COMPLEX(8), DIMENSION(:) :: TAU, WORK
     COMPLEX(8), DIMENSION(:,:) :: A, B, U, V, Q
     INTEGER(8) :: M, P, N, LDA, LDB, K, L, LDU, LDV, LDQ, INFO
     INTEGER(8), DIMENSION(:) :: IWORK
     REAL(8) :: TOLA, TOLB
     REAL(8), DIMENSION(:) :: RWORK

  C INTERFACE
     #include <sunperf.h>

     void zggsvp(char jobu, char jobv, char jobq, int m,  int  p,
               int  n,  doublecomplex  *a, int lda, doublecomplex
               *b, int ldb, double tola, double tolb, int *k, int
               *l,  doublecomplex  *u, int ldu, doublecomplex *v,
               int ldv, doublecomplex *q, int ldq, int *info);

     void zggsvp_64(char jobu, char jobv, char jobq, long m, long
               p,  long n, doublecomplex *a, long lda, doublecom-
               plex *b, long ldb, double tola, double tolb,  long
               *k,  long  *l,  doublecomplex  *u, long ldu, doub-
               lecomplex *v, long  ldv,  doublecomplex  *q,  long
               ldq, long *info);

PURPOSE

     zggsvp computes unitary matrices U, V and Q such that
                   L ( 0     0   A23 )
               M-K-L ( 0     0    0  )

                      N-K-L  K    L
             =     K ( 0    A12  A13 )  if M-K-L < 0;
                 M-K ( 0     0   A23 )

                    N-K-L  K    L
      V'*B*Q =   L ( 0     0   B13 )
               P-L ( 0     0    0  )

     where the K-by-K matrix A12 and L-by-L matrix B13  are  non-
     singular upper triangular; A23 is L-by-L upper triangular if
     M-K-L >= 0, otherwise A23 is (M-K)-by-L  upper  trapezoidal.
     K+L  = the effective numerical rank of the (M+P)-by-N matrix
     (A',B')'.  Z' denotes the conjugate transpose of Z.

     This decomposition is the preprocessing step  for  computing
     the  Generalized  Singular  Value  Decomposition (GSVD), see
     subroutine ZGGSVD.

ARGUMENTS

     JOBU (input)
               = 'U':  Unitary matrix U is computed;
               = 'N':  U is not computed.

     JOBV (input)
               = 'V':  Unitary matrix V is computed;
               = 'N':  V is not computed.

     JOBQ (input)
               = 'Q':  Unitary matrix Q is computed;
               = 'N':  Q is not computed.

     M (input) The number of rows of the matrix A.  M >= 0.

     P (input) The number of rows of the matrix B.  P >= 0.

     N (input) The number of columns of the matrices A and B.   N
               >= 0.

     A (input/output)
               On entry, the M-by-N matrix A.  On  exit,  A  con-
               tains   the  triangular  (or  trapezoidal)  matrix
               described in the Purpose section.

     LDA (input)
               The leading dimension  of  the  array  A.  LDA  >=
               max(1,M).

     B (input/output)
               On entry, the P-by-N matrix B.  On  exit,  B  con-
               tains  the triangular matrix described in the Pur-
               pose section.

     LDB (input)
               The leading dimension  of  the  array  B.  LDB  >=
               max(1,P).

     TOLA (input)
               TOLA and TOLB are the thresholds to determine  the
               effective  numerical  rank  of matrix B and a sub-
               block of A. Generally, they  are  set  to  TOLA  =
               MAX(M,N)*norm(A)*MACHEPS,          TOLB          =
               MAX(P,N)*norm(B)*MACHEPS.  The size  of  TOLA  and
               TOLB may affect the size of backward errors of the
               decomposition.

     TOLB (input)
               See description of TOLA.

     K (output)
               On exit, K and L specify the dimension of the sub-
               blocks  described  in  Purpose  section.   K + L =
               effective numerical rank of (A',B')'.

     L (output)
               See the description of K.

     U (input) If JOBU = 'U', U contains the  unitary  matrix  U.
               If JOBU = 'N', U is not referenced.

     LDU (input)
               The leading dimension  of  the  array  U.  LDU  >=
               max(1,M) if JOBU = 'U'; LDU >= 1 otherwise.

     V (input) If JOBV = 'V', V contains the  unitary  matrix  V.
               If JOBV = 'N', V is not referenced.

     LDV (input)
               The leading dimension  of  the  array  V.  LDV  >=
               max(1,P) if JOBV = 'V'; LDV >= 1 otherwise.

     Q (input) If JOBQ = 'Q', Q contains the  unitary  matrix  Q.
               If JOBQ = 'N', Q is not referenced.

     LDQ (input)
               The leading dimension  of  the  array  Q.  LDQ  >=
               max(1,N) if JOBQ = 'Q'; LDQ >= 1 otherwise.

     IWORK (workspace)
               dimension(N)

     RWORK (workspace)
               dimension(2*N)

     TAU (workspace)
               dimension(N)
     WORK (workspace)
               dimension(MAX(3*N,M,P))

     INFO (output)
               = 0:  successful exit
               < 0:  if INFO = -i, the i-th argument had an ille-
               gal value.

FURTHER DETAILS

     The subroutine uses LAPACK subroutine CGEQPF for the QR fac-
     torization  with  column  pivoting  to  detect the effective
     numerical rank of the a matrix. It  may  be  replaced  by  a
     better rank determination strategy.