LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
ctgsja.f File Reference

Go to the source code of this file.

Functions/Subroutines

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 More...
 

Function/Subroutine Documentation

subroutine ctgsja ( character  JOBU,
character  JOBV,
character  JOBQ,
integer  M,
integer  P,
integer  N,
integer  K,
integer  L,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
real  TOLA,
real  TOLB,
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  NCYCLE,
integer  INFO 
)

CTGSJA

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

Purpose:
 CTGSJA computes the generalized singular value decomposition (GSVD)
 of two 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 CGGSVP
 from a general M-by-N matrix A and P-by-N matrix B:

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

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

            N-K-L  K    L
    B =  L ( 0     0   B13 )
       P-L ( 0     0    0  )

 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,

        U**H *A*Q = D1*( 0 R ),    V**H *B*Q = D2*( 0 R ),

 where U, V and Q are 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,

                     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 ) K
             L (  0    0   R22 ) L

 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.

 R = ( 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 computation of the 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.
Parameters
[in]JOBU
          JOBU is CHARACTER*1
          = 'U':  U must contain a unitary matrix U1 on entry, and
                  the product U1*U is returned;
          = 'I':  U is initialized to the unit matrix, and the
                  unitary matrix U is returned;
          = 'N':  U is not computed.
[in]JOBV
          JOBV is CHARACTER*1
          = 'V':  V must contain a unitary matrix V1 on entry, and
                  the product V1*V is returned;
          = 'I':  V is initialized to the unit matrix, and the
                  unitary matrix V is returned;
          = 'N':  V is not computed.
[in]JOBQ
          JOBQ is CHARACTER*1
          = 'Q':  Q must contain a unitary matrix Q1 on entry, and
                  the product Q1*Q is returned;
          = 'I':  Q is initialized to the unit matrix, and the
                  unitary matrix Q is returned;
          = 'N':  Q is not computed.
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]P
          P is INTEGER
          The number of rows of the matrix B.  P >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrices A and B.  N >= 0.
[in]K
          K is INTEGER
[in]L
          L is INTEGER

          K and L specify the subblocks in the input matrices A and B:
          A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N)
          of A and B, whose GSVD is going to be computed by CTGSJA.
          See Further Details.
[in,out]A
          A is COMPLEX array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, A(N-K+1:N,1:MIN(K+L,M) ) 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, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
          a part of R.  See Purpose for details.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,P).
[in]TOLA
          TOLA is REAL
[in]TOLB
          TOLB is REAL

          TOLA and TOLB are the convergence criteria for the Jacobi-
          Kogbetliantz iteration procedure. Generally, they are the
          same as used in the preprocessing step, say
              TOLA = MAX(M,N)*norm(A)*MACHEPS,
              TOLB = MAX(P,N)*norm(B)*MACHEPS.
[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) = 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
            BETA(K+L+1:N)  = 0.
[in,out]U
          U is COMPLEX array, dimension (LDU,M)
          On entry, if JOBU = 'U', U must contain a matrix U1 (usually
          the unitary matrix returned by CGGSVP).
          On exit,
          if JOBU = 'I', U contains the unitary matrix U;
          if JOBU = 'U', U contains the product U1*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.
[in,out]V
          V is COMPLEX array, dimension (LDV,P)
          On entry, if JOBV = 'V', V must contain a matrix V1 (usually
          the unitary matrix returned by CGGSVP).
          On exit,
          if JOBV = 'I', V contains the unitary matrix V;
          if JOBV = 'V', V contains the product V1*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.
[in,out]Q
          Q is COMPLEX array, dimension (LDQ,N)
          On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
          the unitary matrix returned by CGGSVP).
          On exit,
          if JOBQ = 'I', Q contains the unitary matrix Q;
          if JOBQ = 'Q', Q contains the product Q1*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 (2*N)
[out]NCYCLE
          NCYCLE is INTEGER
          The number of cycles required for convergence.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          = 1:  the procedure does not converge after MAXIT cycles.
Internal Parameters:
  MAXIT   INTEGER
          MAXIT specifies the total loops that the iterative procedure
          may take. If after MAXIT cycles, the routine fails to
          converge, we return INFO = 1.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  CTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce
  min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
  matrix B13 to the form:

           U1**H *A13*Q1 = C1*R1; V1**H *B13*Q1 = S1*R1,

  where U1, V1 and Q1 are unitary matrix.
  C1 and S1 are diagonal matrices satisfying

                C1**2 + S1**2 = I,

  and R1 is an L-by-L nonsingular upper triangular matrix.

Definition at line 381 of file ctgsja.f.

