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