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*> \date November 2011 355* 356*> \ingroup complex16OTHERcomputational 357* 358*> \par Further Details: 359* ===================== 360*> 361*> \verbatim 362*> 363*> ZTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce 364*> min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L 365*> matrix B13 to the form: 366*> 367*> U1**H *A13*Q1 = C1*R1; V1**H *B13*Q1 = S1*R1, 368*> 369*> where U1, V1 and Q1 are unitary matrix. 370*> C1 and S1 are diagonal matrices satisfying 371*> 372*> C1**2 + S1**2 = I, 373*> 374*> and R1 is an L-by-L nonsingular upper triangular matrix. 375*> \endverbatim 376*> 377* ===================================================================== 378 SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, 379 $ LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, 380 $ Q, LDQ, WORK, NCYCLE, INFO ) 381* 382* -- LAPACK computational routine (version 3.4.0) -- 383* -- LAPACK is a software package provided by Univ. of Tennessee, -- 384* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 385* November 2011 386* 387* .. Scalar Arguments .. 388 CHARACTER JOBQ, JOBU, JOBV 389 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, 390 $ NCYCLE, P 391 DOUBLE PRECISION TOLA, TOLB 392* .. 393* .. Array Arguments .. 394 DOUBLE PRECISION ALPHA( * ), BETA( * ) 395 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 396 $ U( LDU, * ), V( LDV, * ), WORK( * ) 397* .. 398* 399* ===================================================================== 400* 401* .. Parameters .. 402 INTEGER MAXIT 403 PARAMETER ( MAXIT = 40 ) 404 DOUBLE PRECISION ZERO, ONE 405 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 406 COMPLEX*16 CZERO, CONE 407 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 408 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 409* .. 410* .. Local Scalars .. 411* 412 LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV 413 INTEGER I, J, KCYCLE 414 DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV, ERROR, GAMMA, 415 $ RWK, SSMIN 416 COMPLEX*16 A2, B2, SNQ, SNU, SNV 417* .. 418* .. External Functions .. 419 LOGICAL LSAME 420 EXTERNAL LSAME 421* .. 422* .. External Subroutines .. 423 EXTERNAL DLARTG, XERBLA, ZCOPY, ZDSCAL, ZLAGS2, ZLAPLL, 424 $ ZLASET, ZROT 425* .. 426* .. Intrinsic Functions .. 427 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN 428* .. 429* .. Executable Statements .. 430* 431* Decode and test the input parameters 432* 433 INITU = LSAME( JOBU, 'I' ) 434 WANTU = INITU .OR. LSAME( JOBU, 'U' ) 435* 436 INITV = LSAME( JOBV, 'I' ) 437 WANTV = INITV .OR. LSAME( JOBV, 'V' ) 438* 439 INITQ = LSAME( JOBQ, 'I' ) 440 WANTQ = INITQ .OR. LSAME( JOBQ, 'Q' ) 441* 442 INFO = 0 443 IF( .NOT.( INITU .OR. WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 444 INFO = -1 445 ELSE IF( .NOT.( INITV .OR. WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 446 INFO = -2 447 ELSE IF( .NOT.( INITQ .OR. WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 448 INFO = -3 449 ELSE IF( M.LT.0 ) THEN 450 INFO = -4 451 ELSE IF( P.LT.0 ) THEN 452 INFO = -5 453 ELSE IF( N.LT.0 ) THEN 454 INFO = -6 455 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 456 INFO = -10 457 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 458 INFO = -12 459 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 460 INFO = -18 461 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 462 INFO = -20 463 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 464 INFO = -22 465 END IF 466 IF( INFO.NE.0 ) THEN 467 CALL XERBLA( 'ZTGSJA', -INFO ) 468 RETURN 469 END IF 470* 471* Initialize U, V and Q, if necessary 472* 473 IF( INITU ) 474 $ CALL ZLASET( 'Full', M, M, CZERO, CONE, U, LDU ) 475 IF( INITV ) 476 $ CALL ZLASET( 'Full', P, P, CZERO, CONE, V, LDV ) 477 IF( INITQ ) 478 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 479* 480* Loop until convergence 481* 482 UPPER = .FALSE. 483 DO 40 KCYCLE = 1, MAXIT 484* 485 UPPER = .NOT.UPPER 486* 487 DO 20 I = 1, L - 1 488 DO 10 J = I + 1, L 489* 490 A1 = ZERO 491 A2 = CZERO 492 A3 = ZERO 493 IF( K+I.LE.M ) 494 $ A1 = DBLE( A( K+I, N-L+I ) ) 495 IF( K+J.LE.M ) 496 $ A3 = DBLE( A( K+J, N-L+J ) ) 497* 498 B1 = DBLE( B( I, N-L+I ) ) 499 B3 = DBLE( B( J, N-L+J ) ) 500* 501 IF( UPPER ) THEN 502 IF( K+I.LE.M ) 503 $ A2 = A( K+I, N-L+J ) 504 B2 = B( I, N-L+J ) 505 ELSE 506 IF( K+J.LE.M ) 507 $ A2 = A( K+J, N-L+I ) 508 B2 = B( J, N-L+I ) 509 END IF 510* 511 CALL ZLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, 512 $ CSV, SNV, CSQ, SNQ ) 513* 514* Update (K+I)-th and (K+J)-th rows of matrix A: U**H *A 515* 516 IF( K+J.LE.M ) 517 $ CALL ZROT( L, A( K+J, N-L+1 ), LDA, A( K+I, N-L+1 ), 518 $ LDA, CSU, DCONJG( SNU ) ) 519* 520* Update I-th and J-th rows of matrix B: V**H *B 521* 522 CALL ZROT( L, B( J, N-L+1 ), LDB, B( I, N-L+1 ), LDB, 523 $ CSV, DCONJG( SNV ) ) 524* 525* Update (N-L+I)-th and (N-L+J)-th columns of matrices 526* A and B: A*Q and B*Q 527* 528 CALL ZROT( MIN( K+L, M ), A( 1, N-L+J ), 1, 529 $ A( 1, N-L+I ), 1, CSQ, SNQ ) 530* 531 CALL ZROT( L, B( 1, N-L+J ), 1, B( 1, N-L+I ), 1, CSQ, 532 $ SNQ ) 533* 534 IF( UPPER ) THEN 535 IF( K+I.LE.M ) 536 $ A( K+I, N-L+J ) = CZERO 537 B( I, N-L+J ) = CZERO 538 ELSE 539 IF( K+J.LE.M ) 540 $ A( K+J, N-L+I ) = CZERO 541 B( J, N-L+I ) = CZERO 542 END IF 543* 544* Ensure that the diagonal elements of A and B are real. 545* 546 IF( K+I.LE.M ) 547 $ A( K+I, N-L+I ) = DBLE( A( K+I, N-L+I ) ) 548 IF( K+J.LE.M ) 549 $ A( K+J, N-L+J ) = DBLE( A( K+J, N-L+J ) ) 550 B( I, N-L+I ) = DBLE( B( I, N-L+I ) ) 551 B( J, N-L+J ) = DBLE( B( J, N-L+J ) ) 552* 553* Update unitary matrices U, V, Q, if desired. 554* 555 IF( WANTU .AND. K+J.LE.M ) 556 $ CALL ZROT( M, U( 1, K+J ), 1, U( 1, K+I ), 1, CSU, 557 $ SNU ) 558* 559 IF( WANTV ) 560 $ CALL ZROT( P, V( 1, J ), 1, V( 1, I ), 1, CSV, SNV ) 561* 562 IF( WANTQ ) 563 $ CALL ZROT( N, Q( 1, N-L+J ), 1, Q( 1, N-L+I ), 1, CSQ, 564 $ SNQ ) 565* 566 10 CONTINUE 567 20 CONTINUE 568* 569 IF( .NOT.UPPER ) THEN 570* 571* The matrices A13 and B13 were lower triangular at the start 572* of the cycle, and are now upper triangular. 573* 574* Convergence test: test the parallelism of the corresponding 575* rows of A and B. 576* 577 ERROR = ZERO 578 DO 30 I = 1, MIN( L, M-K ) 579 CALL ZCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 ) 580 CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 ) 581 CALL ZLAPLL( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN ) 582 ERROR = MAX( ERROR, SSMIN ) 583 30 CONTINUE 584* 585 IF( ABS( ERROR ).LE.MIN( TOLA, TOLB ) ) 586 $ GO TO 50 587 END IF 588* 589* End of cycle loop 590* 591 40 CONTINUE 592* 593* The algorithm has not converged after MAXIT cycles. 594* 595 INFO = 1 596 GO TO 100 597* 598 50 CONTINUE 599* 600* If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged. 601* Compute the generalized singular value pairs (ALPHA, BETA), and 602* set the triangular matrix R to array A. 603* 604 DO 60 I = 1, K 605 ALPHA( I ) = ONE 606 BETA( I ) = ZERO 607 60 CONTINUE 608* 609 DO 70 I = 1, MIN( L, M-K ) 610* 611 A1 = DBLE( A( K+I, N-L+I ) ) 612 B1 = DBLE( B( I, N-L+I ) ) 613* 614 IF( A1.NE.ZERO ) THEN 615 GAMMA = B1 / A1 616* 617 IF( GAMMA.LT.ZERO ) THEN 618 CALL ZDSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) 619 IF( WANTV ) 620 $ CALL ZDSCAL( P, -ONE, V( 1, I ), 1 ) 621 END IF 622* 623 CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), 624 $ RWK ) 625* 626 IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN 627 CALL ZDSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), 628 $ LDA ) 629 ELSE 630 CALL ZDSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), 631 $ LDB ) 632 CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 633 $ LDA ) 634 END IF 635* 636 ELSE 637* 638 ALPHA( K+I ) = ZERO 639 BETA( K+I ) = ONE 640 CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 641 $ LDA ) 642 END IF 643 70 CONTINUE 644* 645* Post-assignment 646* 647 DO 80 I = M + 1, K + L 648 ALPHA( I ) = ZERO 649 BETA( I ) = ONE 650 80 CONTINUE 651* 652 IF( K+L.LT.N ) THEN 653 DO 90 I = K + L + 1, N 654 ALPHA( I ) = ZERO 655 BETA( I ) = ZERO 656 90 CONTINUE 657 END IF 658* 659 100 CONTINUE 660 NCYCLE = KCYCLE 661* 662 RETURN 663* 664* End of ZTGSJA 665* 666 END 667