1*> \brief \b STGSY2 solves the generalized Sylvester equation (unblocked algorithm). 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download STGSY2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsy2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsy2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsy2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE STGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 22* LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, 23* IWORK, PQ, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER TRANS 27* INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N, 28* $ PQ 29* REAL RDSCAL, RDSUM, SCALE 30* .. 31* .. Array Arguments .. 32* INTEGER IWORK( * ) 33* REAL A( LDA, * ), B( LDB, * ), C( LDC, * ), 34* $ D( LDD, * ), E( LDE, * ), F( LDF, * ) 35* .. 36* 37* 38*> \par Purpose: 39* ============= 40*> 41*> \verbatim 42*> 43*> STGSY2 solves the generalized Sylvester equation: 44*> 45*> A * R - L * B = scale * C (1) 46*> D * R - L * E = scale * F, 47*> 48*> using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices, 49*> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M, 50*> N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E) 51*> must be in generalized Schur canonical form, i.e. A, B are upper 52*> quasi triangular and D, E are upper triangular. The solution (R, L) 53*> overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor 54*> chosen to avoid overflow. 55*> 56*> In matrix notation solving equation (1) corresponds to solve 57*> Z*x = scale*b, where Z is defined as 58*> 59*> Z = [ kron(In, A) -kron(B**T, Im) ] (2) 60*> [ kron(In, D) -kron(E**T, Im) ], 61*> 62*> Ik is the identity matrix of size k and X**T is the transpose of X. 63*> kron(X, Y) is the Kronecker product between the matrices X and Y. 64*> In the process of solving (1), we solve a number of such systems 65*> where Dim(In), Dim(In) = 1 or 2. 66*> 67*> If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y, 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 is used to compute an estimate of Dif[(A, D), (B, E)] = 74*> sigma_min(Z) using reverse communication with SLACON. 75*> 76*> STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL 77*> of an upper bound on the separation between to matrix pairs. Then 78*> the input (A, D), (B, E) are sub-pencils of the matrix pair in 79*> STGSYL. See STGSYL for details. 80*> \endverbatim 81* 82* Arguments: 83* ========== 84* 85*> \param[in] TRANS 86*> \verbatim 87*> TRANS is CHARACTER*1 88*> = 'N': solve the generalized Sylvester equation (1). 89*> = 'T': solve the 'transposed' system (3). 90*> \endverbatim 91*> 92*> \param[in] IJOB 93*> \verbatim 94*> IJOB is INTEGER 95*> Specifies what kind of functionality to be performed. 96*> = 0: solve (1) only. 97*> = 1: A contribution from this subsystem to a Frobenius 98*> norm-based estimate of the separation between two matrix 99*> pairs is computed. (look ahead strategy is used). 100*> = 2: A contribution from this subsystem to a Frobenius 101*> norm-based estimate of the separation between two matrix 102*> pairs is computed. (SGECON on sub-systems is used.) 103*> Not referenced if TRANS = 'T'. 104*> \endverbatim 105*> 106*> \param[in] M 107*> \verbatim 108*> M is INTEGER 109*> On entry, M specifies the order of A and D, and the row 110*> dimension of C, F, R and L. 111*> \endverbatim 112*> 113*> \param[in] N 114*> \verbatim 115*> N is INTEGER 116*> On entry, N specifies the order of B and E, and the column 117*> dimension of C, F, R and L. 118*> \endverbatim 119*> 120*> \param[in] A 121*> \verbatim 122*> A is REAL array, dimension (LDA, M) 123*> On entry, A contains an upper quasi triangular matrix. 124*> \endverbatim 125*> 126*> \param[in] LDA 127*> \verbatim 128*> LDA is INTEGER 129*> The leading dimension of the matrix A. LDA >= max(1, M). 130*> \endverbatim 131*> 132*> \param[in] B 133*> \verbatim 134*> B is REAL array, dimension (LDB, N) 135*> On entry, B contains an upper quasi triangular matrix. 136*> \endverbatim 137*> 138*> \param[in] LDB 139*> \verbatim 140*> LDB is INTEGER 141*> The leading dimension of the matrix B. LDB >= max(1, N). 142*> \endverbatim 143*> 144*> \param[in,out] C 145*> \verbatim 146*> C is REAL array, dimension (LDC, N) 147*> On entry, C contains the right-hand-side of the first matrix 148*> equation in (1). 149*> On exit, if IJOB = 0, C has been overwritten by the 150*> solution R. 151*> \endverbatim 152*> 153*> \param[in] LDC 154*> \verbatim 155*> LDC is INTEGER 156*> The leading dimension of the matrix C. LDC >= max(1, M). 157*> \endverbatim 158*> 159*> \param[in] D 160*> \verbatim 161*> D is REAL array, dimension (LDD, M) 162*> On entry, D contains an upper triangular matrix. 163*> \endverbatim 164*> 165*> \param[in] LDD 166*> \verbatim 167*> LDD is INTEGER 168*> The leading dimension of the matrix D. LDD >= max(1, M). 169*> \endverbatim 170*> 171*> \param[in] E 172*> \verbatim 173*> E is REAL array, dimension (LDE, N) 174*> On entry, E contains an upper triangular matrix. 175*> \endverbatim 176*> 177*> \param[in] LDE 178*> \verbatim 179*> LDE is INTEGER 180*> The leading dimension of the matrix E. LDE >= max(1, N). 181*> \endverbatim 182*> 183*> \param[in,out] F 184*> \verbatim 185*> F is REAL array, dimension (LDF, N) 186*> On entry, F contains the right-hand-side of the second matrix 187*> equation in (1). 188*> On exit, if IJOB = 0, F has been overwritten by the 189*> solution L. 190*> \endverbatim 191*> 192*> \param[in] LDF 193*> \verbatim 194*> LDF is INTEGER 195*> The leading dimension of the matrix F. LDF >= max(1, M). 196*> \endverbatim 197*> 198*> \param[out] SCALE 199*> \verbatim 200*> SCALE is REAL 201*> On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions 202*> R and L (C and F on entry) will hold the solutions to a 203*> slightly perturbed system but the input matrices A, B, D and 204*> E have not been changed. If SCALE = 0, R and L will hold the 205*> solutions to the homogeneous system with C = F = 0. Normally, 206*> SCALE = 1. 207*> \endverbatim 208*> 209*> \param[in,out] RDSUM 210*> \verbatim 211*> RDSUM is REAL 212*> On entry, the sum of squares of computed contributions to 213*> the Dif-estimate under computation by STGSYL, where the 214*> scaling factor RDSCAL (see below) has been factored out. 215*> On exit, the corresponding sum of squares updated with the 216*> contributions from the current sub-system. 217*> If TRANS = 'T' RDSUM is not touched. 218*> NOTE: RDSUM only makes sense when STGSY2 is called by STGSYL. 219*> \endverbatim 220*> 221*> \param[in,out] RDSCAL 222*> \verbatim 223*> RDSCAL is REAL 224*> On entry, scaling factor used to prevent overflow in RDSUM. 225*> On exit, RDSCAL is updated w.r.t. the current contributions 226*> in RDSUM. 227*> If TRANS = 'T', RDSCAL is not touched. 228*> NOTE: RDSCAL only makes sense when STGSY2 is called by 229*> STGSYL. 230*> \endverbatim 231*> 232*> \param[out] IWORK 233*> \verbatim 234*> IWORK is INTEGER array, dimension (M+N+2) 235*> \endverbatim 236*> 237*> \param[out] PQ 238*> \verbatim 239*> PQ is INTEGER 240*> On exit, the number of subsystems (of size 2-by-2, 4-by-4 and 241*> 8-by-8) solved by this routine. 242*> \endverbatim 243*> 244*> \param[out] INFO 245*> \verbatim 246*> INFO is INTEGER 247*> On exit, if INFO is set to 248*> =0: Successful exit 249*> <0: If INFO = -i, the i-th argument had an illegal value. 250*> >0: The matrix pairs (A, D) and (B, E) have common or very 251*> close eigenvalues. 252*> \endverbatim 253* 254* Authors: 255* ======== 256* 257*> \author Univ. of Tennessee 258*> \author Univ. of California Berkeley 259*> \author Univ. of Colorado Denver 260*> \author NAG Ltd. 261* 262*> \ingroup realSYauxiliary 263* 264*> \par Contributors: 265* ================== 266*> 267*> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 268*> Umea University, S-901 87 Umea, Sweden. 269* 270* ===================================================================== 271 SUBROUTINE STGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 272 $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, 273 $ IWORK, PQ, INFO ) 274* 275* -- LAPACK auxiliary routine -- 276* -- LAPACK is a software package provided by Univ. of Tennessee, -- 277* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 278* 279* .. Scalar Arguments .. 280 CHARACTER TRANS 281 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N, 282 $ PQ 283 REAL RDSCAL, RDSUM, SCALE 284* .. 285* .. Array Arguments .. 286 INTEGER IWORK( * ) 287 REAL A( LDA, * ), B( LDB, * ), C( LDC, * ), 288 $ D( LDD, * ), E( LDE, * ), F( LDF, * ) 289* .. 290* 291* ===================================================================== 292* Replaced various illegal calls to SCOPY by calls to SLASET. 293* Sven Hammarling, 27/5/02. 294* 295* .. Parameters .. 296 INTEGER LDZ 297 PARAMETER ( LDZ = 8 ) 298 REAL ZERO, ONE 299 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 300* .. 301* .. Local Scalars .. 302 LOGICAL NOTRAN 303 INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1, 304 $ K, MB, NB, P, Q, ZDIM 305 REAL ALPHA, SCALOC 306* .. 307* .. Local Arrays .. 308 INTEGER IPIV( LDZ ), JPIV( LDZ ) 309 REAL RHS( LDZ ), Z( LDZ, LDZ ) 310* .. 311* .. External Functions .. 312 LOGICAL LSAME 313 EXTERNAL LSAME 314* .. 315* .. External Subroutines .. 316 EXTERNAL SAXPY, SCOPY, SGEMM, SGEMV, SGER, SGESC2, 317 $ SGETC2, SSCAL, SLASET, SLATDF, XERBLA 318* .. 319* .. Intrinsic Functions .. 320 INTRINSIC MAX 321* .. 322* .. Executable Statements .. 323* 324* Decode and test input parameters 325* 326 INFO = 0 327 IERR = 0 328 NOTRAN = LSAME( TRANS, 'N' ) 329 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 330 INFO = -1 331 ELSE IF( NOTRAN ) THEN 332 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN 333 INFO = -2 334 END IF 335 END IF 336 IF( INFO.EQ.0 ) THEN 337 IF( M.LE.0 ) THEN 338 INFO = -3 339 ELSE IF( N.LE.0 ) THEN 340 INFO = -4 341 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 342 INFO = -6 343 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 344 INFO = -8 345 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 346 INFO = -10 347 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 348 INFO = -12 349 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 350 INFO = -14 351 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 352 INFO = -16 353 END IF 354 END IF 355 IF( INFO.NE.0 ) THEN 356 CALL XERBLA( 'STGSY2', -INFO ) 357 RETURN 358 END IF 359* 360* Determine block structure of A 361* 362 PQ = 0 363 P = 0 364 I = 1 365 10 CONTINUE 366 IF( I.GT.M ) 367 $ GO TO 20 368 P = P + 1 369 IWORK( P ) = I 370 IF( I.EQ.M ) 371 $ GO TO 20 372 IF( A( I+1, I ).NE.ZERO ) THEN 373 I = I + 2 374 ELSE 375 I = I + 1 376 END IF 377 GO TO 10 378 20 CONTINUE 379 IWORK( P+1 ) = M + 1 380* 381* Determine block structure of B 382* 383 Q = P + 1 384 J = 1 385 30 CONTINUE 386 IF( J.GT.N ) 387 $ GO TO 40 388 Q = Q + 1 389 IWORK( Q ) = J 390 IF( J.EQ.N ) 391 $ GO TO 40 392 IF( B( J+1, J ).NE.ZERO ) THEN 393 J = J + 2 394 ELSE 395 J = J + 1 396 END IF 397 GO TO 30 398 40 CONTINUE 399 IWORK( Q+1 ) = N + 1 400 PQ = P*( Q-P-1 ) 401* 402 IF( NOTRAN ) THEN 403* 404* Solve (I, J) - subsystem 405* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) 406* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) 407* for I = P, P - 1, ..., 1; J = 1, 2, ..., Q 408* 409 SCALE = ONE 410 SCALOC = ONE 411 DO 120 J = P + 2, Q 412 JS = IWORK( J ) 413 JSP1 = JS + 1 414 JE = IWORK( J+1 ) - 1 415 NB = JE - JS + 1 416 DO 110 I = P, 1, -1 417* 418 IS = IWORK( I ) 419 ISP1 = IS + 1 420 IE = IWORK( I+1 ) - 1 421 MB = IE - IS + 1 422 ZDIM = MB*NB*2 423* 424 IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN 425* 426* Build a 2-by-2 system Z * x = RHS 427* 428 Z( 1, 1 ) = A( IS, IS ) 429 Z( 2, 1 ) = D( IS, IS ) 430 Z( 1, 2 ) = -B( JS, JS ) 431 Z( 2, 2 ) = -E( JS, JS ) 432* 433* Set up right hand side(s) 434* 435 RHS( 1 ) = C( IS, JS ) 436 RHS( 2 ) = F( IS, JS ) 437* 438* Solve Z * x = RHS 439* 440 CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 441 IF( IERR.GT.0 ) 442 $ INFO = IERR 443* 444 IF( IJOB.EQ.0 ) THEN 445 CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 446 $ SCALOC ) 447 IF( SCALOC.NE.ONE ) THEN 448 DO 50 K = 1, N 449 CALL SSCAL( M, SCALOC, C( 1, K ), 1 ) 450 CALL SSCAL( M, SCALOC, F( 1, K ), 1 ) 451 50 CONTINUE 452 SCALE = SCALE*SCALOC 453 END IF 454 ELSE 455 CALL SLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 456 $ RDSCAL, IPIV, JPIV ) 457 END IF 458* 459* Unpack solution vector(s) 460* 461 C( IS, JS ) = RHS( 1 ) 462 F( IS, JS ) = RHS( 2 ) 463* 464* Substitute R(I, J) and L(I, J) into remaining 465* equation. 466* 467 IF( I.GT.1 ) THEN 468 ALPHA = -RHS( 1 ) 469 CALL SAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ), 470 $ 1 ) 471 CALL SAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ), 472 $ 1 ) 473 END IF 474 IF( J.LT.Q ) THEN 475 CALL SAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB, 476 $ C( IS, JE+1 ), LDC ) 477 CALL SAXPY( N-JE, RHS( 2 ), E( JS, JE+1 ), LDE, 478 $ F( IS, JE+1 ), LDF ) 479 END IF 480* 481 ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN 482* 483* Build a 4-by-4 system Z * x = RHS 484* 485 Z( 1, 1 ) = A( IS, IS ) 486 Z( 2, 1 ) = ZERO 487 Z( 3, 1 ) = D( IS, IS ) 488 Z( 4, 1 ) = ZERO 489* 490 Z( 1, 2 ) = ZERO 491 Z( 2, 2 ) = A( IS, IS ) 492 Z( 3, 2 ) = ZERO 493 Z( 4, 2 ) = D( IS, IS ) 494* 495 Z( 1, 3 ) = -B( JS, JS ) 496 Z( 2, 3 ) = -B( JS, JSP1 ) 497 Z( 3, 3 ) = -E( JS, JS ) 498 Z( 4, 3 ) = -E( JS, JSP1 ) 499* 500 Z( 1, 4 ) = -B( JSP1, JS ) 501 Z( 2, 4 ) = -B( JSP1, JSP1 ) 502 Z( 3, 4 ) = ZERO 503 Z( 4, 4 ) = -E( JSP1, JSP1 ) 504* 505* Set up right hand side(s) 506* 507 RHS( 1 ) = C( IS, JS ) 508 RHS( 2 ) = C( IS, JSP1 ) 509 RHS( 3 ) = F( IS, JS ) 510 RHS( 4 ) = F( IS, JSP1 ) 511* 512* Solve Z * x = RHS 513* 514 CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 515 IF( IERR.GT.0 ) 516 $ INFO = IERR 517* 518 IF( IJOB.EQ.0 ) THEN 519 CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 520 $ SCALOC ) 521 IF( SCALOC.NE.ONE ) THEN 522 DO 60 K = 1, N 523 CALL SSCAL( M, SCALOC, C( 1, K ), 1 ) 524 CALL SSCAL( M, SCALOC, F( 1, K ), 1 ) 525 60 CONTINUE 526 SCALE = SCALE*SCALOC 527 END IF 528 ELSE 529 CALL SLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 530 $ RDSCAL, IPIV, JPIV ) 531 END IF 532* 533* Unpack solution vector(s) 534* 535 C( IS, JS ) = RHS( 1 ) 536 C( IS, JSP1 ) = RHS( 2 ) 537 F( IS, JS ) = RHS( 3 ) 538 F( IS, JSP1 ) = RHS( 4 ) 539* 540* Substitute R(I, J) and L(I, J) into remaining 541* equation. 542* 543 IF( I.GT.1 ) THEN 544 CALL SGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ), 545 $ 1, C( 1, JS ), LDC ) 546 CALL SGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ), 547 $ 1, F( 1, JS ), LDF ) 548 END IF 549 IF( J.LT.Q ) THEN 550 CALL SAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB, 551 $ C( IS, JE+1 ), LDC ) 552 CALL SAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE, 553 $ F( IS, JE+1 ), LDF ) 554 CALL SAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB, 555 $ C( IS, JE+1 ), LDC ) 556 CALL SAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE, 557 $ F( IS, JE+1 ), LDF ) 558 END IF 559* 560 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN 561* 562* Build a 4-by-4 system Z * x = RHS 563* 564 Z( 1, 1 ) = A( IS, IS ) 565 Z( 2, 1 ) = A( ISP1, IS ) 566 Z( 3, 1 ) = D( IS, IS ) 567 Z( 4, 1 ) = ZERO 568* 569 Z( 1, 2 ) = A( IS, ISP1 ) 570 Z( 2, 2 ) = A( ISP1, ISP1 ) 571 Z( 3, 2 ) = D( IS, ISP1 ) 572 Z( 4, 2 ) = D( ISP1, ISP1 ) 573* 574 Z( 1, 3 ) = -B( JS, JS ) 575 Z( 2, 3 ) = ZERO 576 Z( 3, 3 ) = -E( JS, JS ) 577 Z( 4, 3 ) = ZERO 578* 579 Z( 1, 4 ) = ZERO 580 Z( 2, 4 ) = -B( JS, JS ) 581 Z( 3, 4 ) = ZERO 582 Z( 4, 4 ) = -E( JS, JS ) 583* 584* Set up right hand side(s) 585* 586 RHS( 1 ) = C( IS, JS ) 587 RHS( 2 ) = C( ISP1, JS ) 588 RHS( 3 ) = F( IS, JS ) 589 RHS( 4 ) = F( ISP1, JS ) 590* 591* Solve Z * x = RHS 592* 593 CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 594 IF( IERR.GT.0 ) 595 $ INFO = IERR 596 IF( IJOB.EQ.0 ) THEN 597 CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 598 $ SCALOC ) 599 IF( SCALOC.NE.ONE ) THEN 600 DO 70 K = 1, N 601 CALL SSCAL( M, SCALOC, C( 1, K ), 1 ) 602 CALL SSCAL( M, SCALOC, F( 1, K ), 1 ) 603 70 CONTINUE 604 SCALE = SCALE*SCALOC 605 END IF 606 ELSE 607 CALL SLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 608 $ RDSCAL, IPIV, JPIV ) 609 END IF 610* 611* Unpack solution vector(s) 612* 613 C( IS, JS ) = RHS( 1 ) 614 C( ISP1, JS ) = RHS( 2 ) 615 F( IS, JS ) = RHS( 3 ) 616 F( ISP1, JS ) = RHS( 4 ) 617* 618* Substitute R(I, J) and L(I, J) into remaining 619* equation. 620* 621 IF( I.GT.1 ) THEN 622 CALL SGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA, 623 $ RHS( 1 ), 1, ONE, C( 1, JS ), 1 ) 624 CALL SGEMV( 'N', IS-1, MB, -ONE, D( 1, IS ), LDD, 625 $ RHS( 1 ), 1, ONE, F( 1, JS ), 1 ) 626 END IF 627 IF( J.LT.Q ) THEN 628 CALL SGER( MB, N-JE, ONE, RHS( 3 ), 1, 629 $ B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC ) 630 CALL SGER( MB, N-JE, ONE, RHS( 3 ), 1, 631 $ E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF ) 632 END IF 633* 634 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN 635* 636* Build an 8-by-8 system Z * x = RHS 637* 638 CALL SLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ ) 639* 640 Z( 1, 1 ) = A( IS, IS ) 641 Z( 2, 1 ) = A( ISP1, IS ) 642 Z( 5, 1 ) = D( IS, IS ) 643* 644 Z( 1, 2 ) = A( IS, ISP1 ) 645 Z( 2, 2 ) = A( ISP1, ISP1 ) 646 Z( 5, 2 ) = D( IS, ISP1 ) 647 Z( 6, 2 ) = D( ISP1, ISP1 ) 648* 649 Z( 3, 3 ) = A( IS, IS ) 650 Z( 4, 3 ) = A( ISP1, IS ) 651 Z( 7, 3 ) = D( IS, IS ) 652* 653 Z( 3, 4 ) = A( IS, ISP1 ) 654 Z( 4, 4 ) = A( ISP1, ISP1 ) 655 Z( 7, 4 ) = D( IS, ISP1 ) 656 Z( 8, 4 ) = D( ISP1, ISP1 ) 657* 658 Z( 1, 5 ) = -B( JS, JS ) 659 Z( 3, 5 ) = -B( JS, JSP1 ) 660 Z( 5, 5 ) = -E( JS, JS ) 661 Z( 7, 5 ) = -E( JS, JSP1 ) 662* 663 Z( 2, 6 ) = -B( JS, JS ) 664 Z( 4, 6 ) = -B( JS, JSP1 ) 665 Z( 6, 6 ) = -E( JS, JS ) 666 Z( 8, 6 ) = -E( JS, JSP1 ) 667* 668 Z( 1, 7 ) = -B( JSP1, JS ) 669 Z( 3, 7 ) = -B( JSP1, JSP1 ) 670 Z( 7, 7 ) = -E( JSP1, JSP1 ) 671* 672 Z( 2, 8 ) = -B( JSP1, JS ) 673 Z( 4, 8 ) = -B( JSP1, JSP1 ) 674 Z( 8, 8 ) = -E( JSP1, JSP1 ) 675* 676* Set up right hand side(s) 677* 678 K = 1 679 II = MB*NB + 1 680 DO 80 JJ = 0, NB - 1 681 CALL SCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 ) 682 CALL SCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 ) 683 K = K + MB 684 II = II + MB 685 80 CONTINUE 686* 687* Solve Z * x = RHS 688* 689 CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 690 IF( IERR.GT.0 ) 691 $ INFO = IERR 692 IF( IJOB.EQ.0 ) THEN 693 CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 694 $ SCALOC ) 695 IF( SCALOC.NE.ONE ) THEN 696 DO 90 K = 1, N 697 CALL SSCAL( M, SCALOC, C( 1, K ), 1 ) 698 CALL SSCAL( M, SCALOC, F( 1, K ), 1 ) 699 90 CONTINUE 700 SCALE = SCALE*SCALOC 701 END IF 702 ELSE 703 CALL SLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 704 $ RDSCAL, IPIV, JPIV ) 705 END IF 706* 707* Unpack solution vector(s) 708* 709 K = 1 710 II = MB*NB + 1 711 DO 100 JJ = 0, NB - 1 712 CALL SCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 ) 713 CALL SCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 ) 714 K = K + MB 715 II = II + MB 716 100 CONTINUE 717* 718* Substitute R(I, J) and L(I, J) into remaining 719* equation. 720* 721 IF( I.GT.1 ) THEN 722 CALL SGEMM( 'N', 'N', IS-1, NB, MB, -ONE, 723 $ A( 1, IS ), LDA, RHS( 1 ), MB, ONE, 724 $ C( 1, JS ), LDC ) 725 CALL SGEMM( 'N', 'N', IS-1, NB, MB, -ONE, 726 $ D( 1, IS ), LDD, RHS( 1 ), MB, ONE, 727 $ F( 1, JS ), LDF ) 728 END IF 729 IF( J.LT.Q ) THEN 730 K = MB*NB + 1 731 CALL SGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ), 732 $ MB, B( JS, JE+1 ), LDB, ONE, 733 $ C( IS, JE+1 ), LDC ) 734 CALL SGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ), 735 $ MB, E( JS, JE+1 ), LDE, ONE, 736 $ F( IS, JE+1 ), LDF ) 737 END IF 738* 739 END IF 740* 741 110 CONTINUE 742 120 CONTINUE 743 ELSE 744* 745* Solve (I, J) - subsystem 746* A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J) 747* R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 748* for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1 749* 750 SCALE = ONE 751 SCALOC = ONE 752 DO 200 I = 1, P 753* 754 IS = IWORK( I ) 755 ISP1 = IS + 1 756 IE = IWORK( I+1 ) - 1 757 MB = IE - IS + 1 758 DO 190 J = Q, P + 2, -1 759* 760 JS = IWORK( J ) 761 JSP1 = JS + 1 762 JE = IWORK( J+1 ) - 1 763 NB = JE - JS + 1 764 ZDIM = MB*NB*2 765 IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN 766* 767* Build a 2-by-2 system Z**T * x = RHS 768* 769 Z( 1, 1 ) = A( IS, IS ) 770 Z( 2, 1 ) = -B( JS, JS ) 771 Z( 1, 2 ) = D( IS, IS ) 772 Z( 2, 2 ) = -E( JS, JS ) 773* 774* Set up right hand side(s) 775* 776 RHS( 1 ) = C( IS, JS ) 777 RHS( 2 ) = F( IS, JS ) 778* 779* Solve Z**T * x = RHS 780* 781 CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 782 IF( IERR.GT.0 ) 783 $ INFO = IERR 784* 785 CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 786 IF( SCALOC.NE.ONE ) THEN 787 DO 130 K = 1, N 788 CALL SSCAL( M, SCALOC, C( 1, K ), 1 ) 789 CALL SSCAL( M, SCALOC, F( 1, K ), 1 ) 790 130 CONTINUE 791 SCALE = SCALE*SCALOC 792 END IF 793* 794* Unpack solution vector(s) 795* 796 C( IS, JS ) = RHS( 1 ) 797 F( IS, JS ) = RHS( 2 ) 798* 799* Substitute R(I, J) and L(I, J) into remaining 800* equation. 801* 802 IF( J.GT.P+2 ) THEN 803 ALPHA = RHS( 1 ) 804 CALL SAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ), 805 $ LDF ) 806 ALPHA = RHS( 2 ) 807 CALL SAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ), 808 $ LDF ) 809 END IF 810 IF( I.LT.P ) THEN 811 ALPHA = -RHS( 1 ) 812 CALL SAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA, 813 $ C( IE+1, JS ), 1 ) 814 ALPHA = -RHS( 2 ) 815 CALL SAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD, 816 $ C( IE+1, JS ), 1 ) 817 END IF 818* 819 ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN 820* 821* Build a 4-by-4 system Z**T * x = RHS 822* 823 Z( 1, 1 ) = A( IS, IS ) 824 Z( 2, 1 ) = ZERO 825 Z( 3, 1 ) = -B( JS, JS ) 826 Z( 4, 1 ) = -B( JSP1, JS ) 827* 828 Z( 1, 2 ) = ZERO 829 Z( 2, 2 ) = A( IS, IS ) 830 Z( 3, 2 ) = -B( JS, JSP1 ) 831 Z( 4, 2 ) = -B( JSP1, JSP1 ) 832* 833 Z( 1, 3 ) = D( IS, IS ) 834 Z( 2, 3 ) = ZERO 835 Z( 3, 3 ) = -E( JS, JS ) 836 Z( 4, 3 ) = ZERO 837* 838 Z( 1, 4 ) = ZERO 839 Z( 2, 4 ) = D( IS, IS ) 840 Z( 3, 4 ) = -E( JS, JSP1 ) 841 Z( 4, 4 ) = -E( JSP1, JSP1 ) 842* 843* Set up right hand side(s) 844* 845 RHS( 1 ) = C( IS, JS ) 846 RHS( 2 ) = C( IS, JSP1 ) 847 RHS( 3 ) = F( IS, JS ) 848 RHS( 4 ) = F( IS, JSP1 ) 849* 850* Solve Z**T * x = RHS 851* 852 CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 853 IF( IERR.GT.0 ) 854 $ INFO = IERR 855 CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 856 IF( SCALOC.NE.ONE ) THEN 857 DO 140 K = 1, N 858 CALL SSCAL( M, SCALOC, C( 1, K ), 1 ) 859 CALL SSCAL( M, SCALOC, F( 1, K ), 1 ) 860 140 CONTINUE 861 SCALE = SCALE*SCALOC 862 END IF 863* 864* Unpack solution vector(s) 865* 866 C( IS, JS ) = RHS( 1 ) 867 C( IS, JSP1 ) = RHS( 2 ) 868 F( IS, JS ) = RHS( 3 ) 869 F( IS, JSP1 ) = RHS( 4 ) 870* 871* Substitute R(I, J) and L(I, J) into remaining 872* equation. 873* 874 IF( J.GT.P+2 ) THEN 875 CALL SAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1, 876 $ F( IS, 1 ), LDF ) 877 CALL SAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1, 878 $ F( IS, 1 ), LDF ) 879 CALL SAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1, 880 $ F( IS, 1 ), LDF ) 881 CALL SAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1, 882 $ F( IS, 1 ), LDF ) 883 END IF 884 IF( I.LT.P ) THEN 885 CALL SGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA, 886 $ RHS( 1 ), 1, C( IE+1, JS ), LDC ) 887 CALL SGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD, 888 $ RHS( 3 ), 1, C( IE+1, JS ), LDC ) 889 END IF 890* 891 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN 892* 893* Build a 4-by-4 system Z**T * x = RHS 894* 895 Z( 1, 1 ) = A( IS, IS ) 896 Z( 2, 1 ) = A( IS, ISP1 ) 897 Z( 3, 1 ) = -B( JS, JS ) 898 Z( 4, 1 ) = ZERO 899* 900 Z( 1, 2 ) = A( ISP1, IS ) 901 Z( 2, 2 ) = A( ISP1, ISP1 ) 902 Z( 3, 2 ) = ZERO 903 Z( 4, 2 ) = -B( JS, JS ) 904* 905 Z( 1, 3 ) = D( IS, IS ) 906 Z( 2, 3 ) = D( IS, ISP1 ) 907 Z( 3, 3 ) = -E( JS, JS ) 908 Z( 4, 3 ) = ZERO 909* 910 Z( 1, 4 ) = ZERO 911 Z( 2, 4 ) = D( ISP1, ISP1 ) 912 Z( 3, 4 ) = ZERO 913 Z( 4, 4 ) = -E( JS, JS ) 914* 915* Set up right hand side(s) 916* 917 RHS( 1 ) = C( IS, JS ) 918 RHS( 2 ) = C( ISP1, JS ) 919 RHS( 3 ) = F( IS, JS ) 920 RHS( 4 ) = F( ISP1, JS ) 921* 922* Solve Z**T * x = RHS 923* 924 CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 925 IF( IERR.GT.0 ) 926 $ INFO = IERR 927* 928 CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 929 IF( SCALOC.NE.ONE ) THEN 930 DO 150 K = 1, N 931 CALL SSCAL( M, SCALOC, C( 1, K ), 1 ) 932 CALL SSCAL( M, SCALOC, F( 1, K ), 1 ) 933 150 CONTINUE 934 SCALE = SCALE*SCALOC 935 END IF 936* 937* Unpack solution vector(s) 938* 939 C( IS, JS ) = RHS( 1 ) 940 C( ISP1, JS ) = RHS( 2 ) 941 F( IS, JS ) = RHS( 3 ) 942 F( ISP1, JS ) = RHS( 4 ) 943* 944* Substitute R(I, J) and L(I, J) into remaining 945* equation. 946* 947 IF( J.GT.P+2 ) THEN 948 CALL SGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ), 949 $ 1, F( IS, 1 ), LDF ) 950 CALL SGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ), 951 $ 1, F( IS, 1 ), LDF ) 952 END IF 953 IF( I.LT.P ) THEN 954 CALL SGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ), 955 $ LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ), 956 $ 1 ) 957 CALL SGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ), 958 $ LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ), 959 $ 1 ) 960 END IF 961* 962 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN 963* 964* Build an 8-by-8 system Z**T * x = RHS 965* 966 CALL SLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ ) 967* 968 Z( 1, 1 ) = A( IS, IS ) 969 Z( 2, 1 ) = A( IS, ISP1 ) 970 Z( 5, 1 ) = -B( JS, JS ) 971 Z( 7, 1 ) = -B( JSP1, JS ) 972* 973 Z( 1, 2 ) = A( ISP1, IS ) 974 Z( 2, 2 ) = A( ISP1, ISP1 ) 975 Z( 6, 2 ) = -B( JS, JS ) 976 Z( 8, 2 ) = -B( JSP1, JS ) 977* 978 Z( 3, 3 ) = A( IS, IS ) 979 Z( 4, 3 ) = A( IS, ISP1 ) 980 Z( 5, 3 ) = -B( JS, JSP1 ) 981 Z( 7, 3 ) = -B( JSP1, JSP1 ) 982* 983 Z( 3, 4 ) = A( ISP1, IS ) 984 Z( 4, 4 ) = A( ISP1, ISP1 ) 985 Z( 6, 4 ) = -B( JS, JSP1 ) 986 Z( 8, 4 ) = -B( JSP1, JSP1 ) 987* 988 Z( 1, 5 ) = D( IS, IS ) 989 Z( 2, 5 ) = D( IS, ISP1 ) 990 Z( 5, 5 ) = -E( JS, JS ) 991* 992 Z( 2, 6 ) = D( ISP1, ISP1 ) 993 Z( 6, 6 ) = -E( JS, JS ) 994* 995 Z( 3, 7 ) = D( IS, IS ) 996 Z( 4, 7 ) = D( IS, ISP1 ) 997 Z( 5, 7 ) = -E( JS, JSP1 ) 998 Z( 7, 7 ) = -E( JSP1, JSP1 ) 999* 1000 Z( 4, 8 ) = D( ISP1, ISP1 ) 1001 Z( 6, 8 ) = -E( JS, JSP1 ) 1002 Z( 8, 8 ) = -E( JSP1, JSP1 ) 1003* 1004* Set up right hand side(s) 1005* 1006 K = 1 1007 II = MB*NB + 1 1008 DO 160 JJ = 0, NB - 1 1009 CALL SCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 ) 1010 CALL SCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 ) 1011 K = K + MB 1012 II = II + MB 1013 160 CONTINUE 1014* 1015* 1016* Solve Z**T * x = RHS 1017* 1018 CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 1019 IF( IERR.GT.0 ) 1020 $ INFO = IERR 1021* 1022 CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 1023 IF( SCALOC.NE.ONE ) THEN 1024 DO 170 K = 1, N 1025 CALL SSCAL( M, SCALOC, C( 1, K ), 1 ) 1026 CALL SSCAL( M, SCALOC, F( 1, K ), 1 ) 1027 170 CONTINUE 1028 SCALE = SCALE*SCALOC 1029 END IF 1030* 1031* Unpack solution vector(s) 1032* 1033 K = 1 1034 II = MB*NB + 1 1035 DO 180 JJ = 0, NB - 1 1036 CALL SCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 ) 1037 CALL SCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 ) 1038 K = K + MB 1039 II = II + MB 1040 180 CONTINUE 1041* 1042* Substitute R(I, J) and L(I, J) into remaining 1043* equation. 1044* 1045 IF( J.GT.P+2 ) THEN 1046 CALL SGEMM( 'N', 'T', MB, JS-1, NB, ONE, 1047 $ C( IS, JS ), LDC, B( 1, JS ), LDB, ONE, 1048 $ F( IS, 1 ), LDF ) 1049 CALL SGEMM( 'N', 'T', MB, JS-1, NB, ONE, 1050 $ F( IS, JS ), LDF, E( 1, JS ), LDE, ONE, 1051 $ F( IS, 1 ), LDF ) 1052 END IF 1053 IF( I.LT.P ) THEN 1054 CALL SGEMM( 'T', 'N', M-IE, NB, MB, -ONE, 1055 $ A( IS, IE+1 ), LDA, C( IS, JS ), LDC, 1056 $ ONE, C( IE+1, JS ), LDC ) 1057 CALL SGEMM( 'T', 'N', M-IE, NB, MB, -ONE, 1058 $ D( IS, IE+1 ), LDD, F( IS, JS ), LDF, 1059 $ ONE, C( IE+1, JS ), LDC ) 1060 END IF 1061* 1062 END IF 1063* 1064 190 CONTINUE 1065 200 CONTINUE 1066* 1067 END IF 1068 RETURN 1069* 1070* End of STGSY2 1071* 1072 END 1073