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 December 2016 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.7.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* December 2016 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, HUGENUM 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, HUGE 423 PARAMETER ( HUGENUM = HUGE(ZERO) ) 424* .. 425* .. Executable Statements .. 426* 427* Decode and test the input parameters 428* 429 INITU = LSAME( JOBU, 'I' ) 430 WANTU = INITU .OR. LSAME( JOBU, 'U' ) 431* 432 INITV = LSAME( JOBV, 'I' ) 433 WANTV = INITV .OR. LSAME( JOBV, 'V' ) 434* 435 INITQ = LSAME( JOBQ, 'I' ) 436 WANTQ = INITQ .OR. LSAME( JOBQ, 'Q' ) 437* 438 INFO = 0 439 IF( .NOT.( INITU .OR. WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 440 INFO = -1 441 ELSE IF( .NOT.( INITV .OR. WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 442 INFO = -2 443 ELSE IF( .NOT.( INITQ .OR. WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 444 INFO = -3 445 ELSE IF( M.LT.0 ) THEN 446 INFO = -4 447 ELSE IF( P.LT.0 ) THEN 448 INFO = -5 449 ELSE IF( N.LT.0 ) THEN 450 INFO = -6 451 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 452 INFO = -10 453 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 454 INFO = -12 455 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 456 INFO = -18 457 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 458 INFO = -20 459 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 460 INFO = -22 461 END IF 462 IF( INFO.NE.0 ) THEN 463 CALL XERBLA( 'STGSJA', -INFO ) 464 RETURN 465 END IF 466* 467* Initialize U, V and Q, if necessary 468* 469 IF( INITU ) 470 $ CALL SLASET( 'Full', M, M, ZERO, ONE, U, LDU ) 471 IF( INITV ) 472 $ CALL SLASET( 'Full', P, P, ZERO, ONE, V, LDV ) 473 IF( INITQ ) 474 $ CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 475* 476* Loop until convergence 477* 478 UPPER = .FALSE. 479 DO 40 KCYCLE = 1, MAXIT 480* 481 UPPER = .NOT.UPPER 482* 483 DO 20 I = 1, L - 1 484 DO 10 J = I + 1, L 485* 486 A1 = ZERO 487 A2 = ZERO 488 A3 = ZERO 489 IF( K+I.LE.M ) 490 $ A1 = A( K+I, N-L+I ) 491 IF( K+J.LE.M ) 492 $ A3 = A( K+J, N-L+J ) 493* 494 B1 = B( I, N-L+I ) 495 B3 = B( J, N-L+J ) 496* 497 IF( UPPER ) THEN 498 IF( K+I.LE.M ) 499 $ A2 = A( K+I, N-L+J ) 500 B2 = B( I, N-L+J ) 501 ELSE 502 IF( K+J.LE.M ) 503 $ A2 = A( K+J, N-L+I ) 504 B2 = B( J, N-L+I ) 505 END IF 506* 507 CALL SLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, 508 $ CSV, SNV, CSQ, SNQ ) 509* 510* Update (K+I)-th and (K+J)-th rows of matrix A: U**T *A 511* 512 IF( K+J.LE.M ) 513 $ CALL SROT( L, A( K+J, N-L+1 ), LDA, A( K+I, N-L+1 ), 514 $ LDA, CSU, SNU ) 515* 516* Update I-th and J-th rows of matrix B: V**T *B 517* 518 CALL SROT( L, B( J, N-L+1 ), LDB, B( I, N-L+1 ), LDB, 519 $ CSV, SNV ) 520* 521* Update (N-L+I)-th and (N-L+J)-th columns of matrices 522* A and B: A*Q and B*Q 523* 524 CALL SROT( MIN( K+L, M ), A( 1, N-L+J ), 1, 525 $ A( 1, N-L+I ), 1, CSQ, SNQ ) 526* 527 CALL SROT( L, B( 1, N-L+J ), 1, B( 1, N-L+I ), 1, CSQ, 528 $ SNQ ) 529* 530 IF( UPPER ) THEN 531 IF( K+I.LE.M ) 532 $ A( K+I, N-L+J ) = ZERO 533 B( I, N-L+J ) = ZERO 534 ELSE 535 IF( K+J.LE.M ) 536 $ A( K+J, N-L+I ) = ZERO 537 B( J, N-L+I ) = ZERO 538 END IF 539* 540* Update orthogonal matrices U, V, Q, if desired. 541* 542 IF( WANTU .AND. K+J.LE.M ) 543 $ CALL SROT( M, U( 1, K+J ), 1, U( 1, K+I ), 1, CSU, 544 $ SNU ) 545* 546 IF( WANTV ) 547 $ CALL SROT( P, V( 1, J ), 1, V( 1, I ), 1, CSV, SNV ) 548* 549 IF( WANTQ ) 550 $ CALL SROT( N, Q( 1, N-L+J ), 1, Q( 1, N-L+I ), 1, CSQ, 551 $ SNQ ) 552* 553 10 CONTINUE 554 20 CONTINUE 555* 556 IF( .NOT.UPPER ) THEN 557* 558* The matrices A13 and B13 were lower triangular at the start 559* of the cycle, and are now upper triangular. 560* 561* Convergence test: test the parallelism of the corresponding 562* rows of A and B. 563* 564 ERROR = ZERO 565 DO 30 I = 1, MIN( L, M-K ) 566 CALL SCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 ) 567 CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 ) 568 CALL SLAPLL( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN ) 569 ERROR = MAX( ERROR, SSMIN ) 570 30 CONTINUE 571* 572 IF( ABS( ERROR ).LE.MIN( TOLA, TOLB ) ) 573 $ GO TO 50 574 END IF 575* 576* End of cycle loop 577* 578 40 CONTINUE 579* 580* The algorithm has not converged after MAXIT cycles. 581* 582 INFO = 1 583 GO TO 100 584* 585 50 CONTINUE 586* 587* If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged. 588* Compute the generalized singular value pairs (ALPHA, BETA), and 589* set the triangular matrix R to array A. 590* 591 DO 60 I = 1, K 592 ALPHA( I ) = ONE 593 BETA( I ) = ZERO 594 60 CONTINUE 595* 596 DO 70 I = 1, MIN( L, M-K ) 597* 598 A1 = A( K+I, N-L+I ) 599 B1 = B( I, N-L+I ) 600 GAMMA = B1 / A1 601* 602 IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN 603* 604* change sign if necessary 605* 606 IF( GAMMA.LT.ZERO ) THEN 607 CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) 608 IF( WANTV ) 609 $ CALL SSCAL( P, -ONE, V( 1, I ), 1 ) 610 END IF 611* 612 CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), 613 $ RWK ) 614* 615 IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN 616 CALL SSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), 617 $ LDA ) 618 ELSE 619 CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), 620 $ LDB ) 621 CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 622 $ LDA ) 623 END IF 624* 625 ELSE 626* 627 ALPHA( K+I ) = ZERO 628 BETA( K+I ) = ONE 629 CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 630 $ LDA ) 631* 632 END IF 633* 634 70 CONTINUE 635* 636* Post-assignment 637* 638 DO 80 I = M + 1, K + L 639 ALPHA( I ) = ZERO 640 BETA( I ) = ONE 641 80 CONTINUE 642* 643 IF( K+L.LT.N ) THEN 644 DO 90 I = K + L + 1, N 645 ALPHA( I ) = ZERO 646 BETA( I ) = ZERO 647 90 CONTINUE 648 END IF 649* 650 100 CONTINUE 651 NCYCLE = KCYCLE 652 RETURN 653* 654* End of STGSJA 655* 656 END 657