LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zggsvd()

subroutine zggsvd ( character jobu,
character jobv,
character jobq,
integer m,
integer n,
integer p,
integer k,
integer l,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) alpha,
double precision, dimension( * ) beta,
complex*16, dimension( ldu, * ) u,
integer ldu,
complex*16, dimension( ldv, * ) v,
integer ldv,
complex*16, dimension( ldq, * ) q,
integer ldq,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer info )

ZGGSVD computes the singular value decomposition (SVD) for OTHER matrices

Download ZGGSVD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> This routine is deprecated and has been replaced by routine ZGGSVD3.
!>
!> ZGGSVD computes the generalized singular value decomposition (GSVD)
!> of an M-by-N complex matrix A and P-by-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)-by-(K+L) nonsingular upper
!> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) 
!> 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+1), ... , ALPHA(K+L) ),
!>   S = diag( BETA(K+1),  ... , BETA(K+L) ),
!>   C**2 + S**2 = I.
!>
!>   R is stored in A(1:K+L,N-K-L+1:N) 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+1), ... , ALPHA(M) ),
!>   S = diag( BETA(K+1),  ... , BETA(M) ),
!>   C**2 + S**2 = I.
!>
!>   (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
!>   ( 0  R22 R23 )
!>   in B(M-K+1:L,N+M-K-L+1:N) 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-by-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 orthogonal and X is nonsingular, and 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
!>          JOBU is CHARACTER*1
!>          = 'U':  Unitary matrix U is computed;
!>          = 'N':  U is not computed.
!>
[in]JOBV
!>          JOBV is CHARACTER*1
!>          = 'V':  Unitary matrix V is computed;
!>          = 'N':  V is not computed.
!>
[in]JOBQ
!>          JOBQ is CHARACTER*1
!>          = 'Q':  Unitary matrix Q is computed;
!>          = 'N':  Q is not computed.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrices A and B.  N >= 0.
!>
[in]P
!>          P is INTEGER
!>          The number of rows of the matrix B.  P >= 0.
!>
[out]K
!>          K is INTEGER
!>
[out]L
!>          L is INTEGER
!>
!>          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,out]A
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit, A contains the triangular matrix R, or part of R.
!>          See Purpose for details.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,M).
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,N)
!>          On entry, the P-by-N matrix B.
!>          On exit, B contains part of the triangular matrix R if
!>          M-K-L < 0.  See Purpose for details.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,P).
!>
[out]ALPHA
!>          ALPHA is DOUBLE PRECISION array, dimension (N)
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>
!>          On exit, ALPHA and BETA 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) = C,
!>            BETA(K+1:K+L)  = 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
!>          and
!>            ALPHA(K+L+1:N) = 0
!>            BETA(K+L+1:N)  = 0
!>
[out]U
!>          U is COMPLEX*16 array, dimension (LDU,M)
!>          If JOBU = 'U', U contains the M-by-M unitary matrix U.
!>          If JOBU = 'N', U is not referenced.
!>
[in]LDU
!>          LDU is INTEGER
!>          The leading dimension of the array U. LDU >= max(1,M) if
!>          JOBU = 'U'; LDU >= 1 otherwise.
!>
[out]V
!>          V is COMPLEX*16 array, dimension (LDV,P)
!>          If JOBV = 'V', V contains the P-by-P unitary matrix V.
!>          If JOBV = 'N', V is not referenced.
!>
[in]LDV
!>          LDV is INTEGER
!>          The leading dimension of the array V. LDV >= max(1,P) if
!>          JOBV = 'V'; LDV >= 1 otherwise.
!>
[out]Q
!>          Q is COMPLEX*16 array, dimension (LDQ,N)
!>          If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q.
!>          If JOBQ = 'N', Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q. LDQ >= max(1,N) if
!>          JOBQ = 'Q'; LDQ >= 1 otherwise.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (max(3*N,M,P)+N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (2*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>          On exit, IWORK stores the sorting information. More
!>          precisely, the following loop will sort ALPHA
!>             for I = K+1, min(M,K+L)
!>                 swap ALPHA(I) and ALPHA(IWORK(I))
!>             endfor
!>          such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if INFO = 1, the Jacobi-type procedure failed to
!>                converge.  For further details, see subroutine ZTGSJA.
!>
Internal Parameters:
!>  TOLA    DOUBLE PRECISION
!>  TOLB    DOUBLE PRECISION
!>          TOLA and TOLB are the thresholds to determine the effective
!>          rank of (A**H,B**H)**H. Generally, they are set to
!>                   TOLA = MAX(M,N)*norm(A)*MAZHEPS,
!>                   TOLB = MAX(P,N)*norm(B)*MAZHEPS.
!>          The size of TOLA and TOLB may affect the size of backward
!>          errors of the decomposition.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 332 of file zggsvd.f.

335*
336* -- LAPACK driver routine --
337* -- LAPACK is a software package provided by Univ. of Tennessee, --
338* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
339*
340* .. Scalar Arguments ..
341 CHARACTER JOBQ, JOBU, JOBV
342 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
343* ..
344* .. Array Arguments ..
345 INTEGER IWORK( * )
346 DOUBLE PRECISION ALPHA( * ), BETA( * ), RWORK( * )
347 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
348 $ U( LDU, * ), V( LDV, * ), WORK( * )
349* ..
350*
351* =====================================================================
352*
353* .. Local Scalars ..
354 LOGICAL WANTQ, WANTU, WANTV
355 INTEGER I, IBND, ISUB, J, NCYCLE
356 DOUBLE PRECISION ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
357* ..
358* .. External Functions ..
359 LOGICAL LSAME
360 DOUBLE PRECISION DLAMCH, ZLANGE
361 EXTERNAL lsame, dlamch, zlange
362* ..
363* .. External Subroutines ..
364 EXTERNAL dcopy, xerbla, zggsvp, ztgsja
365* ..
366* .. Intrinsic Functions ..
367 INTRINSIC max, min
368* ..
369* .. Executable Statements ..
370*
371* Decode and test the input parameters
372*
373 wantu = lsame( jobu, 'U' )
374 wantv = lsame( jobv, 'V' )
375 wantq = lsame( jobq, 'Q' )
376*
377 info = 0
378 IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
379 info = -1
380 ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
381 info = -2
382 ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
383 info = -3
384 ELSE IF( m.LT.0 ) THEN
385 info = -4
386 ELSE IF( n.LT.0 ) THEN
387 info = -5
388 ELSE IF( p.LT.0 ) THEN
389 info = -6
390 ELSE IF( lda.LT.max( 1, m ) ) THEN
391 info = -10
392 ELSE IF( ldb.LT.max( 1, p ) ) THEN
393 info = -12
394 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
395 info = -16
396 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
397 info = -18
398 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
399 info = -20
400 END IF
401 IF( info.NE.0 ) THEN
402 CALL xerbla( 'ZGGSVD', -info )
403 RETURN
404 END IF
405*
406* Compute the Frobenius norm of matrices A and B
407*
408 anorm = zlange( '1', m, n, a, lda, rwork )
409 bnorm = zlange( '1', p, n, b, ldb, rwork )
410*
411* Get machine precision and set up threshold for determining
412* the effective numerical rank of the matrices A and B.
413*
414 ulp = dlamch( 'Precision' )
415 unfl = dlamch( 'Safe Minimum' )
416 tola = max( m, n )*max( anorm, unfl )*ulp
417 tolb = max( p, n )*max( bnorm, unfl )*ulp
418*
419 CALL zggsvp( jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola,
420 $ tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, rwork,
421 $ work, work( n+1 ), info )
422*
423* Compute the GSVD of two upper "triangular" matrices
424*
425 CALL ztgsja( jobu, jobv, jobq, m, p, n, k, l, a, lda, b, ldb,
426 $ tola, tolb, alpha, beta, u, ldu, v, ldv, q, ldq,
427 $ work, ncycle, info )
428*
429* Sort the singular values and store the pivot indices in IWORK
430* Copy ALPHA to RWORK, then sort ALPHA in RWORK
431*
432 CALL dcopy( n, alpha, 1, rwork, 1 )
433 ibnd = min( l, m-k )
434 DO 20 i = 1, ibnd
435*
436* Scan for largest ALPHA(K+I)
437*
438 isub = i
439 smax = rwork( k+i )
440 DO 10 j = i + 1, ibnd
441 temp = rwork( k+j )
442 IF( temp.GT.smax ) THEN
443 isub = j
444 smax = temp
445 END IF
446 10 CONTINUE
447 IF( isub.NE.i ) THEN
448 rwork( k+isub ) = rwork( k+i )
449 rwork( k+i ) = smax
450 iwork( k+i ) = k + isub
451 ELSE
452 iwork( k+i ) = k + i
453 END IF
454 20 CONTINUE
455*
456 RETURN
457*
458* End of ZGGSVD
459*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlange
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:113
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ztgsja
subroutine 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)
ZTGSJA
Definition ztgsja.f:377
zggsvp
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)
ZGGSVP
Definition zggsvp.f:263
Here is the call graph for this function: