stgsja.f

*> \brief \b STGSJA
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download STGSJA + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsja.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsja.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsja.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE STGSJA( 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 )
*
*       .. 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( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOBU
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] JOBV
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] JOBQ
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix A.  M >= 0.
*> \endverbatim
*>
*> \param[in] P
*> \verbatim
*>          P is INTEGER
*>          The number of rows of the matrix B.  P >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrices A and B.  N >= 0.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*>          K is INTEGER
*> \endverbatim
*>
*> \param[in] L
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B. LDB >= max(1,P).
*> \endverbatim
*>
*> \param[in] TOLA
*> \verbatim
*>          TOLA is REAL
*> \endverbatim
*>
*> \param[in] TOLB
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[out] ALPHA
*> \verbatim
*>          ALPHA is REAL array, dimension (N)
*> \endverbatim
*>
*> \param[out] BETA
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in,out] U
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER
*>          The leading dimension of the array U. LDU >= max(1,M) if
*>          JOBU = 'U'; LDU >= 1 otherwise.
*> \endverbatim
*>
*> \param[in,out] V
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDV
*> \verbatim
*>          LDV is INTEGER
*>          The leading dimension of the array V. LDV >= max(1,P) if
*>          JOBV = 'V'; LDV >= 1 otherwise.
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q. LDQ >= max(1,N) if
*>          JOBQ = 'Q'; LDQ >= 1 otherwise.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (2*N)
*> \endverbatim
*>
*> \param[out] NCYCLE
*> \verbatim
*>          NCYCLE is INTEGER
*>          The number of cycles required for convergence.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \verbatim
*>  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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup tgsja
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  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.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE stgsja( 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 )
*
*  -- 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
*
      END