zgejsv zgesdd zgesvd zgesvdq zgesvdx zggsvd3

zggsvd3() void zggsvd3 ( char  jobu, char  jobv, char  jobq, int  m, int  n, int  p, int *  k, int *  l, int  lda, doublecomplex  a[], int  ldb, doublecomplex  b[], double  alpha[], double  beta[], int  ldu, doublecomplex  u[], int  ldv, doublecomplex  v[], int  ldq, doublecomplex  q[], doublecomplex  work[], int  lwork, double  rwork[], int  iwork[], int *  info  ) Generalized singular value decomposition (GSVD) of complex matrices Purpose This routine computes the generalized singular value decomposition (GSVD) of an m x n complex matrix A and p x n complex matrix B: U^H*A*Q = D1*(0 R), V^H*B*Q = D2*(0 R) where, U, V and Q are unitary matrices. Let k+l = the effective numerical rank of the matrix (A^H, B^H)^H, then R is a k+l x k+l nonsingular upper triangular matrix, D1 and D2 are m x k+l and p x k+l "diagonal" matrices and of the following structures, respectively:    If m-k-l >= 0, k l D1 = k (I 0) l (0 C) m-k-l (0 0) k l D2 = l (0 S) p-l (0 0) n-k-l k l (0 R) = k ( 0 R11 R12) l ( 0 0 R22) where C = diag(alpha[k], ... , alpha[k+l-1]), S = diag(beta[k]), ... , beta[k+l-1]), C^2 + S^2 = I. R is stored in a[n-k-l to n-1][0 to k+l-1] on exit. If m-k-l < 0, k m-k k+l-m D1 = k (I 0 0) m-k (0 C 0) k m-k k+l-m D2 = m-k (0 S 0) k+l-m (0 0 I) p-l (0 0 0) n-k-l k m-k k+l-m (0 R) = k ( 0 R11 R12 R13) m-k ( 0 0 R22 R23) k+l-m ( 0 0 0 R33) where C = diag(alpha[k], ... , alpha[m-1]) S = diag(beta[k], ... , beta[m-1]) C^2 + S^2 = I (R11 R12 R13) is stored in a[n-k-l to n-1][0 to m-1], and R33 is stored ( 0 R22 R23) in b[n+m-k-l to n-1][m-k to l-1] on exit. The routine computes C, S, R, and optionally the unitary transformation matrices U, V and Q. In particular, if B is an n x n nonsingular matrix, then the GSVD of A and B implicitly gives the SVD of A*inv(B). A*inv(B) = U*(D1*inv(D2))*V^H If (A^H, B^H)^H has orthonormal columns, then the GSVD of A and B is also equal to the CS decomposition of A and B. Furthermore, the GSVD can be used to derive the solution of the eigenvalue problem: A^H*A x = lambda B^H*B x In some literature, the GSVD of A and B is presented in the form U^H*A*X = (0 D1), V^H*B*X = (0 D2) where U and V are unitary and X is nonsingular, D1 and D2 are "diagonal". The former GSVD form can be converted to the latter form by taking the nonsingular matrix X as X = Q*(I 0 ) (0 inv(R)) Parameters [in] jobu = 'U': Unitary matrix U is computed. = 'N': U is not computed. [in] jobv = 'V': Unitary matrix V is computed. = 'N': V is not computed. [in] jobq = "Q": Unitary matrix Q is computed. = 'N': Q is not computed. [in] m The number of rows of the matrix A. (m >= 0) [in] n The number of columns of the matrices A and B. (n >= 0) [in] p The number of rows of the matrix B. (p >= 0) [out] k [out] l On exit, k and l specify the dimension of the subblocks described in Purpose. (k + l = effective numerical rank of (A^H, B^H)^H) [in] lda Leading dimension of the two dimensional array a[][]. (lda >= max(1, m)) [in,out] a[][] Array a[la][lda] (la >= n) [in] m x n matrix A [out] a[][] contains the triangular matrix R, or part of R. See Purpose for details. [in] ldb Leading dimension of the two dimensional array b[][]. (ldb >= max(1, p)) [in,out] b[][] Array b[lb][ldb] (lb >= n) [in] p x n matrix B. [out] b[][] contains the triangular matrix R if m-k-l < 0. See Purpose for details. [out] alpha[] Array alpha[lalpha] (lalpha >= n) [out] beta[] Array beta[lbeta] (lbeta >= n) alpha[] and beta[] contain the generalized singular value pairs of A and B;   alpha[0 to k-1] = 1,   beta[0 to k-1] = 0, and if m-k-l >= 0,   alpha[k to k+l-1] = C,   beta[k to k+l-1] = S, or if m-k-l < 0,   alpha[k to m-1] = C, alpha[m to k+l-1] = 0,   beta[k to m-1] = S, beta[m to k+l-1] = 1, and   alpha[k+l to n-1] = 0,   beta[k+l to n-1] = 0. [in] ldu Leading dimension of the two dimensional array u[][]. (ldu >= 1 if jobu = 'N', ldu >= max(1, m) if jobu = 'U') [out] u[][] Array u[lu][ldu] (lu >= m) jobu = 'U': u[][] contains the m x m unitary matrix U. jobu = 'N': u[][] is not referenced. [in] ldv Leading dimension of the two dimensional array v[][]. (ldv >= 1 if jobv = 'N', ldu >= max(1, p) if jobv = 'V') [out] v[][] Array v[lv][ldv] (lv >= p) jobv = 'V': v[][] contains the p x p unitary matrix V. jobv = 'N': v[][] is not referenced. [in] ldq Leading dimension of the two dimensional array q[][]. (ldq >= 1 if jobq = 'N', ldq >= max(1, n) if jobq = 'Q') [out] q[][] Array q[lq][ldq] (lq >= n) jobq = 'Q': q[][] contains the n x n unitary matrix Q. jobq = 'N': q[][] is not referenced. [out] work[] Array work[max(1, lwork)] Work array. If info = 0, work[0] returns the optimal lwork. [in] lwork The dimension of the array work[]. If lwork = -1, then a workspace query is assumed. The routine only calculates the optimal size of the work[] array, and returns the value in work[0]. [out] rwork[] Array rwork[lrwork] (lrwork >= 2*n) Work array. [out] iwork[] Array iwork[liwork] (liwork >= n) Work array. On exit, iwork[] stores the sorting information. More precisely, the following loop will sort alpha[]. for i = k, min(m, k+l)-1 swap alpha[i] and alpha[iwork[i]-1]. end for such that alpha[0] >= alpha[1] >= ... >= alpha[n-1]. [out] info = 0: Successful exit = -1: The argument jobu had an illegal value (jobu != 'U' nor 'N') = -2: The argument jobv had an illegal value (jobv != 'V' nor 'N') = -3: The argument jobq had an illegal value (jobq != 'Q' nor 'N') = -4: The argument m had an illegal value (m < 0) = -5: The argument n had an illegal value (n < 0) = -6: The argument p had an illegal value (p < 0) = -9: The argument lda had an illegal value (lda < max(1, m)) = -11: The argument ldb had an illegal value (ldb < max(1, p)) = -15: The argument ldu had an illegal value (ldu too small) = -17: The argument ldv had an illegal value (ldv too small) = -19: The argument ldq had an illegal value (ldq too small) = -24: The argument lwork had an illegal value (lwork too small) = 1: The Jacobi-type procedure failed to converge Reference LAPACK