stgsja()

subroutine stgsja	(	character	jobu,
		character	jobv,
		character	jobq,
		integer	m,
		integer	p,
		integer	n,
		integer	k,
		integer	l,
		real, dimension( lda, * )	a,
		integer	lda,
		real, dimension( ldb, * )	b,
		integer	ldb,
		real	tola,
		real	tolb,
		real, dimension( * )	alpha,
		real, dimension( * )	beta,
		real, dimension( ldu, * )	u,
		integer	ldu,
		real, dimension( ldv, * )	v,
		integer	ldv,
		real, dimension( ldq, * )	q,
		integer	ldq,
		real, dimension( * )	work,
		integer	ncycle,
		integer	info )

STGSJA

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

Purpose:

!>
!> STGSJA computes the generalized singular value decomposition (GSVD)
!> of two real 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 SGGSVP
!> 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**T *A*Q = D1*( 0 R ),    V**T *B*Q = D2*( 0 R ),
!>
!> where U, V and Q are orthogonal 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 orthogonal 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 an orthogonal matrix U1 on entry, and !> the product U1*U is returned; !> = 'I': U is initialized to the unit matrix, and the !> orthogonal matrix U is returned; !> = 'N': U is not computed. !>
[in]	JOBV	!> JOBV is CHARACTER*1 !> = 'V': V must contain an orthogonal matrix V1 on entry, and !> the product V1*V is returned; !> = 'I': V is initialized to the unit matrix, and the !> orthogonal matrix V is returned; !> = 'N': V is not computed. !>
[in]	JOBQ	!> JOBQ is CHARACTER*1 !> = 'Q': Q must contain an orthogonal matrix Q1 on entry, and !> the product Q1*Q is returned; !> = 'I': Q is initialized to the unit matrix, and the !> orthogonal 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 STGSJA. !> See Further Details. !>
[in,out]	A	!> A is REAL 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 REAL 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 and !> BETA(K+L+1:N) = 0. !>
[in,out]	U	!> U is REAL array, dimension (LDU,M) !> On entry, if JOBU = 'U', U must contain a matrix U1 (usually !> the orthogonal matrix returned by SGGSVP). !> On exit, !> if JOBU = 'I', U contains the orthogonal 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 REAL array, dimension (LDV,P) !> On entry, if JOBV = 'V', V must contain a matrix V1 (usually !> the orthogonal matrix returned by SGGSVP). !> On exit, !> if JOBV = 'I', V contains the orthogonal 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 REAL array, dimension (LDQ,N) !> On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually !> the orthogonal matrix returned by SGGSVP). !> On exit, !> if JOBQ = 'I', Q contains the orthogonal 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 REAL 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.

Further Details:

!>
!>  STGSJA 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**T *A13*Q1 = C1*R1; V1**T *B13*Q1 = S1*R1,
!>
!>  where U1, V1 and Q1 are orthogonal matrix, and Z**T is the transpose
!>  of Z.  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 373 of file stgsja.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          JOBQ, JOBU, JOBV
      INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N,
     $                   NCYCLE, P
      REAL               TOLA, TOLB
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * ), ALPHA( * ), B( LDB, * ),
     $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
     $                   V( LDV, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            MAXIT
      parameter( maxit = 40 )
      REAL               ZERO, ONE, HUGENUM
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
*
      LOGICAL            INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
      INTEGER            I, J, KCYCLE
      REAL               A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR,
     $                   GAMMA, RWK, SNQ, SNU, SNV, SSMIN
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, slags2, slapll, slartg, slaset,
     $                   srot,
     $                   sscal, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, huge
      parameter( hugenum = huge(zero) )
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input parameters
*
      initu = lsame( jobu, 'I' )
      wantu = initu .OR. lsame( jobu, 'U' )
*
      initv = lsame( jobv, 'I' )
      wantv = initv .OR. lsame( jobv, 'V' )
*
      initq = lsame( jobq, 'I' )
      wantq = initq .OR. lsame( jobq, 'Q' )
*
      info = 0
      IF( .NOT.( initu .OR. wantu .OR. lsame( jobu, 'N' ) ) ) THEN
         info = -1
      ELSE IF( .NOT.( initv .OR.
     $         wantv .OR.
     $         lsame( jobv, 'N' ) ) ) THEN
         info = -2
      ELSE IF( .NOT.( initq .OR.
     $         wantq .OR.
     $         lsame( jobq, 'N' ) ) ) THEN
         info = -3
      ELSE IF( m.LT.0 ) THEN
         info = -4
      ELSE IF( p.LT.0 ) THEN
         info = -5
      ELSE IF( n.LT.0 ) THEN
         info = -6
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -10
      ELSE IF( ldb.LT.max( 1, p ) ) THEN
         info = -12
      ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
         info = -18
      ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
         info = -20
      ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
         info = -22
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'STGSJA', -info )
         RETURN
      END IF
*
*     Initialize U, V and Q, if necessary
*
      IF( initu )
     $   CALL slaset( 'Full', m, m, zero, one, u, ldu )
      IF( initv )
     $   CALL slaset( 'Full', p, p, zero, one, v, ldv )
      IF( initq )
     $   CALL slaset( 'Full', n, n, zero, one, q, ldq )
*
*     Loop until convergence
*
      upper = .false.
      DO 40 kcycle = 1, maxit
*
         upper = .NOT.upper
*
         DO 20 i = 1, l - 1
            DO 10 j = i + 1, l
*
               a1 = zero
               a2 = zero
               a3 = zero
               IF( k+i.LE.m )
     $            a1 = a( k+i, n-l+i )
               IF( k+j.LE.m )
     $            a3 = a( k+j, n-l+j )
