1*> \brief \b ZTGSYL 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZTGSYL + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsyl.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsyl.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsyl.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 22* LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, 23* IWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER TRANS 27* INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, 28* $ LWORK, M, N 29* DOUBLE PRECISION DIF, SCALE 30* .. 31* .. Array Arguments .. 32* INTEGER IWORK( * ) 33* COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ), 34* $ D( LDD, * ), E( LDE, * ), F( LDF, * ), 35* $ WORK( * ) 36* .. 37* 38* 39*> \par Purpose: 40* ============= 41*> 42*> \verbatim 43*> 44*> ZTGSYL solves the generalized Sylvester equation: 45*> 46*> A * R - L * B = scale * C (1) 47*> D * R - L * E = scale * F 48*> 49*> where R and L are unknown m-by-n matrices, (A, D), (B, E) and 50*> (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n, 51*> respectively, with complex entries. A, B, D and E are upper 52*> triangular (i.e., (A,D) and (B,E) in generalized Schur form). 53*> 54*> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 55*> is an output scaling factor chosen to avoid overflow. 56*> 57*> In matrix notation (1) is equivalent to solve Zx = scale*b, where Z 58*> is defined as 59*> 60*> Z = [ kron(In, A) -kron(B**H, Im) ] (2) 61*> [ kron(In, D) -kron(E**H, Im) ], 62*> 63*> Here Ix is the identity matrix of size x and X**H is the conjugate 64*> transpose of X. Kron(X, Y) is the Kronecker product between the 65*> matrices X and Y. 66*> 67*> If TRANS = 'C', y in the conjugate transposed system Z**H *y = scale*b 68*> is solved for, which is equivalent to solve for R and L in 69*> 70*> A**H * R + D**H * L = scale * C (3) 71*> R * B**H + L * E**H = scale * -F 72*> 73*> This case (TRANS = 'C') is used to compute an one-norm-based estimate 74*> of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D) 75*> and (B,E), using ZLACON. 76*> 77*> If IJOB >= 1, ZTGSYL computes a Frobenius norm-based estimate of 78*> Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the 79*> reciprocal of the smallest singular value of Z. 80*> 81*> This is a level-3 BLAS algorithm. 82*> \endverbatim 83* 84* Arguments: 85* ========== 86* 87*> \param[in] TRANS 88*> \verbatim 89*> TRANS is CHARACTER*1 90*> = 'N': solve the generalized sylvester equation (1). 91*> = 'C': solve the "conjugate transposed" system (3). 92*> \endverbatim 93*> 94*> \param[in] IJOB 95*> \verbatim 96*> IJOB is INTEGER 97*> Specifies what kind of functionality to be performed. 98*> =0: solve (1) only. 99*> =1: The functionality of 0 and 3. 100*> =2: The functionality of 0 and 4. 101*> =3: Only an estimate of Dif[(A,D), (B,E)] is computed. 102*> (look ahead strategy is used). 103*> =4: Only an estimate of Dif[(A,D), (B,E)] is computed. 104*> (ZGECON on sub-systems is used). 105*> Not referenced if TRANS = 'C'. 106*> \endverbatim 107*> 108*> \param[in] M 109*> \verbatim 110*> M is INTEGER 111*> The order of the matrices A and D, and the row dimension of 112*> the matrices C, F, R and L. 113*> \endverbatim 114*> 115*> \param[in] N 116*> \verbatim 117*> N is INTEGER 118*> The order of the matrices B and E, and the column dimension 119*> of the matrices C, F, R and L. 120*> \endverbatim 121*> 122*> \param[in] A 123*> \verbatim 124*> A is COMPLEX*16 array, dimension (LDA, M) 125*> The upper triangular matrix A. 126*> \endverbatim 127*> 128*> \param[in] LDA 129*> \verbatim 130*> LDA is INTEGER 131*> The leading dimension of the array A. LDA >= max(1, M). 132*> \endverbatim 133*> 134*> \param[in] B 135*> \verbatim 136*> B is COMPLEX*16 array, dimension (LDB, N) 137*> The upper triangular matrix B. 138*> \endverbatim 139*> 140*> \param[in] LDB 141*> \verbatim 142*> LDB is INTEGER 143*> The leading dimension of the array B. LDB >= max(1, N). 144*> \endverbatim 145*> 146*> \param[in,out] C 147*> \verbatim 148*> C is COMPLEX*16 array, dimension (LDC, N) 149*> On entry, C contains the right-hand-side of the first matrix 150*> equation in (1) or (3). 151*> On exit, if IJOB = 0, 1 or 2, C has been overwritten by 152*> the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R, 153*> the solution achieved during the computation of the 154*> Dif-estimate. 155*> \endverbatim 156*> 157*> \param[in] LDC 158*> \verbatim 159*> LDC is INTEGER 160*> The leading dimension of the array C. LDC >= max(1, M). 161*> \endverbatim 162*> 163*> \param[in] D 164*> \verbatim 165*> D is COMPLEX*16 array, dimension (LDD, M) 166*> The upper triangular matrix D. 167*> \endverbatim 168*> 169*> \param[in] LDD 170*> \verbatim 171*> LDD is INTEGER 172*> The leading dimension of the array D. LDD >= max(1, M). 173*> \endverbatim 174*> 175*> \param[in] E 176*> \verbatim 177*> E is COMPLEX*16 array, dimension (LDE, N) 178*> The upper triangular matrix E. 179*> \endverbatim 180*> 181*> \param[in] LDE 182*> \verbatim 183*> LDE is INTEGER 184*> The leading dimension of the array E. LDE >= max(1, N). 185*> \endverbatim 186*> 187*> \param[in,out] F 188*> \verbatim 189*> F is COMPLEX*16 array, dimension (LDF, N) 190*> On entry, F contains the right-hand-side of the second matrix 191*> equation in (1) or (3). 192*> On exit, if IJOB = 0, 1 or 2, F has been overwritten by 193*> the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L, 194*> the solution achieved during the computation of the 195*> Dif-estimate. 196*> \endverbatim 197*> 198*> \param[in] LDF 199*> \verbatim 200*> LDF is INTEGER 201*> The leading dimension of the array F. LDF >= max(1, M). 202*> \endverbatim 203*> 204*> \param[out] DIF 205*> \verbatim 206*> DIF is DOUBLE PRECISION 207*> On exit DIF is the reciprocal of a lower bound of the 208*> reciprocal of the Dif-function, i.e. DIF is an upper bound of 209*> Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2). 210*> IF IJOB = 0 or TRANS = 'C', DIF is not referenced. 211*> \endverbatim 212*> 213*> \param[out] SCALE 214*> \verbatim 215*> SCALE is DOUBLE PRECISION 216*> On exit SCALE is the scaling factor in (1) or (3). 217*> If 0 < SCALE < 1, C and F hold the solutions R and L, resp., 218*> to a slightly perturbed system but the input matrices A, B, 219*> D and E have not been changed. If SCALE = 0, R and L will 220*> hold the solutions to the homogenious system with C = F = 0. 221*> \endverbatim 222*> 223*> \param[out] WORK 224*> \verbatim 225*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 226*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 227*> \endverbatim 228*> 229*> \param[in] LWORK 230*> \verbatim 231*> LWORK is INTEGER 232*> The dimension of the array WORK. LWORK > = 1. 233*> If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N). 234*> 235*> If LWORK = -1, then a workspace query is assumed; the routine 236*> only calculates the optimal size of the WORK array, returns 237*> this value as the first entry of the WORK array, and no error 238*> message related to LWORK is issued by XERBLA. 239*> \endverbatim 240*> 241*> \param[out] IWORK 242*> \verbatim 243*> IWORK is INTEGER array, dimension (M+N+2) 244*> \endverbatim 245*> 246*> \param[out] INFO 247*> \verbatim 248*> INFO is INTEGER 249*> =0: successful exit 250*> <0: If INFO = -i, the i-th argument had an illegal value. 251*> >0: (A, D) and (B, E) have common or very close 252*> eigenvalues. 253*> \endverbatim 254* 255* Authors: 256* ======== 257* 258*> \author Univ. of Tennessee 259*> \author Univ. of California Berkeley 260*> \author Univ. of Colorado Denver 261*> \author NAG Ltd. 262* 263*> \date November 2011 264* 265*> \ingroup complex16SYcomputational 266* 267*> \par Contributors: 268* ================== 269*> 270*> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 271*> Umea University, S-901 87 Umea, Sweden. 272* 273*> \par References: 274* ================ 275*> 276*> [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 277*> for Solving the Generalized Sylvester Equation and Estimating the 278*> Separation between Regular Matrix Pairs, Report UMINF - 93.23, 279*> Department of Computing Science, Umea University, S-901 87 Umea, 280*> Sweden, December 1993, Revised April 1994, Also as LAPACK Working 281*> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, 282*> No 1, 1996. 283*> \n 284*> [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester 285*> Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal. 286*> Appl., 15(4):1045-1060, 1994. 287*> \n 288*> [3] B. Kagstrom and L. Westin, Generalized Schur Methods with 289*> Condition Estimators for Solving the Generalized Sylvester 290*> Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7, 291*> July 1989, pp 745-751. 292*> 293* ===================================================================== 294 SUBROUTINE ZTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 295 $ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, 296 $ IWORK, INFO ) 297* 298* -- LAPACK computational routine (version 3.4.0) -- 299* -- LAPACK is a software package provided by Univ. of Tennessee, -- 300* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 301* November 2011 302* 303* .. Scalar Arguments .. 304 CHARACTER TRANS 305 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, 306 $ LWORK, M, N 307 DOUBLE PRECISION DIF, SCALE 308* .. 309* .. Array Arguments .. 310 INTEGER IWORK( * ) 311 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ), 312 $ D( LDD, * ), E( LDE, * ), F( LDF, * ), 313 $ WORK( * ) 314* .. 315* 316* ===================================================================== 317* Replaced various illegal calls to CCOPY by calls to CLASET. 318* Sven Hammarling, 1/5/02. 319* 320* .. Parameters .. 321 DOUBLE PRECISION ZERO, ONE 322 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 323 COMPLEX*16 CZERO 324 PARAMETER ( CZERO = (0.0D+0, 0.0D+0) ) 325* .. 326* .. Local Scalars .. 327 LOGICAL LQUERY, NOTRAN 328 INTEGER I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K, 329 $ LINFO, LWMIN, MB, NB, P, PQ, Q 330 DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC 331* .. 332* .. External Functions .. 333 LOGICAL LSAME 334 INTEGER ILAENV 335 EXTERNAL LSAME, ILAENV 336* .. 337* .. External Subroutines .. 338 EXTERNAL XERBLA, ZGEMM, ZLACPY, ZLASET, ZSCAL, ZTGSY2 339* .. 340* .. Intrinsic Functions .. 341 INTRINSIC DBLE, DCMPLX, MAX, SQRT 342* .. 343* .. Executable Statements .. 344* 345* Decode and test input parameters 346* 347 INFO = 0 348 NOTRAN = LSAME( TRANS, 'N' ) 349 LQUERY = ( LWORK.EQ.-1 ) 350* 351 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 352 INFO = -1 353 ELSE IF( NOTRAN ) THEN 354 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN 355 INFO = -2 356 END IF 357 END IF 358 IF( INFO.EQ.0 ) THEN 359 IF( M.LE.0 ) THEN 360 INFO = -3 361 ELSE IF( N.LE.0 ) THEN 362 INFO = -4 363 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 364 INFO = -6 365 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 366 INFO = -8 367 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 368 INFO = -10 369 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 370 INFO = -12 371 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 372 INFO = -14 373 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 374 INFO = -16 375 END IF 376 END IF 377* 378 IF( INFO.EQ.0 ) THEN 379 IF( NOTRAN ) THEN 380 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN 381 LWMIN = MAX( 1, 2*M*N ) 382 ELSE 383 LWMIN = 1 384 END IF 385 ELSE 386 LWMIN = 1 387 END IF 388 WORK( 1 ) = LWMIN 389* 390 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 391 INFO = -20 392 END IF 393 END IF 394* 395 IF( INFO.NE.0 ) THEN 396 CALL XERBLA( 'ZTGSYL', -INFO ) 397 RETURN 398 ELSE IF( LQUERY ) THEN 399 RETURN 400 END IF 401* 402* Quick return if possible 403* 404 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 405 SCALE = 1 406 IF( NOTRAN ) THEN 407 IF( IJOB.NE.0 ) THEN 408 DIF = 0 409 END IF 410 END IF 411 RETURN 412 END IF 413* 414* Determine optimal block sizes MB and NB 415* 416 MB = ILAENV( 2, 'ZTGSYL', TRANS, M, N, -1, -1 ) 417 NB = ILAENV( 5, 'ZTGSYL', TRANS, M, N, -1, -1 ) 418* 419 ISOLVE = 1 420 IFUNC = 0 421 IF( NOTRAN ) THEN 422 IF( IJOB.GE.3 ) THEN 423 IFUNC = IJOB - 2 424 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 425 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF ) 426 ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN 427 ISOLVE = 2 428 END IF 429 END IF 430* 431 IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) ) 432 $ THEN 433* 434* Use unblocked Level 2 solver 435* 436 DO 30 IROUND = 1, ISOLVE 437* 438 SCALE = ONE 439 DSCALE = ZERO 440 DSUM = ONE 441 PQ = M*N 442 CALL ZTGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D, 443 $ LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE, 444 $ INFO ) 445 IF( DSCALE.NE.ZERO ) THEN 446 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN 447 DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) ) 448 ELSE 449 DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) ) 450 END IF 451 END IF 452 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN 453 IF( NOTRAN ) THEN 454 IFUNC = IJOB 455 END IF 456 SCALE2 = SCALE 457 CALL ZLACPY( 'F', M, N, C, LDC, WORK, M ) 458 CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M ) 459 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 460 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF ) 461 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN 462 CALL ZLACPY( 'F', M, N, WORK, M, C, LDC ) 463 CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF ) 464 SCALE = SCALE2 465 END IF 466 30 CONTINUE 467* 468 RETURN 469* 470 END IF 471* 472* Determine block structure of A 473* 474 P = 0 475 I = 1 476 40 CONTINUE 477 IF( I.GT.M ) 478 $ GO TO 50 479 P = P + 1 480 IWORK( P ) = I 481 I = I + MB 482 IF( I.GE.M ) 483 $ GO TO 50 484 GO TO 40 485 50 CONTINUE 486 IWORK( P+1 ) = M + 1 487 IF( IWORK( P ).EQ.IWORK( P+1 ) ) 488 $ P = P - 1 489* 490* Determine block structure of B 491* 492 Q = P + 1 493 J = 1 494 60 CONTINUE 495 IF( J.GT.N ) 496 $ GO TO 70 497* 498 Q = Q + 1 499 IWORK( Q ) = J 500 J = J + NB 501 IF( J.GE.N ) 502 $ GO TO 70 503 GO TO 60 504* 505 70 CONTINUE 506 IWORK( Q+1 ) = N + 1 507 IF( IWORK( Q ).EQ.IWORK( Q+1 ) ) 508 $ Q = Q - 1 509* 510 IF( NOTRAN ) THEN 511 DO 150 IROUND = 1, ISOLVE 512* 513* Solve (I, J) - subsystem 514* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) 515* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) 516* for I = P, P - 1, ..., 1; J = 1, 2, ..., Q 517* 518 PQ = 0 519 SCALE = ONE 520 DSCALE = ZERO 521 DSUM = ONE 522 DO 130 J = P + 2, Q 523 JS = IWORK( J ) 524 JE = IWORK( J+1 ) - 1 525 NB = JE - JS + 1 526 DO 120 I = P, 1, -1 527 IS = IWORK( I ) 528 IE = IWORK( I+1 ) - 1 529 MB = IE - IS + 1 530 CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA, 531 $ B( JS, JS ), LDB, C( IS, JS ), LDC, 532 $ D( IS, IS ), LDD, E( JS, JS ), LDE, 533 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE, 534 $ LINFO ) 535 IF( LINFO.GT.0 ) 536 $ INFO = LINFO 537 PQ = PQ + MB*NB 538 IF( SCALOC.NE.ONE ) THEN 539 DO 80 K = 1, JS - 1 540 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 541 $ C( 1, K ), 1 ) 542 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 543 $ F( 1, K ), 1 ) 544 80 CONTINUE 545 DO 90 K = JS, JE 546 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ), 547 $ C( 1, K ), 1 ) 548 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ), 549 $ F( 1, K ), 1 ) 550 90 CONTINUE 551 DO 100 K = JS, JE 552 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ), 553 $ C( IE+1, K ), 1 ) 554 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ), 555 $ F( IE+1, K ), 1 ) 556 100 CONTINUE 557 DO 110 K = JE + 1, N 558 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 559 $ C( 1, K ), 1 ) 560 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 561 $ F( 1, K ), 1 ) 562 110 CONTINUE 563 SCALE = SCALE*SCALOC 564 END IF 565* 566* Substitute R(I,J) and L(I,J) into remaining equation. 567* 568 IF( I.GT.1 ) THEN 569 CALL ZGEMM( 'N', 'N', IS-1, NB, MB, 570 $ DCMPLX( -ONE, ZERO ), A( 1, IS ), LDA, 571 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ), 572 $ C( 1, JS ), LDC ) 573 CALL ZGEMM( 'N', 'N', IS-1, NB, MB, 574 $ DCMPLX( -ONE, ZERO ), D( 1, IS ), LDD, 575 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ), 576 $ F( 1, JS ), LDF ) 577 END IF 578 IF( J.LT.Q ) THEN 579 CALL ZGEMM( 'N', 'N', MB, N-JE, NB, 580 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF, 581 $ B( JS, JE+1 ), LDB, 582 $ DCMPLX( ONE, ZERO ), C( IS, JE+1 ), 583 $ LDC ) 584 CALL ZGEMM( 'N', 'N', MB, N-JE, NB, 585 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF, 586 $ E( JS, JE+1 ), LDE, 587 $ DCMPLX( ONE, ZERO ), F( IS, JE+1 ), 588 $ LDF ) 589 END IF 590 120 CONTINUE 591 130 CONTINUE 592 IF( DSCALE.NE.ZERO ) THEN 593 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN 594 DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) ) 595 ELSE 596 DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) ) 597 END IF 598 END IF 599 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN 600 IF( NOTRAN ) THEN 601 IFUNC = IJOB 602 END IF 603 SCALE2 = SCALE 604 CALL ZLACPY( 'F', M, N, C, LDC, WORK, M ) 605 CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M ) 606 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 607 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF ) 608 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN 609 CALL ZLACPY( 'F', M, N, WORK, M, C, LDC ) 610 CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF ) 611 SCALE = SCALE2 612 END IF 613 150 CONTINUE 614 ELSE 615* 616* Solve transposed (I, J)-subsystem 617* A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J) 618* R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 619* for I = 1,2,..., P; J = Q, Q-1,..., 1 620* 621 SCALE = ONE 622 DO 210 I = 1, P 623 IS = IWORK( I ) 624 IE = IWORK( I+1 ) - 1 625 MB = IE - IS + 1 626 DO 200 J = Q, P + 2, -1 627 JS = IWORK( J ) 628 JE = IWORK( J+1 ) - 1 629 NB = JE - JS + 1 630 CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA, 631 $ B( JS, JS ), LDB, C( IS, JS ), LDC, 632 $ D( IS, IS ), LDD, E( JS, JS ), LDE, 633 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE, 634 $ LINFO ) 635 IF( LINFO.GT.0 ) 636 $ INFO = LINFO 637 IF( SCALOC.NE.ONE ) THEN 638 DO 160 K = 1, JS - 1 639 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ), 640 $ 1 ) 641 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ), 642 $ 1 ) 643 160 CONTINUE 644 DO 170 K = JS, JE 645 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ), 646 $ C( 1, K ), 1 ) 647 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ), 648 $ F( 1, K ), 1 ) 649 170 CONTINUE 650 DO 180 K = JS, JE 651 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ), 652 $ C( IE+1, K ), 1 ) 653 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ), 654 $ F( IE+1, K ), 1 ) 655 180 CONTINUE 656 DO 190 K = JE + 1, N 657 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ), 658 $ 1 ) 659 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ), 660 $ 1 ) 661 190 CONTINUE 662 SCALE = SCALE*SCALOC 663 END IF 664* 665* Substitute R(I,J) and L(I,J) into remaining equation. 666* 667 IF( J.GT.P+2 ) THEN 668 CALL ZGEMM( 'N', 'C', MB, JS-1, NB, 669 $ DCMPLX( ONE, ZERO ), C( IS, JS ), LDC, 670 $ B( 1, JS ), LDB, DCMPLX( ONE, ZERO ), 671 $ F( IS, 1 ), LDF ) 672 CALL ZGEMM( 'N', 'C', MB, JS-1, NB, 673 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF, 674 $ E( 1, JS ), LDE, DCMPLX( ONE, ZERO ), 675 $ F( IS, 1 ), LDF ) 676 END IF 677 IF( I.LT.P ) THEN 678 CALL ZGEMM( 'C', 'N', M-IE, NB, MB, 679 $ DCMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA, 680 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ), 681 $ C( IE+1, JS ), LDC ) 682 CALL ZGEMM( 'C', 'N', M-IE, NB, MB, 683 $ DCMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD, 684 $ F( IS, JS ), LDF, DCMPLX( ONE, ZERO ), 685 $ C( IE+1, JS ), LDC ) 686 END IF 687 200 CONTINUE 688 210 CONTINUE 689 END IF 690* 691 WORK( 1 ) = LWMIN 692* 693 RETURN 694* 695* End of ZTGSYL 696* 697 END 698