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