*
               b1 = b( i, n-l+i )
               b3 = b( j, n-l+j )
*
               IF( upper ) THEN
                  IF( k+i.LE.m )
     $               a2 = a( k+i, n-l+j )
                  b2 = b( i, n-l+j )
               ELSE
                  IF( k+j.LE.m )
     $               a2 = a( k+j, n-l+i )
                  b2 = b( j, n-l+i )
               END IF
*
               CALL slags2( upper, a1, a2, a3, b1, b2, b3, csu, snu,
     $                      csv, snv, csq, snq )
*
*              Update (K+I)-th and (K+J)-th rows of matrix A: U**T *A
*
               IF( k+j.LE.m )
     $            CALL srot( l, a( k+j, n-l+1 ), lda, a( k+i,
     $                       n-l+1 ),
     $                       lda, csu, snu )
*
*              Update I-th and J-th rows of matrix B: V**T *B
*
               CALL srot( l, b( j, n-l+1 ), ldb, b( i, n-l+1 ), ldb,
     $                    csv, snv )
*
*              Update (N-L+I)-th and (N-L+J)-th columns of matrices
*              A and B: A*Q and B*Q
*
               CALL srot( min( k+l, m ), a( 1, n-l+j ), 1,
     $                    a( 1, n-l+i ), 1, csq, snq )
*
               CALL srot( l, b( 1, n-l+j ), 1, b( 1, n-l+i ), 1, csq,
     $                    snq )
*
               IF( upper ) THEN
                  IF( k+i.LE.m )
     $               a( k+i, n-l+j ) = zero
                  b( i, n-l+j ) = zero
               ELSE
                  IF( k+j.LE.m )
     $               a( k+j, n-l+i ) = zero
                  b( j, n-l+i ) = zero
               END IF
*
*              Update orthogonal matrices U, V, Q, if desired.
*
               IF( wantu .AND. k+j.LE.m )
     $            CALL srot( m, u( 1, k+j ), 1, u( 1, k+i ), 1, csu,
     $                       snu )
*
               IF( wantv )
     $            CALL srot( p, v( 1, j ), 1, v( 1, i ), 1, csv,
     $                       snv )
*
               IF( wantq )
     $            CALL srot( n, q( 1, n-l+j ), 1, q( 1, n-l+i ), 1,
     $                       csq,
     $                       snq )
*
   10       CONTINUE
   20    CONTINUE
*
         IF( .NOT.upper ) THEN
*
*           The matrices A13 and B13 were lower triangular at the start
*           of the cycle, and are now upper triangular.
*
*           Convergence test: test the parallelism of the corresponding
*           rows of A and B.
*
            error = zero
            DO 30 i = 1, min( l, m-k )
               CALL scopy( l-i+1, a( k+i, n-l+i ), lda, work, 1 )
               CALL scopy( l-i+1, b( i, n-l+i ), ldb, work( l+1 ),
     $                     1 )
               CALL slapll( l-i+1, work, 1, work( l+1 ), 1, ssmin )
               error = max( error, ssmin )
   30       CONTINUE
*
            IF( abs( error ).LE.min( tola, tolb ) )
     $         GO TO 50
         END IF
*
*        End of cycle loop
*
   40 CONTINUE
*
*     The algorithm has not converged after MAXIT cycles.
*
      info = 1
      GO TO 100
*
   50 CONTINUE
*
*     If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged.
*     Compute the generalized singular value pairs (ALPHA, BETA), and
*     set the triangular matrix R to array A.
*
      DO 60 i = 1, k
         alpha( i ) = one
         beta( i ) = zero
   60 CONTINUE
*
      DO 70 i = 1, min( l, m-k )
*
         a1 = a( k+i, n-l+i )
         b1 = b( i, n-l+i )
         gamma = b1 / a1
*
         IF( (gamma.LE.hugenum).AND.(gamma.GE.-hugenum) ) THEN
*
*           change sign if necessary
*
            IF( gamma.LT.zero ) THEN
               CALL sscal( l-i+1, -one, b( i, n-l+i ), ldb )
               IF( wantv )
     $            CALL sscal( p, -one, v( 1, i ), 1 )
            END IF
*
            CALL slartg( abs( gamma ), one, beta( k+i ),
     $                   alpha( k+i ),
     $                   rwk )
*
            IF( alpha( k+i ).GE.beta( k+i ) ) THEN
               CALL sscal( l-i+1, one / alpha( k+i ), a( k+i,
     $                     n-l+i ),
     $                     lda )
            ELSE
               CALL sscal( l-i+1, one / beta( k+i ), b( i, n-l+i ),
     $                     ldb )
               CALL scopy( l-i+1, b( i, n-l+i ), ldb, a( k+i,
     $                     n-l+i ),
     $                     lda )
            END IF
*
         ELSE
*
            alpha( k+i ) = zero
            beta( k+i ) = one
            CALL scopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
     $                  lda )
*
         END IF
*
   70 CONTINUE
*
*     Post-assignment
*
      DO 80 i = m + 1, k + l
         alpha( i ) = zero
         beta( i ) = one
   80 CONTINUE
*
      IF( k+l.LT.n ) THEN
         DO 90 i = k + l + 1, n
            alpha( i ) = zero
            beta( i ) = zero
   90    CONTINUE
      END IF
*
  100 CONTINUE
      ncycle = kcycle
      RETURN
*
*     End of STGSJA
*

Here is the call graph for this function:

Here is the caller graph for this function: