SRC/stgsja.f

   1 *> \brief \b STGSJA
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download STGSJA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsja.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsja.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsja.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
  22 *                          LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV,
  23 *                          Q, LDQ, WORK, NCYCLE, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBQ, JOBU, JOBV
  27 *       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N,
  28 *      $                   NCYCLE, P
  29 *       REAL               TOLA, TOLB
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       REAL               A( LDA, * ), ALPHA( * ), B( LDB, * ),
  33 *      $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
  34 *      $                   V( LDV, * ), WORK( * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> STGSJA computes the generalized singular value decomposition (GSVD)
  44 *> of two real upper triangular (or trapezoidal) matrices A and B.
  45 *>
  46 *> On entry, it is assumed that matrices A and B have the following
  47 *> forms, which may be obtained by the preprocessing subroutine SGGSVP
  48 *> from a general M-by-N matrix A and P-by-N matrix B:
  49 *>
  50 *>              N-K-L  K    L
  51 *>    A =    K ( 0    A12  A13 ) if M-K-L >= 0;
  52 *>           L ( 0     0   A23 )
  53 *>       M-K-L ( 0     0    0  )
  54 *>
  55 *>            N-K-L  K    L
  56 *>    A =  K ( 0    A12  A13 ) if M-K-L < 0;
  57 *>       M-K ( 0     0   A23 )
  58 *>
  59 *>            N-K-L  K    L
  60 *>    B =  L ( 0     0   B13 )
  61 *>       P-L ( 0     0    0  )
  62 *>
  63 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
  64 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
  65 *> otherwise A23 is (M-K)-by-L upper trapezoidal.
  66 *>
  67 *> On exit,
  68 *>
  69 *>        U**T *A*Q = D1*( 0 R ),    V**T *B*Q = D2*( 0 R ),
  70 *>
  71 *> where U, V and Q are orthogonal matrices.
  72 *> R is a nonsingular upper triangular matrix, and D1 and D2 are
  73 *> ``diagonal'' matrices, which are of the following structures:
  74 *>
  75 *> If M-K-L >= 0,
  76 *>
  77 *>                     K  L
  78 *>        D1 =     K ( I  0 )
  79 *>                 L ( 0  C )
  80 *>             M-K-L ( 0  0 )
  81 *>
  82 *>                   K  L
  83 *>        D2 = L   ( 0  S )
  84 *>             P-L ( 0  0 )
  85 *>
  86 *>                N-K-L  K    L
  87 *>   ( 0 R ) = K (  0   R11  R12 ) K
  88 *>             L (  0    0   R22 ) L
  89 *>
  90 *> where
  91 *>
  92 *>   C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
  93 *>   S = diag( BETA(K+1),  ... , BETA(K+L) ),
  94 *>   C**2 + S**2 = I.
  95 *>
  96 *>   R is stored in A(1:K+L,N-K-L+1:N) on exit.
  97 *>
  98 *> If M-K-L < 0,
  99 *>
 100 *>                K M-K K+L-M
 101 *>     D1 =   K ( I  0    0   )
 102 *>          M-K ( 0  C    0   )
 103 *>
 104 *>                  K M-K K+L-M
 105 *>     D2 =   M-K ( 0  S    0   )
 106 *>          K+L-M ( 0  0    I   )
 107 *>            P-L ( 0  0    0   )
 108 *>
 109 *>                N-K-L  K   M-K  K+L-M
 110 *> ( 0 R ) =    K ( 0    R11  R12  R13  )
 111 *>           M-K ( 0     0   R22  R23  )
 112 *>         K+L-M ( 0     0    0   R33  )
 113 *>
 114 *> where
 115 *> C = diag( ALPHA(K+1), ... , ALPHA(M) ),
 116 *> S = diag( BETA(K+1),  ... , BETA(M) ),
 117 *> C**2 + S**2 = I.
 118 *>
 119 *> R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
 120 *>     (  0  R22 R23 )
 121 *> in B(M-K+1:L,N+M-K-L+1:N) on exit.
 122 *>
 123 *> The computation of the orthogonal transformation matrices U, V or Q
 124 *> is optional.  These matrices may either be formed explicitly, or they
 125 *> may be postmultiplied into input matrices U1, V1, or Q1.
 126 *> \endverbatim
 127 *
 128 *  Arguments:
 129 *  ==========
 130 *
 131 *> \param[in] JOBU
 132 *> \verbatim
 133 *>          JOBU is CHARACTER*1
 134 *>          = 'U':  U must contain an orthogonal matrix U1 on entry, and
 135 *>                  the product U1*U is returned;
 136 *>          = 'I':  U is initialized to the unit matrix, and the
 137 *>                  orthogonal matrix U is returned;
 138 *>          = 'N':  U is not computed.
 139 *> \endverbatim
 140 *>
 141 *> \param[in] JOBV
 142 *> \verbatim
 143 *>          JOBV is CHARACTER*1
 144 *>          = 'V':  V must contain an orthogonal matrix V1 on entry, and
 145 *>                  the product V1*V is returned;
 146 *>          = 'I':  V is initialized to the unit matrix, and the
 147 *>                  orthogonal matrix V is returned;
 148 *>          = 'N':  V is not computed.
 149 *> \endverbatim
 150 *>
 151 *> \param[in] JOBQ
 152 *> \verbatim
 153 *>          JOBQ is CHARACTER*1
 154 *>          = 'Q':  Q must contain an orthogonal matrix Q1 on entry, and
 155 *>                  the product Q1*Q is returned;
 156 *>          = 'I':  Q is initialized to the unit matrix, and the
 157 *>                  orthogonal matrix Q is returned;
 158 *>          = 'N':  Q is not computed.
 159 *> \endverbatim
 160 *>
 161 *> \param[in] M
 162 *> \verbatim
 163 *>          M is INTEGER
 164 *>          The number of rows of the matrix A.  M >= 0.
 165 *> \endverbatim
 166 *>
 167 *> \param[in] P
 168 *> \verbatim
 169 *>          P is INTEGER
 170 *>          The number of rows of the matrix B.  P >= 0.
 171 *> \endverbatim
 172 *>
 173 *> \param[in] N
 174 *> \verbatim
 175 *>          N is INTEGER
 176 *>          The number of columns of the matrices A and B.  N >= 0.
 177 *> \endverbatim
 178 *>
 179 *> \param[in] K
 180 *> \verbatim
 181 *>          K is INTEGER
 182 *> \endverbatim
 183 *>
 184 *> \param[in] L
 185 *> \verbatim
 186 *>          L is INTEGER
 187 *>
 188 *>          K and L specify the subblocks in the input matrices A and B:
 189 *>          A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N)
 190 *>          of A and B, whose GSVD is going to be computed by STGSJA.
 191 *>          See Further Details.
 192 *> \endverbatim
 193 *>
 194 *> \param[in,out] A
 195 *> \verbatim
 196 *>          A is REAL array, dimension (LDA,N)
 197 *>          On entry, the M-by-N matrix A.
 198 *>          On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular
 199 *>          matrix R or part of R.  See Purpose for details.
 200 *> \endverbatim
 201 *>
 202 *> \param[in] LDA
 203 *> \verbatim
 204 *>          LDA is INTEGER
 205 *>          The leading dimension of the array A. LDA >= max(1,M).
 206 *> \endverbatim
 207 *>
 208 *> \param[in,out] B
 209 *> \verbatim
 210 *>          B is REAL array, dimension (LDB,N)
 211 *>          On entry, the P-by-N matrix B.
 212 *>          On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
 213 *>          a part of R.  See Purpose for details.
 214 *> \endverbatim
 215 *>
 216 *> \param[in] LDB
 217 *> \verbatim
 218 *>          LDB is INTEGER
 219 *>          The leading dimension of the array B. LDB >= max(1,P).
 220 *> \endverbatim
 221 *>
 222 *> \param[in] TOLA
 223 *> \verbatim
 224 *>          TOLA is REAL
 225 *> \endverbatim
 226 *>
 227 *> \param[in] TOLB
 228 *> \verbatim
 229 *>          TOLB is REAL
 230 *>
 231 *>          TOLA and TOLB are the convergence criteria for the Jacobi-
 232 *>          Kogbetliantz iteration procedure. Generally, they are the
 233 *>          same as used in the preprocessing step, say
 234 *>              TOLA = max(M,N)*norm(A)*MACHEPS,
 235 *>              TOLB = max(P,N)*norm(B)*MACHEPS.
 236 *> \endverbatim
 237 *>
 238 *> \param[out] ALPHA
 239 *> \verbatim
 240 *>          ALPHA is REAL array, dimension (N)
 241 *> \endverbatim
 242 *>
 243 *> \param[out] BETA
 244 *> \verbatim
 245 *>          BETA is REAL array, dimension (N)
 246 *>
 247 *>          On exit, ALPHA and BETA contain the generalized singular
 248 *>          value pairs of A and B;
 249 *>            ALPHA(1:K) = 1,
 250 *>            BETA(1:K)  = 0,
 251 *>          and if M-K-L >= 0,
 252 *>            ALPHA(K+1:K+L) = diag(C),
 253 *>            BETA(K+1:K+L)  = diag(S),
 254 *>          or if M-K-L < 0,
 255 *>            ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0
 256 *>            BETA(K+1:M) = S, BETA(M+1:K+L) = 1.
 257 *>          Furthermore, if K+L < N,
 258 *>            ALPHA(K+L+1:N) = 0 and
 259 *>            BETA(K+L+1:N)  = 0.
 260 *> \endverbatim
 261 *>
 262 *> \param[in,out] U
 263 *> \verbatim
 264 *>          U is REAL array, dimension (LDU,M)
 265 *>          On entry, if JOBU = 'U', U must contain a matrix U1 (usually
 266 *>          the orthogonal matrix returned by SGGSVP).
 267 *>          On exit,
 268 *>          if JOBU = 'I', U contains the orthogonal matrix U;
 269 *>          if JOBU = 'U', U contains the product U1*U.
 270 *>          If JOBU = 'N', U is not referenced.
 271 *> \endverbatim
 272 *>
 273 *> \param[in] LDU
 274 *> \verbatim
 275 *>          LDU is INTEGER
 276 *>          The leading dimension of the array U. LDU >= max(1,M) if
 277 *>          JOBU = 'U'; LDU >= 1 otherwise.
 278 *> \endverbatim
 279 *>
 280 *> \param[in,out] V
 281 *> \verbatim
 282 *>          V is REAL array, dimension (LDV,P)
 283 *>          On entry, if JOBV = 'V', V must contain a matrix V1 (usually
 284 *>          the orthogonal matrix returned by SGGSVP).
 285 *>          On exit,
 286 *>          if JOBV = 'I', V contains the orthogonal matrix V;
 287 *>          if JOBV = 'V', V contains the product V1*V.
 288 *>          If JOBV = 'N', V is not referenced.
 289 *> \endverbatim
 290 *>
 291 *> \param[in] LDV
 292 *> \verbatim
 293 *>          LDV is INTEGER
 294 *>          The leading dimension of the array V. LDV >= max(1,P) if
 295 *>          JOBV = 'V'; LDV >= 1 otherwise.
 296 *> \endverbatim
 297 *>
 298 *> \param[in,out] Q
 299 *> \verbatim
 300 *>          Q is REAL array, dimension (LDQ,N)
 301 *>          On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
 302 *>          the orthogonal matrix returned by SGGSVP).
 303 *>          On exit,
 304 *>          if JOBQ = 'I', Q contains the orthogonal matrix Q;
 305 *>          if JOBQ = 'Q', Q contains the product Q1*Q.
 306 *>          If JOBQ = 'N', Q is not referenced.
 307 *> \endverbatim
 308 *>
 309 *> \param[in] LDQ
 310 *> \verbatim
 311 *>          LDQ is INTEGER
 312 *>          The leading dimension of the array Q. LDQ >= max(1,N) if
 313 *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
 314 *> \endverbatim
 315 *>
 316 *> \param[out] WORK
 317 *> \verbatim
 318 *>          WORK is REAL array, dimension (2*N)
 319 *> \endverbatim
 320 *>
 321 *> \param[out] NCYCLE
 322 *> \verbatim
 323 *>          NCYCLE is INTEGER
 324 *>          The number of cycles required for convergence.
 325 *> \endverbatim
 326 *>
 327 *> \param[out] INFO
 328 *> \verbatim
 329 *>          INFO is INTEGER
 330 *>          = 0:  successful exit
 331 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 332 *>          = 1:  the procedure does not converge after MAXIT cycles.
 333 *> \endverbatim
 334 *>
 335 *> \verbatim
 336 *>  Internal Parameters
 337 *>  ===================
 338 *>
 339 *>  MAXIT   INTEGER
 340 *>          MAXIT specifies the total loops that the iterative procedure
 341 *>          may take. If after MAXIT cycles, the routine fails to
 342 *>          converge, we return INFO = 1.
 343 *> \endverbatim
 344 *
 345 *  Authors:
 346 *  ========
 347 *
 348 *> \author Univ. of Tennessee
 349 *> \author Univ. of California Berkeley
 350 *> \author Univ. of Colorado Denver
 351 *> \author NAG Ltd.
 352 *
 353 *> \date November 2011
 354 *
 355 *> \ingroup realOTHERcomputational
 356 *
 357 *> \par Further Details:
 358 *  =====================
 359 *>
 360 *> \verbatim
 361 *>
 362 *>  STGSJA essentially uses a variant of Kogbetliantz algorithm to reduce
 363 *>  min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
 364 *>  matrix B13 to the form:
 365 *>
 366 *>           U1**T *A13*Q1 = C1*R1; V1**T *B13*Q1 = S1*R1,
 367 *>
 368 *>  where U1, V1 and Q1 are orthogonal matrix, and Z**T is the transpose
 369 *>  of Z.  C1 and S1 are diagonal matrices satisfying
 370 *>
 371 *>                C1**2 + S1**2 = I,
 372 *>
 373 *>  and R1 is an L-by-L nonsingular upper triangular matrix.
 374 *> \endverbatim
 375 *>
 376 *  =====================================================================
 377       SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
 378      $                   LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV,
 379      $                   Q, LDQ, WORK, NCYCLE, INFO )
 380 *
 381 *  -- LAPACK computational routine (version 3.4.0) --
 382 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 383 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 384 *     November 2011
 385 *
 386 *     .. Scalar Arguments ..
 387       CHARACTER          JOBQ, JOBU, JOBV
 388       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N,
 389      $                   NCYCLE, P
 390       REAL               TOLA, TOLB
 391 *     ..
 392 *     .. Array Arguments ..
 393       REAL               A( LDA, * ), ALPHA( * ), B( LDB, * ),
 394      $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
 395      $                   V( LDV, * ), WORK( * )
 396 *     ..
 397 *
 398 *  =====================================================================
 399 *
 400 *     .. Parameters ..
 401       INTEGER            MAXIT
 402       PARAMETER          ( MAXIT = 40 )
 403       REAL               ZERO, ONE
 404       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 405 *     ..
 406 *     .. Local Scalars ..
 407 *
 408       LOGICAL            INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
 409       INTEGER            I, J, KCYCLE
 410       REAL               A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR,
 411      $                   GAMMA, RWK, SNQ, SNU, SNV, SSMIN
 412 *     ..
 413 *     .. External Functions ..
 414       LOGICAL            LSAME
 415       EXTERNAL           LSAME
 416 *     ..
 417 *     .. External Subroutines ..
 418       EXTERNAL           SCOPY, SLAGS2, SLAPLL, SLARTG, SLASET, SROT,
 419      $                   SSCAL, XERBLA
 420 *     ..
 421 *     .. Intrinsic Functions ..
 422       INTRINSIC          ABS, MAX, MIN
 423 *     ..
 424 *     .. Executable Statements ..
 425 *
 426 *     Decode and test the input parameters
 427 *
 428       INITU = LSAME( JOBU, 'I' )
 429       WANTU = INITU .OR. LSAME( JOBU, 'U' )
 430 *
 431       INITV = LSAME( JOBV, 'I' )
 432       WANTV = INITV .OR. LSAME( JOBV, 'V' )
 433 *
 434       INITQ = LSAME( JOBQ, 'I' )
 435       WANTQ = INITQ .OR. LSAME( JOBQ, 'Q' )
 436 *
 437       INFO = 0
 438       IF( .NOT.( INITU .OR. WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
 439          INFO = -1
 440       ELSE IF( .NOT.( INITV .OR. WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
 441          INFO = -2
 442       ELSE IF( .NOT.( INITQ .OR. WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
 443          INFO = -3
 444       ELSE IF( M.LT.0 ) THEN
 445          INFO = -4
 446       ELSE IF( P.LT.0 ) THEN
 447          INFO = -5
 448       ELSE IF( N.LT.0 ) THEN
 449          INFO = -6
 450       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 451          INFO = -10
 452       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
 453          INFO = -12
 454       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
 455          INFO = -18
 456       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
 457          INFO = -20
 458       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
 459          INFO = -22
 460       END IF
 461       IF( INFO.NE.0 ) THEN
 462          CALL XERBLA( 'STGSJA', -INFO )
 463          RETURN
 464       END IF
 465 *
 466 *     Initialize U, V and Q, if necessary
 467 *
 468       IF( INITU )
 469      $   CALL SLASET( 'Full', M, M, ZERO, ONE, U, LDU )
 470       IF( INITV )
 471      $   CALL SLASET( 'Full', P, P, ZERO, ONE, V, LDV )
 472       IF( INITQ )
 473      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
 474 *
 475 *     Loop until convergence
 476 *
 477       UPPER = .FALSE.
 478       DO 40 KCYCLE = 1, MAXIT
 479 *
 480          UPPER = .NOT.UPPER
 481 *
 482          DO 20 I = 1, L - 1
 483             DO 10 J = I + 1, L
 484 *
 485                A1 = ZERO
 486                A2 = ZERO
 487                A3 = ZERO
 488                IF( K+I.LE.M )
 489      $            A1 = A( K+I, N-L+I )
 490                IF( K+J.LE.M )
 491      $            A3 = A( K+J, N-L+J )
 492 *
 493                B1 = B( I, N-L+I )
 494                B3 = B( J, N-L+J )
 495 *
 496                IF( UPPER ) THEN
 497                   IF( K+I.LE.M )
 498      $               A2 = A( K+I, N-L+J )
 499                   B2 = B( I, N-L+J )
 500                ELSE
 501                   IF( K+J.LE.M )
 502      $               A2 = A( K+J, N-L+I )
 503                   B2 = B( J, N-L+I )
 504                END IF
 505 *
 506                CALL SLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU,
 507      $                      CSV, SNV, CSQ, SNQ )
 508 *
 509 *              Update (K+I)-th and (K+J)-th rows of matrix A: U**T *A
 510 *
 511                IF( K+J.LE.M )
 512      $            CALL SROT( L, A( K+J, N-L+1 ), LDA, A( K+I, N-L+1 ),
 513      $                       LDA, CSU, SNU )
 514 *
 515 *              Update I-th and J-th rows of matrix B: V**T *B
 516 *
 517                CALL SROT( L, B( J, N-L+1 ), LDB, B( I, N-L+1 ), LDB,
 518      $                    CSV, SNV )
 519 *
 520 *              Update (N-L+I)-th and (N-L+J)-th columns of matrices
 521 *              A and B: A*Q and B*Q
 522 *
 523                CALL SROT( MIN( K+L, M ), A( 1, N-L+J ), 1,
 524      $                    A( 1, N-L+I ), 1, CSQ, SNQ )
 525 *
 526                CALL SROT( L, B( 1, N-L+J ), 1, B( 1, N-L+I ), 1, CSQ,
 527      $                    SNQ )
 528 *
 529                IF( UPPER ) THEN
 530                   IF( K+I.LE.M )
 531      $               A( K+I, N-L+J ) = ZERO
 532                   B( I, N-L+J ) = ZERO
 533                ELSE
 534                   IF( K+J.LE.M )
 535      $               A( K+J, N-L+I ) = ZERO
 536                   B( J, N-L+I ) = ZERO
 537                END IF
 538 *
 539 *              Update orthogonal matrices U, V, Q, if desired.
 540 *
 541                IF( WANTU .AND. K+J.LE.M )
 542      $            CALL SROT( M, U( 1, K+J ), 1, U( 1, K+I ), 1, CSU,
 543      $                       SNU )
 544 *
 545                IF( WANTV )
 546      $            CALL SROT( P, V( 1, J ), 1, V( 1, I ), 1, CSV, SNV )
 547 *
 548                IF( WANTQ )
 549      $            CALL SROT( N, Q( 1, N-L+J ), 1, Q( 1, N-L+I ), 1, CSQ,
 550      $                       SNQ )
 551 *
 552    10       CONTINUE
 553    20    CONTINUE
 554 *
 555          IF( .NOT.UPPER ) THEN
 556 *
 557 *           The matrices A13 and B13 were lower triangular at the start
 558 *           of the cycle, and are now upper triangular.
 559 *
 560 *           Convergence test: test the parallelism of the corresponding
 561 *           rows of A and B.
 562 *
 563             ERROR = ZERO
 564             DO 30 I = 1, MIN( L, M-K )
 565                CALL SCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 )
 566                CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 )
 567                CALL SLAPLL( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN )
 568                ERROR = MAX( ERROR, SSMIN )
 569    30       CONTINUE
 570 *
 571             IF( ABS( ERROR ).LE.MIN( TOLA, TOLB ) )
 572      $         GO TO 50
 573          END IF
 574 *
 575 *        End of cycle loop
 576 *
 577    40 CONTINUE
 578 *
 579 *     The algorithm has not converged after MAXIT cycles.
 580 *
 581       INFO = 1
 582       GO TO 100
 583 *
 584    50 CONTINUE
 585 *
 586 *     If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged.
 587 *     Compute the generalized singular value pairs (ALPHA, BETA), and
 588 *     set the triangular matrix R to array A.
 589 *
 590       DO 60 I = 1, K
 591          ALPHA( I ) = ONE
 592          BETA( I ) = ZERO
 593    60 CONTINUE
 594 *
 595       DO 70 I = 1, MIN( L, M-K )
 596 *
 597          A1 = A( K+I, N-L+I )
 598          B1 = B( I, N-L+I )
 599 *
 600          IF( A1.NE.ZERO ) THEN
 601             GAMMA = B1 / A1
 602 *
 603 *           change sign if necessary
 604 *
 605             IF( GAMMA.LT.ZERO ) THEN
 606                CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
 607                IF( WANTV )
 608      $            CALL SSCAL( P, -ONE, V( 1, I ), 1 )
 609             END IF
 610 *
 611             CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
 612      $                   RWK )
 613 *
 614             IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
 615                CALL SSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ),
 616      $                     LDA )
 617             ELSE
 618                CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
 619      $                     LDB )
 620                CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
 621      $                     LDA )
 622             END IF
 623 *
 624          ELSE
 625 *
 626             ALPHA( K+I ) = ZERO
 627             BETA( K+I ) = ONE
 628             CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
 629      $                  LDA )
 630 *
 631          END IF
 632 *
 633    70 CONTINUE
 634 *
 635 *     Post-assignment
 636 *
 637       DO 80 I = M + 1, K + L
 638          ALPHA( I ) = ZERO
 639          BETA( I ) = ONE
 640    80 CONTINUE
 641 *
 642       IF( K+L.LT.N ) THEN
 643          DO 90 I = K + L + 1, N
 644             ALPHA( I ) = ZERO
 645             BETA( I ) = ZERO
 646    90    CONTINUE
 647       END IF
 648 *
 649   100 CONTINUE
 650       NCYCLE = KCYCLE
 651       RETURN
 652 *
 653 *     End of STGSJA
 654 *
 655       END