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

◆ cggsvd3()

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

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

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

Purpose:
!>
!> CGGSVD3 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 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 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 REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL 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 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 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 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 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= 1.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL 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 CTGSJA.
!>
Internal Parameters:
!>  TOLA    REAL
!>  TOLB    REAL
!>          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)*MACHEPS,
!>                   TOLB = MAX(P,N)*norm(B)*MACHEPS.
!>          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
Further Details:
CGGSVD3 replaces the deprecated subroutine CGGSVD.

Definition at line 349 of file cggsvd3.f.

352*
353* -- LAPACK driver routine --
354* -- LAPACK is a software package provided by Univ. of Tennessee, --
355* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
356*
357* .. Scalar Arguments ..
358 CHARACTER JOBQ, JOBU, JOBV
359 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P,
360 $ LWORK
361* ..
362* .. Array Arguments ..
363 INTEGER IWORK( * )
364 REAL ALPHA( * ), BETA( * ), RWORK( * )
365 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
366 $ U( LDU, * ), V( LDV, * ), WORK( * )
367* ..
368*
369* =====================================================================
370*
371* .. Local Scalars ..
372 LOGICAL WANTQ, WANTU, WANTV, LQUERY
373 INTEGER I, IBND, ISUB, J, NCYCLE, LWKOPT
374 REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
375* ..
376* .. External Functions ..
377 LOGICAL LSAME
378 REAL CLANGE, SLAMCH
379 EXTERNAL lsame, clange, slamch
380* ..
381* .. External Subroutines ..
382 EXTERNAL cggsvp3, ctgsja, scopy, xerbla
383* ..
384* .. Intrinsic Functions ..
385 INTRINSIC max, min
386* ..
387* .. Executable Statements ..
388*
389* Decode and test the input parameters
390*
391 wantu = lsame( jobu, 'U' )
392 wantv = lsame( jobv, 'V' )
393 wantq = lsame( jobq, 'Q' )
394 lquery = ( lwork.EQ.-1 )
395 lwkopt = 1
396*
397* Test the input arguments
398*
399 info = 0
400 IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
401 info = -1
402 ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
403 info = -2
404 ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
405 info = -3
406 ELSE IF( m.LT.0 ) THEN
407 info = -4
408 ELSE IF( n.LT.0 ) THEN
409 info = -5
410 ELSE IF( p.LT.0 ) THEN
411 info = -6
412 ELSE IF( lda.LT.max( 1, m ) ) THEN
413 info = -10
414 ELSE IF( ldb.LT.max( 1, p ) ) THEN
415 info = -12
416 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
417 info = -16
418 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
419 info = -18
420 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
421 info = -20
422 ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
423 info = -24
424 END IF
425*
426* Compute workspace
427*
428 IF( info.EQ.0 ) THEN
429 CALL cggsvp3( jobu, jobv, jobq, m, p, n, a, lda, b, ldb,
430 $ tola,
431 $ tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, rwork,
432 $ work, work, -1, info )
433 lwkopt = n + int( work( 1 ) )
434 lwkopt = max( 2*n, lwkopt )
435 lwkopt = max( 1, lwkopt )
436 work( 1 ) = cmplx( lwkopt )
437 END IF
438*
439 IF( info.NE.0 ) THEN
440 CALL xerbla( 'CGGSVD3', -info )
441 RETURN
442 END IF
443 IF( lquery ) THEN
444 RETURN
445 ENDIF
446*
447* Compute the Frobenius norm of matrices A and B
448*
449 anorm = clange( '1', m, n, a, lda, rwork )
450 bnorm = clange( '1', p, n, b, ldb, rwork )
451*
452* Get machine precision and set up threshold for determining
453* the effective numerical rank of the matrices A and B.
454*
455 ulp = slamch( 'Precision' )
456 unfl = slamch( 'Safe Minimum' )
457 tola = real( max( m, n ) )*max( anorm, unfl )*ulp
458 tolb = real( max( p, n ) )*max( bnorm, unfl )*ulp
459*
460 CALL cggsvp3( jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola,
461 $ tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, rwork,
462 $ work, work( n+1 ), lwork-n, info )
463*
464* Compute the GSVD of two upper "triangular" matrices
465*
466 CALL ctgsja( jobu, jobv, jobq, m, p, n, k, l, a, lda, b, ldb,
467 $ tola, tolb, alpha, beta, u, ldu, v, ldv, q, ldq,
468 $ work, ncycle, info )
469*
470* Sort the singular values and store the pivot indices in IWORK
471* Copy ALPHA to RWORK, then sort ALPHA in RWORK
472*
473 CALL scopy( n, alpha, 1, rwork, 1 )
474 ibnd = min( l, m-k )
475 DO 20 i = 1, ibnd
476*
477* Scan for largest ALPHA(K+I)
478*
479 isub = i
480 smax = rwork( k+i )
481 DO 10 j = i + 1, ibnd
482 temp = rwork( k+j )
483 IF( temp.GT.smax ) THEN
484 isub = j
485 smax = temp
486 END IF
487 10 CONTINUE
488 IF( isub.NE.i ) THEN
489 rwork( k+isub ) = rwork( k+i )
490 rwork( k+i ) = smax
491 iwork( k+i ) = k + isub
492 ELSE
493 iwork( k+i ) = k + i
494 END IF
495 20 CONTINUE
496*
497 work( 1 ) = cmplx( lwkopt )
498 RETURN
499*
500* End of CGGSVD3
501*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
cggsvp3
subroutine cggsvp3(jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola, tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, rwork, tau, work, lwork, info)
CGGSVP3
Definition cggsvp3.f:276
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:113
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ctgsja
subroutine 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)
CTGSJA
Definition ctgsja.f:377
Here is the call graph for this function:
Here is the caller graph for this function: