1*> \brief \b DTRSYL 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DTRSYL + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsyl.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsyl.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsyl.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 22* LDC, SCALE, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER TRANA, TRANB 26* INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 27* DOUBLE PRECISION SCALE 28* .. 29* .. Array Arguments .. 30* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DTRSYL solves the real Sylvester matrix equation: 40*> 41*> op(A)*X + X*op(B) = scale*C or 42*> op(A)*X - X*op(B) = scale*C, 43*> 44*> where op(A) = A or A**T, and A and B are both upper quasi- 45*> triangular. A is M-by-M and B is N-by-N; the right hand side C and 46*> the solution X are M-by-N; and scale is an output scale factor, set 47*> <= 1 to avoid overflow in X. 48*> 49*> A and B must be in Schur canonical form (as returned by DHSEQR), that 50*> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; 51*> each 2-by-2 diagonal block has its diagonal elements equal and its 52*> off-diagonal elements of opposite sign. 53*> \endverbatim 54* 55* Arguments: 56* ========== 57* 58*> \param[in] TRANA 59*> \verbatim 60*> TRANA is CHARACTER*1 61*> Specifies the option op(A): 62*> = 'N': op(A) = A (No transpose) 63*> = 'T': op(A) = A**T (Transpose) 64*> = 'C': op(A) = A**H (Conjugate transpose = Transpose) 65*> \endverbatim 66*> 67*> \param[in] TRANB 68*> \verbatim 69*> TRANB is CHARACTER*1 70*> Specifies the option op(B): 71*> = 'N': op(B) = B (No transpose) 72*> = 'T': op(B) = B**T (Transpose) 73*> = 'C': op(B) = B**H (Conjugate transpose = Transpose) 74*> \endverbatim 75*> 76*> \param[in] ISGN 77*> \verbatim 78*> ISGN is INTEGER 79*> Specifies the sign in the equation: 80*> = +1: solve op(A)*X + X*op(B) = scale*C 81*> = -1: solve op(A)*X - X*op(B) = scale*C 82*> \endverbatim 83*> 84*> \param[in] M 85*> \verbatim 86*> M is INTEGER 87*> The order of the matrix A, and the number of rows in the 88*> matrices X and C. M >= 0. 89*> \endverbatim 90*> 91*> \param[in] N 92*> \verbatim 93*> N is INTEGER 94*> The order of the matrix B, and the number of columns in the 95*> matrices X and C. N >= 0. 96*> \endverbatim 97*> 98*> \param[in] A 99*> \verbatim 100*> A is DOUBLE PRECISION array, dimension (LDA,M) 101*> The upper quasi-triangular matrix A, in Schur canonical form. 102*> \endverbatim 103*> 104*> \param[in] LDA 105*> \verbatim 106*> LDA is INTEGER 107*> The leading dimension of the array A. LDA >= max(1,M). 108*> \endverbatim 109*> 110*> \param[in] B 111*> \verbatim 112*> B is DOUBLE PRECISION array, dimension (LDB,N) 113*> The upper quasi-triangular matrix B, in Schur canonical form. 114*> \endverbatim 115*> 116*> \param[in] LDB 117*> \verbatim 118*> LDB is INTEGER 119*> The leading dimension of the array B. LDB >= max(1,N). 120*> \endverbatim 121*> 122*> \param[in,out] C 123*> \verbatim 124*> C is DOUBLE PRECISION array, dimension (LDC,N) 125*> On entry, the M-by-N right hand side matrix C. 126*> On exit, C is overwritten by the solution matrix X. 127*> \endverbatim 128*> 129*> \param[in] LDC 130*> \verbatim 131*> LDC is INTEGER 132*> The leading dimension of the array C. LDC >= max(1,M) 133*> \endverbatim 134*> 135*> \param[out] SCALE 136*> \verbatim 137*> SCALE is DOUBLE PRECISION 138*> The scale factor, scale, set <= 1 to avoid overflow in X. 139*> \endverbatim 140*> 141*> \param[out] INFO 142*> \verbatim 143*> INFO is INTEGER 144*> = 0: successful exit 145*> < 0: if INFO = -i, the i-th argument had an illegal value 146*> = 1: A and B have common or very close eigenvalues; perturbed 147*> values were used to solve the equation (but the matrices 148*> A and B are unchanged). 149*> \endverbatim 150* 151* Authors: 152* ======== 153* 154*> \author Univ. of Tennessee 155*> \author Univ. of California Berkeley 156*> \author Univ. of Colorado Denver 157*> \author NAG Ltd. 158* 159*> \ingroup doubleSYcomputational 160* 161* ===================================================================== 162 SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 163 $ LDC, SCALE, INFO ) 164* 165* -- LAPACK computational routine -- 166* -- LAPACK is a software package provided by Univ. of Tennessee, -- 167* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 168* 169* .. Scalar Arguments .. 170 CHARACTER TRANA, TRANB 171 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 172 DOUBLE PRECISION SCALE 173* .. 174* .. Array Arguments .. 175 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) 176* .. 177* 178* ===================================================================== 179* 180* .. Parameters .. 181 DOUBLE PRECISION ZERO, ONE 182 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 183* .. 184* .. Local Scalars .. 185 LOGICAL NOTRNA, NOTRNB 186 INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT 187 DOUBLE PRECISION A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, 188 $ SMLNUM, SUML, SUMR, XNORM 189* .. 190* .. Local Arrays .. 191 DOUBLE PRECISION DUM( 1 ), VEC( 2, 2 ), X( 2, 2 ) 192* .. 193* .. External Functions .. 194 LOGICAL LSAME 195 DOUBLE PRECISION DDOT, DLAMCH, DLANGE 196 EXTERNAL LSAME, DDOT, DLAMCH, DLANGE 197* .. 198* .. External Subroutines .. 199 EXTERNAL DLABAD, DLALN2, DLASY2, DSCAL, XERBLA 200* .. 201* .. Intrinsic Functions .. 202 INTRINSIC ABS, DBLE, MAX, MIN 203* .. 204* .. Executable Statements .. 205* 206* Decode and Test input parameters 207* 208 NOTRNA = LSAME( TRANA, 'N' ) 209 NOTRNB = LSAME( TRANB, 'N' ) 210* 211 INFO = 0 212 IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT. 213 $ LSAME( TRANA, 'C' ) ) THEN 214 INFO = -1 215 ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT. 216 $ LSAME( TRANB, 'C' ) ) THEN 217 INFO = -2 218 ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN 219 INFO = -3 220 ELSE IF( M.LT.0 ) THEN 221 INFO = -4 222 ELSE IF( N.LT.0 ) THEN 223 INFO = -5 224 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 225 INFO = -7 226 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 227 INFO = -9 228 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 229 INFO = -11 230 END IF 231 IF( INFO.NE.0 ) THEN 232 CALL XERBLA( 'DTRSYL', -INFO ) 233 RETURN 234 END IF 235* 236* Quick return if possible 237* 238 SCALE = ONE 239 IF( M.EQ.0 .OR. N.EQ.0 ) 240 $ RETURN 241* 242* Set constants to control overflow 243* 244 EPS = DLAMCH( 'P' ) 245 SMLNUM = DLAMCH( 'S' ) 246 BIGNUM = ONE / SMLNUM 247 CALL DLABAD( SMLNUM, BIGNUM ) 248 SMLNUM = SMLNUM*DBLE( M*N ) / EPS 249 BIGNUM = ONE / SMLNUM 250* 251 SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ), 252 $ EPS*DLANGE( 'M', N, N, B, LDB, DUM ) ) 253* 254 SGN = ISGN 255* 256 IF( NOTRNA .AND. NOTRNB ) THEN 257* 258* Solve A*X + ISGN*X*B = scale*C. 259* 260* The (K,L)th block of X is determined starting from 261* bottom-left corner column by column by 262* 263* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 264* 265* Where 266* M L-1 267* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]. 268* I=K+1 J=1 269* 270* Start column loop (index = L) 271* L1 (L2) : column index of the first (first) row of X(K,L). 272* 273 LNEXT = 1 274 DO 60 L = 1, N 275 IF( L.LT.LNEXT ) 276 $ GO TO 60 277 IF( L.EQ.N ) THEN 278 L1 = L 279 L2 = L 280 ELSE 281 IF( B( L+1, L ).NE.ZERO ) THEN 282 L1 = L 283 L2 = L + 1 284 LNEXT = L + 2 285 ELSE 286 L1 = L 287 L2 = L 288 LNEXT = L + 1 289 END IF 290 END IF 291* 292* Start row loop (index = K) 293* K1 (K2): row index of the first (last) row of X(K,L). 294* 295 KNEXT = M 296 DO 50 K = M, 1, -1 297 IF( K.GT.KNEXT ) 298 $ GO TO 50 299 IF( K.EQ.1 ) THEN 300 K1 = K 301 K2 = K 302 ELSE 303 IF( A( K, K-1 ).NE.ZERO ) THEN 304 K1 = K - 1 305 K2 = K 306 KNEXT = K - 2 307 ELSE 308 K1 = K 309 K2 = K 310 KNEXT = K - 1 311 END IF 312 END IF 313* 314 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 315 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 316 $ C( MIN( K1+1, M ), L1 ), 1 ) 317 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 318 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 319 SCALOC = ONE 320* 321 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 322 DA11 = ABS( A11 ) 323 IF( DA11.LE.SMIN ) THEN 324 A11 = SMIN 325 DA11 = SMIN 326 INFO = 1 327 END IF 328 DB = ABS( VEC( 1, 1 ) ) 329 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 330 IF( DB.GT.BIGNUM*DA11 ) 331 $ SCALOC = ONE / DB 332 END IF 333 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 334* 335 IF( SCALOC.NE.ONE ) THEN 336 DO 10 J = 1, N 337 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 338 10 CONTINUE 339 SCALE = SCALE*SCALOC 340 END IF 341 C( K1, L1 ) = X( 1, 1 ) 342* 343 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 344* 345 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 346 $ C( MIN( K2+1, M ), L1 ), 1 ) 347 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 348 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 349* 350 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 351 $ C( MIN( K2+1, M ), L1 ), 1 ) 352 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 353 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 354* 355 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ), 356 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 357 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 358 IF( IERR.NE.0 ) 359 $ INFO = 1 360* 361 IF( SCALOC.NE.ONE ) THEN 362 DO 20 J = 1, N 363 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 364 20 CONTINUE 365 SCALE = SCALE*SCALOC 366 END IF 367 C( K1, L1 ) = X( 1, 1 ) 368 C( K2, L1 ) = X( 2, 1 ) 369* 370 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 371* 372 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 373 $ C( MIN( K1+1, M ), L1 ), 1 ) 374 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 375 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 376* 377 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 378 $ C( MIN( K1+1, M ), L2 ), 1 ) 379 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 380 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 381* 382 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ), 383 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 384 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 385 IF( IERR.NE.0 ) 386 $ INFO = 1 387* 388 IF( SCALOC.NE.ONE ) THEN 389 DO 30 J = 1, N 390 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 391 30 CONTINUE 392 SCALE = SCALE*SCALOC 393 END IF 394 C( K1, L1 ) = X( 1, 1 ) 395 C( K1, L2 ) = X( 2, 1 ) 396* 397 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 398* 399 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 400 $ C( MIN( K2+1, M ), L1 ), 1 ) 401 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 402 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 403* 404 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 405 $ C( MIN( K2+1, M ), L2 ), 1 ) 406 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 407 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 408* 409 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 410 $ C( MIN( K2+1, M ), L1 ), 1 ) 411 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 412 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 413* 414 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 415 $ C( MIN( K2+1, M ), L2 ), 1 ) 416 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 417 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 418* 419 CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2, 420 $ A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC, 421 $ 2, SCALOC, X, 2, XNORM, IERR ) 422 IF( IERR.NE.0 ) 423 $ INFO = 1 424* 425 IF( SCALOC.NE.ONE ) THEN 426 DO 40 J = 1, N 427 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 428 40 CONTINUE 429 SCALE = SCALE*SCALOC 430 END IF 431 C( K1, L1 ) = X( 1, 1 ) 432 C( K1, L2 ) = X( 1, 2 ) 433 C( K2, L1 ) = X( 2, 1 ) 434 C( K2, L2 ) = X( 2, 2 ) 435 END IF 436* 437 50 CONTINUE 438* 439 60 CONTINUE 440* 441 ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN 442* 443* Solve A**T *X + ISGN*X*B = scale*C. 444* 445* The (K,L)th block of X is determined starting from 446* upper-left corner column by column by 447* 448* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 449* 450* Where 451* K-1 T L-1 452* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)] 453* I=1 J=1 454* 455* Start column loop (index = L) 456* L1 (L2): column index of the first (last) row of X(K,L) 457* 458 LNEXT = 1 459 DO 120 L = 1, N 460 IF( L.LT.LNEXT ) 461 $ GO TO 120 462 IF( L.EQ.N ) THEN 463 L1 = L 464 L2 = L 465 ELSE 466 IF( B( L+1, L ).NE.ZERO ) THEN 467 L1 = L 468 L2 = L + 1 469 LNEXT = L + 2 470 ELSE 471 L1 = L 472 L2 = L 473 LNEXT = L + 1 474 END IF 475 END IF 476* 477* Start row loop (index = K) 478* K1 (K2): row index of the first (last) row of X(K,L) 479* 480 KNEXT = 1 481 DO 110 K = 1, M 482 IF( K.LT.KNEXT ) 483 $ GO TO 110 484 IF( K.EQ.M ) THEN 485 K1 = K 486 K2 = K 487 ELSE 488 IF( A( K+1, K ).NE.ZERO ) THEN 489 K1 = K 490 K2 = K + 1 491 KNEXT = K + 2 492 ELSE 493 K1 = K 494 K2 = K 495 KNEXT = K + 1 496 END IF 497 END IF 498* 499 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 500 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 501 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 502 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 503 SCALOC = ONE 504* 505 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 506 DA11 = ABS( A11 ) 507 IF( DA11.LE.SMIN ) THEN 508 A11 = SMIN 509 DA11 = SMIN 510 INFO = 1 511 END IF 512 DB = ABS( VEC( 1, 1 ) ) 513 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 514 IF( DB.GT.BIGNUM*DA11 ) 515 $ SCALOC = ONE / DB 516 END IF 517 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 518* 519 IF( SCALOC.NE.ONE ) THEN 520 DO 70 J = 1, N 521 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 522 70 CONTINUE 523 SCALE = SCALE*SCALOC 524 END IF 525 C( K1, L1 ) = X( 1, 1 ) 526* 527 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 528* 529 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 530 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 531 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 532* 533 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 534 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 535 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 536* 537 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ), 538 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 539 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 540 IF( IERR.NE.0 ) 541 $ INFO = 1 542* 543 IF( SCALOC.NE.ONE ) THEN 544 DO 80 J = 1, N 545 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 546 80 CONTINUE 547 SCALE = SCALE*SCALOC 548 END IF 549 C( K1, L1 ) = X( 1, 1 ) 550 C( K2, L1 ) = X( 2, 1 ) 551* 552 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 553* 554 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 555 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 556 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 557* 558 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 559 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 560 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 561* 562 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ), 563 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 564 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 565 IF( IERR.NE.0 ) 566 $ INFO = 1 567* 568 IF( SCALOC.NE.ONE ) THEN 569 DO 90 J = 1, N 570 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 571 90 CONTINUE 572 SCALE = SCALE*SCALOC 573 END IF 574 C( K1, L1 ) = X( 1, 1 ) 575 C( K1, L2 ) = X( 2, 1 ) 576* 577 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 578* 579 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 580 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 581 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 582* 583 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 584 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 585 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 586* 587 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 588 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 589 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 590* 591 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 592 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 593 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 594* 595 CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ), 596 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 597 $ 2, XNORM, IERR ) 598 IF( IERR.NE.0 ) 599 $ INFO = 1 600* 601 IF( SCALOC.NE.ONE ) THEN 602 DO 100 J = 1, N 603 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 604 100 CONTINUE 605 SCALE = SCALE*SCALOC 606 END IF 607 C( K1, L1 ) = X( 1, 1 ) 608 C( K1, L2 ) = X( 1, 2 ) 609 C( K2, L1 ) = X( 2, 1 ) 610 C( K2, L2 ) = X( 2, 2 ) 611 END IF 612* 613 110 CONTINUE 614 120 CONTINUE 615* 616 ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN 617* 618* Solve A**T*X + ISGN*X*B**T = scale*C. 619* 620* The (K,L)th block of X is determined starting from 621* top-right corner column by column by 622* 623* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) 624* 625* Where 626* K-1 N 627* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. 628* I=1 J=L+1 629* 630* Start column loop (index = L) 631* L1 (L2): column index of the first (last) row of X(K,L) 632* 633 LNEXT = N 634 DO 180 L = N, 1, -1 635 IF( L.GT.LNEXT ) 636 $ GO TO 180 637 IF( L.EQ.1 ) THEN 638 L1 = L 639 L2 = L 640 ELSE 641 IF( B( L, L-1 ).NE.ZERO ) THEN 642 L1 = L - 1 643 L2 = L 644 LNEXT = L - 2 645 ELSE 646 L1 = L 647 L2 = L 648 LNEXT = L - 1 649 END IF 650 END IF 651* 652* Start row loop (index = K) 653* K1 (K2): row index of the first (last) row of X(K,L) 654* 655 KNEXT = 1 656 DO 170 K = 1, M 657 IF( K.LT.KNEXT ) 658 $ GO TO 170 659 IF( K.EQ.M ) THEN 660 K1 = K 661 K2 = K 662 ELSE 663 IF( A( K+1, K ).NE.ZERO ) THEN 664 K1 = K 665 K2 = K + 1 666 KNEXT = K + 2 667 ELSE 668 K1 = K 669 K2 = K 670 KNEXT = K + 1 671 END IF 672 END IF 673* 674 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 675 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 676 SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC, 677 $ B( L1, MIN( L1+1, N ) ), LDB ) 678 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 679 SCALOC = ONE 680* 681 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 682 DA11 = ABS( A11 ) 683 IF( DA11.LE.SMIN ) THEN 684 A11 = SMIN 685 DA11 = SMIN 686 INFO = 1 687 END IF 688 DB = ABS( VEC( 1, 1 ) ) 689 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 690 IF( DB.GT.BIGNUM*DA11 ) 691 $ SCALOC = ONE / DB 692 END IF 693 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 694* 695 IF( SCALOC.NE.ONE ) THEN 696 DO 130 J = 1, N 697 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 698 130 CONTINUE 699 SCALE = SCALE*SCALOC 700 END IF 701 C( K1, L1 ) = X( 1, 1 ) 702* 703 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 704* 705 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 706 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 707 $ B( L1, MIN( L2+1, N ) ), LDB ) 708 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 709* 710 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 711 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 712 $ B( L1, MIN( L2+1, N ) ), LDB ) 713 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 714* 715 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ), 716 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 717 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 718 IF( IERR.NE.0 ) 719 $ INFO = 1 720* 721 IF( SCALOC.NE.ONE ) THEN 722 DO 140 J = 1, N 723 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 724 140 CONTINUE 725 SCALE = SCALE*SCALOC 726 END IF 727 C( K1, L1 ) = X( 1, 1 ) 728 C( K2, L1 ) = X( 2, 1 ) 729* 730 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 731* 732 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 733 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 734 $ B( L1, MIN( L2+1, N ) ), LDB ) 735 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 736* 737 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 738 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 739 $ B( L2, MIN( L2+1, N ) ), LDB ) 740 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 741* 742 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ), 743 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 744 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 745 IF( IERR.NE.0 ) 746 $ INFO = 1 747* 748 IF( SCALOC.NE.ONE ) THEN 749 DO 150 J = 1, N 750 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 751 150 CONTINUE 752 SCALE = SCALE*SCALOC 753 END IF 754 C( K1, L1 ) = X( 1, 1 ) 755 C( K1, L2 ) = X( 2, 1 ) 756* 757 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 758* 759 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 760 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 761 $ B( L1, MIN( L2+1, N ) ), LDB ) 762 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 763* 764 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 765 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 766 $ B( L2, MIN( L2+1, N ) ), LDB ) 767 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 768* 769 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 770 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 771 $ B( L1, MIN( L2+1, N ) ), LDB ) 772 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 773* 774 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 775 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 776 $ B( L2, MIN( L2+1, N ) ), LDB ) 777 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 778* 779 CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ), 780 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 781 $ 2, XNORM, IERR ) 782 IF( IERR.NE.0 ) 783 $ INFO = 1 784* 785 IF( SCALOC.NE.ONE ) THEN 786 DO 160 J = 1, N 787 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 788 160 CONTINUE 789 SCALE = SCALE*SCALOC 790 END IF 791 C( K1, L1 ) = X( 1, 1 ) 792 C( K1, L2 ) = X( 1, 2 ) 793 C( K2, L1 ) = X( 2, 1 ) 794 C( K2, L2 ) = X( 2, 2 ) 795 END IF 796* 797 170 CONTINUE 798 180 CONTINUE 799* 800 ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN 801* 802* Solve A*X + ISGN*X*B**T = scale*C. 803* 804* The (K,L)th block of X is determined starting from 805* bottom-right corner column by column by 806* 807* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) 808* 809* Where 810* M N 811* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. 812* I=K+1 J=L+1 813* 814* Start column loop (index = L) 815* L1 (L2): column index of the first (last) row of X(K,L) 816* 817 LNEXT = N 818 DO 240 L = N, 1, -1 819 IF( L.GT.LNEXT ) 820 $ GO TO 240 821 IF( L.EQ.1 ) THEN 822 L1 = L 823 L2 = L 824 ELSE 825 IF( B( L, L-1 ).NE.ZERO ) THEN 826 L1 = L - 1 827 L2 = L 828 LNEXT = L - 2 829 ELSE 830 L1 = L 831 L2 = L 832 LNEXT = L - 1 833 END IF 834 END IF 835* 836* Start row loop (index = K) 837* K1 (K2): row index of the first (last) row of X(K,L) 838* 839 KNEXT = M 840 DO 230 K = M, 1, -1 841 IF( K.GT.KNEXT ) 842 $ GO TO 230 843 IF( K.EQ.1 ) THEN 844 K1 = K 845 K2 = K 846 ELSE 847 IF( A( K, K-1 ).NE.ZERO ) THEN 848 K1 = K - 1 849 K2 = K 850 KNEXT = K - 2 851 ELSE 852 K1 = K 853 K2 = K 854 KNEXT = K - 1 855 END IF 856 END IF 857* 858 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 859 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 860 $ C( MIN( K1+1, M ), L1 ), 1 ) 861 SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC, 862 $ B( L1, MIN( L1+1, N ) ), LDB ) 863 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 864 SCALOC = ONE 865* 866 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 867 DA11 = ABS( A11 ) 868 IF( DA11.LE.SMIN ) THEN 869 A11 = SMIN 870 DA11 = SMIN 871 INFO = 1 872 END IF 873 DB = ABS( VEC( 1, 1 ) ) 874 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 875 IF( DB.GT.BIGNUM*DA11 ) 876 $ SCALOC = ONE / DB 877 END IF 878 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 879* 880 IF( SCALOC.NE.ONE ) THEN 881 DO 190 J = 1, N 882 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 883 190 CONTINUE 884 SCALE = SCALE*SCALOC 885 END IF 886 C( K1, L1 ) = X( 1, 1 ) 887* 888 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 889* 890 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 891 $ C( MIN( K2+1, M ), L1 ), 1 ) 892 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 893 $ B( L1, MIN( L2+1, N ) ), LDB ) 894 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 895* 896 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 897 $ C( MIN( K2+1, M ), L1 ), 1 ) 898 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 899 $ B( L1, MIN( L2+1, N ) ), LDB ) 900 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 901* 902 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ), 903 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 904 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 905 IF( IERR.NE.0 ) 906 $ INFO = 1 907* 908 IF( SCALOC.NE.ONE ) THEN 909 DO 200 J = 1, N 910 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 911 200 CONTINUE 912 SCALE = SCALE*SCALOC 913 END IF 914 C( K1, L1 ) = X( 1, 1 ) 915 C( K2, L1 ) = X( 2, 1 ) 916* 917 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 918* 919 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 920 $ C( MIN( K1+1, M ), L1 ), 1 ) 921 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 922 $ B( L1, MIN( L2+1, N ) ), LDB ) 923 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 924* 925 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 926 $ C( MIN( K1+1, M ), L2 ), 1 ) 927 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 928 $ B( L2, MIN( L2+1, N ) ), LDB ) 929 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 930* 931 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ), 932 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 933 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 934 IF( IERR.NE.0 ) 935 $ INFO = 1 936* 937 IF( SCALOC.NE.ONE ) THEN 938 DO 210 J = 1, N 939 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 940 210 CONTINUE 941 SCALE = SCALE*SCALOC 942 END IF 943 C( K1, L1 ) = X( 1, 1 ) 944 C( K1, L2 ) = X( 2, 1 ) 945* 946 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 947* 948 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 949 $ C( MIN( K2+1, M ), L1 ), 1 ) 950 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 951 $ B( L1, MIN( L2+1, N ) ), LDB ) 952 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 953* 954 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 955 $ C( MIN( K2+1, M ), L2 ), 1 ) 956 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 957 $ B( L2, MIN( L2+1, N ) ), LDB ) 958 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 959* 960 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 961 $ C( MIN( K2+1, M ), L1 ), 1 ) 962 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 963 $ B( L1, MIN( L2+1, N ) ), LDB ) 964 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 965* 966 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 967 $ C( MIN( K2+1, M ), L2 ), 1 ) 968 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 969 $ B( L2, MIN( L2+1, N ) ), LDB ) 970 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 971* 972 CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ), 973 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 974 $ 2, XNORM, IERR ) 975 IF( IERR.NE.0 ) 976 $ INFO = 1 977* 978 IF( SCALOC.NE.ONE ) THEN 979 DO 220 J = 1, N 980 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 981 220 CONTINUE 982 SCALE = SCALE*SCALOC 983 END IF 984 C( K1, L1 ) = X( 1, 1 ) 985 C( K1, L2 ) = X( 1, 2 ) 986 C( K2, L1 ) = X( 2, 1 ) 987 C( K2, L2 ) = X( 2, 2 ) 988 END IF 989* 990 230 CONTINUE 991 240 CONTINUE 992* 993 END IF 994* 995 RETURN 996* 997* End of DTRSYL 998* 999 END 1000