381 *
382 * -- LAPACK computational routine (version 3.4.0) --
383 * -- LAPACK is a software package provided by Univ. of Tennessee, --
384 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
385 * November 2011
386 *
387 * .. Scalar Arguments ..
388  CHARACTER jobq, jobu, jobv
389  INTEGER info, k, l, lda, ldb, ldq, ldu, ldv, m, n,
390  $ ncycle, p
391  REAL tola, tolb
392 * ..
393 * .. Array Arguments ..
394  REAL alpha( * ), beta( * )
395  COMPLEX a( lda, * ), b( ldb, * ), q( ldq, * ),
396  $ u( ldu, * ), v( ldv, * ), work( * )
397 * ..
398 *
399 * =====================================================================
400 *
401 * .. Parameters ..
402  INTEGER maxit
403  parameter( maxit = 40 )
404  REAL zero, one
405  parameter( zero = 0.0e+0, one = 1.0e+0 )
406  COMPLEX czero, cone
407  parameter( czero = ( 0.0e+0, 0.0e+0 ),
408  $ cone = ( 1.0e+0, 0.0e+0 ) )
409 * ..
410 * .. Local Scalars ..
411 *
412  LOGICAL initq, initu, initv, upper, wantq, wantu, wantv
413  INTEGER i, j, kcycle
414  REAL a1, a3, b1, b3, csq, csu, csv, error, gamma,
415  $ rwk, ssmin
416  COMPLEX a2, b2, snq, snu, snv
417 * ..
418 * .. External Functions ..
419  LOGICAL lsame
420  EXTERNAL lsame
421 * ..
422 * .. External Subroutines ..
423  EXTERNAL ccopy, clags2, clapll, claset, crot, csscal,
424  $ slartg, xerbla
425 * ..
426 * .. Intrinsic Functions ..
427  INTRINSIC abs, conjg, max, min, real
428 * ..
429 * .. Executable Statements ..
430 *
431 * Decode and test the input parameters
432 *
433  initu = lsame( jobu, 'I' )
434  wantu = initu .OR. lsame( jobu, 'U' )
435 *
436  initv = lsame( jobv, 'I' )
437  wantv = initv .OR. lsame( jobv, 'V' )
438 *
439  initq = lsame( jobq, 'I' )
440  wantq = initq .OR. lsame( jobq, 'Q' )
441 *
442  info = 0
443  IF( .NOT.( initu .OR. wantu .OR. lsame( jobu, 'N' ) ) ) THEN
444  info = -1
445  ELSE IF( .NOT.( initv .OR. wantv .OR. lsame( jobv, 'N' ) ) ) THEN
446  info = -2
447  ELSE IF( .NOT.( initq .OR. wantq .OR. lsame( jobq, 'N' ) ) ) THEN
448  info = -3
449  ELSE IF( m.LT.0 ) THEN
450  info = -4
451  ELSE IF( p.LT.0 ) THEN
452  info = -5
453  ELSE IF( n.LT.0 ) THEN
454  info = -6
455  ELSE IF( lda.LT.max( 1, m ) ) THEN
456  info = -10
457  ELSE IF( ldb.LT.max( 1, p ) ) THEN
458  info = -12
459  ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
460  info = -18
461  ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
462  info = -20
463  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
464  info = -22
465  END IF
466  IF( info.NE.0 ) THEN
467  CALL xerbla( 'CTGSJA', -info )
468  RETURN
469  END IF
470 *
471 * Initialize U, V and Q, if necessary
472 *
473  IF( initu )
474  $ CALL claset( 'Full', m, m, czero, cone, u, ldu )
475  IF( initv )
476  $ CALL claset( 'Full', p, p, czero, cone, v, ldv )
477  IF( initq )
478  $ CALL claset( 'Full', n, n, czero, cone, q, ldq )
479 *
480 * Loop until convergence
481 *
482  upper = .false.
483  DO 40 kcycle = 1, maxit
484 *
485  upper = .NOT.upper
486 *
487  DO 20 i = 1, l - 1
488  DO 10 j = i + 1, l
489 *
490  a1 = zero
491  a2 = czero
492  a3 = zero
493  IF( k+i.LE.m )
494  $ a1 = REAL( A( K+I, N-L+I ) )
495  IF( k+j.LE.m )
496  $ a3 = REAL( A( K+J, N-L+J ) )
497 *
498  b1 = REAL( B( I, N-L+I ) )
499  b3 = REAL( B( J, N-L+J ) )
500 *
501  IF( upper ) THEN
502  IF( k+i.LE.m )
503  $ a2 = a( k+i, n-l+j )
504  b2 = b( i, n-l+j )
505  ELSE
506  IF( k+j.LE.m )
507  $ a2 = a( k+j, n-l+i )
508  b2 = b( j, n-l+i )
509  END IF
510 *
511  CALL clags2( upper, a1, a2, a3, b1, b2, b3, csu, snu,
512  $ csv, snv, csq, snq )
513 *
514 * Update (K+I)-th and (K+J)-th rows of matrix A: U**H *A
515 *
516  IF( k+j.LE.m )
517  $ CALL crot( l, a( k+j, n-l+1 ), lda, a( k+i, n-l+1 ),
518  $ lda, csu, conjg( snu ) )
519 *
520 * Update I-th and J-th rows of matrix B: V**H *B
521 *
522  CALL crot( l, b( j, n-l+1 ), ldb, b( i, n-l+1 ), ldb,
523  $ csv, conjg( snv ) )
524 *
525 * Update (N-L+I)-th and (N-L+J)-th columns of matrices
526 * A and B: A*Q and B*Q
527 *
528  CALL crot( min( k+l, m ), a( 1, n-l+j ), 1,
529  $ a( 1, n-l+i ), 1, csq, snq )
530 *
531  CALL crot( l, b( 1, n-l+j ), 1, b( 1, n-l+i ), 1, csq,
532  $ snq )
533 *
534  IF( upper ) THEN
535  IF( k+i.LE.m )
536  $ a( k+i, n-l+j ) = czero
537  b( i, n-l+j ) = czero
538  ELSE
539  IF( k+j.LE.m )
540  $ a( k+j, n-l+i ) = czero
541  b( j, n-l+i ) = czero
542  END IF
543 *
544 * Ensure that the diagonal elements of A and B are real.
545 *
546  IF( k+i.LE.m )
547  $ a( k+i, n-l+i ) = REAL( A( K+I, N-L+I ) )
548  IF( k+j.LE.m )
549  $ a( k+j, n-l+j ) = REAL( A( K+J, N-L+J ) )
550  b( i, n-l+i ) = REAL( B( I, N-L+I ) )
551  b( j, n-l+j ) = REAL( B( J, N-L+J ) )
552 *
553 * Update unitary matrices U, V, Q, if desired.
554 *
555  IF( wantu .AND. k+j.LE.m )
556  $ CALL crot( m, u( 1, k+j ), 1, u( 1, k+i ), 1, csu,
557  $ snu )
558 *
559  IF( wantv )
560  $ CALL crot( p, v( 1, j ), 1, v( 1, i ), 1, csv, snv )
561 *
562  IF( wantq )
563  $ CALL crot( n, q( 1, n-l+j ), 1, q( 1, n-l+i ), 1, csq,
564  $ snq )
565 *
566  10 CONTINUE
567  20 CONTINUE
568 *
569  IF( .NOT.upper ) THEN
570 *
571 * The matrices A13 and B13 were lower triangular at the start
572 * of the cycle, and are now upper triangular.
573 *
574 * Convergence test: test the parallelism of the corresponding
575 * rows of A and B.
576 *
577  error = zero
578  DO 30 i = 1, min( l, m-k )
579  CALL ccopy( l-i+1, a( k+i, n-l+i ), lda, work, 1 )
580  CALL ccopy( l-i+1, b( i, n-l+i ), ldb, work( l+1 ), 1 )
581  CALL clapll( l-i+1, work, 1, work( l+1 ), 1, ssmin )
582  error = max( error, ssmin )
583  30 CONTINUE
584 *
585  IF( abs( error ).LE.min( tola, tolb ) )
586  $ GO TO 50
587  END IF
588 *
589 * End of cycle loop
590 *
591  40 CONTINUE
592 *
593 * The algorithm has not converged after MAXIT cycles.
594 *
595  info = 1
596  GO TO 100
597 *
598  50 CONTINUE
599 *
600 * If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged.
601 * Compute the generalized singular value pairs (ALPHA, BETA), and
602 * set the triangular matrix R to array A.
603 *
604  DO 60 i = 1, k
605  alpha( i ) = one
606  beta( i ) = zero
607  60 CONTINUE
608 *
609  DO 70 i = 1, min( l, m-k )
610 *
611  a1 = REAL( A( K+I, N-L+I ) )
612  b1 = REAL( B( I, N-L+I ) )
613 *
614  IF( a1.NE.zero ) THEN
615  gamma = b1 / a1
616 *
617  IF( gamma.LT.zero ) THEN
618  CALL csscal( l-i+1, -one, b( i, n-l+i ), ldb )
619  IF( wantv )
620  $ CALL csscal( p, -one, v( 1, i ), 1 )
621  END IF
622 *
623  CALL slartg( abs( gamma ), one, beta( k+i ), alpha( k+i ),
624  $ rwk )
625 *
626  IF( alpha( k+i ).GE.beta( k+i ) ) THEN
627  CALL csscal( l-i+1, one / alpha( k+i ), a( k+i, n-l+i ),
628  $ lda )
629  ELSE
630  CALL csscal( l-i+1, one / beta( k+i ), b( i, n-l+i ),
631  $ ldb )
632  CALL ccopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
633  $ lda )
634  END IF
635 *
636  ELSE
637  alpha( k+i ) = zero
638  beta( k+i ) = one
639  CALL ccopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
640  $ lda )
641  END IF
642  70 CONTINUE
643 *
644 * Post-assignment
645 *
646  DO 80 i = m + 1, k + l
647  alpha( i ) = zero
648  beta( i ) = one
649  80 CONTINUE
650 *
651  IF( k+l.LT.n ) THEN
652  DO 90 i = k + l + 1, n
653  alpha( i ) = zero
654  beta( i ) = zero
655  90 CONTINUE
656  END IF
657 *
658  100 CONTINUE
659  ncycle = kcycle
660 *
661  RETURN
662 *
663 * End of CTGSJA
664 *
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
Definition: crot.f:105
subroutine clags2(UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)
CLAGS2
Definition: clags2.f:160
subroutine clapll(N, X, INCX, Y, INCY, SSMIN)
CLAPLL measures the linear dependence of two vectors.
Definition: clapll.f:102
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99

Here is the call graph for this function:

Here is the caller graph for this function: