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*> \date December 2016 160* 161*> \ingroup doubleSYcomputational 162* 163* ===================================================================== 164 SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 165 $ LDC, SCALE, INFO ) 166* 167* -- LAPACK computational routine (version 3.7.0) -- 168* -- LAPACK is a software package provided by Univ. of Tennessee, -- 169* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 170* December 2016 171* 172* .. Scalar Arguments .. 173 CHARACTER TRANA, TRANB 174 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 175 DOUBLE PRECISION SCALE 176* .. 177* .. Array Arguments .. 178 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) 179* .. 180* 181* ===================================================================== 182* 183* .. Parameters .. 184 DOUBLE PRECISION ZERO, ONE 185 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 186* .. 187* .. Local Scalars .. 188 LOGICAL NOTRNA, NOTRNB 189 INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT 190 DOUBLE PRECISION A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, 191 $ SMLNUM, SUML, SUMR, XNORM 192* .. 193* .. Local Arrays .. 194 DOUBLE PRECISION DUM( 1 ), VEC( 2, 2 ), X( 2, 2 ) 195* .. 196* .. External Functions .. 197 LOGICAL LSAME 198 DOUBLE PRECISION DDOT, DLAMCH, DLANGE 199 EXTERNAL LSAME, DDOT, DLAMCH, DLANGE 200* .. 201* .. External Subroutines .. 202 EXTERNAL DLABAD, DLALN2, DLASY2, DSCAL, XERBLA 203* .. 204* .. Intrinsic Functions .. 205 INTRINSIC ABS, DBLE, MAX, MIN 206* .. 207* .. Executable Statements .. 208* 209* Decode and Test input parameters 210* 211 NOTRNA = LSAME( TRANA, 'N' ) 212 NOTRNB = LSAME( TRANB, 'N' ) 213* 214 INFO = 0 215 IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT. 216 $ LSAME( TRANA, 'C' ) ) THEN 217 INFO = -1 218 ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT. 219 $ LSAME( TRANB, 'C' ) ) THEN 220 INFO = -2 221 ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN 222 INFO = -3 223 ELSE IF( M.LT.0 ) THEN 224 INFO = -4 225 ELSE IF( N.LT.0 ) THEN 226 INFO = -5 227 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 228 INFO = -7 229 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 230 INFO = -9 231 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 232 INFO = -11 233 END IF 234 IF( INFO.NE.0 ) THEN 235 CALL XERBLA( 'DTRSYL', -INFO ) 236 RETURN 237 END IF 238* 239* Quick return if possible 240* 241 SCALE = ONE 242 IF( M.EQ.0 .OR. N.EQ.0 ) 243 $ RETURN 244* 245* Set constants to control overflow 246* 247 EPS = DLAMCH( 'P' ) 248 SMLNUM = DLAMCH( 'S' ) 249 BIGNUM = ONE / SMLNUM 250 CALL DLABAD( SMLNUM, BIGNUM ) 251 SMLNUM = SMLNUM*DBLE( M*N ) / EPS 252 BIGNUM = ONE / SMLNUM 253* 254 SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ), 255 $ EPS*DLANGE( 'M', N, N, B, LDB, DUM ) ) 256* 257 SGN = ISGN 258* 259 IF( NOTRNA .AND. NOTRNB ) THEN 260* 261* Solve A*X + ISGN*X*B = scale*C. 262* 263* The (K,L)th block of X is determined starting from 264* bottom-left corner column by column by 265* 266* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 267* 268* Where 269* M L-1 270* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]. 271* I=K+1 J=1 272* 273* Start column loop (index = L) 274* L1 (L2) : column index of the first (first) row of X(K,L). 275* 276 LNEXT = 1 277 DO 60 L = 1, N 278 IF( L.LT.LNEXT ) 279 $ GO TO 60 280 IF( L.EQ.N ) THEN 281 L1 = L 282 L2 = L 283 ELSE 284 IF( B( L+1, L ).NE.ZERO ) THEN 285 L1 = L 286 L2 = L + 1 287 LNEXT = L + 2 288 ELSE 289 L1 = L 290 L2 = L 291 LNEXT = L + 1 292 END IF 293 END IF 294* 295* Start row loop (index = K) 296* K1 (K2): row index of the first (last) row of X(K,L). 297* 298 KNEXT = M 299 DO 50 K = M, 1, -1 300 IF( K.GT.KNEXT ) 301 $ GO TO 50 302 IF( K.EQ.1 ) THEN 303 K1 = K 304 K2 = K 305 ELSE 306 IF( A( K, K-1 ).NE.ZERO ) THEN 307 K1 = K - 1 308 K2 = K 309 KNEXT = K - 2 310 ELSE 311 K1 = K 312 K2 = K 313 KNEXT = K - 1 314 END IF 315 END IF 316* 317 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 318 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 319 $ C( MIN( K1+1, M ), L1 ), 1 ) 320 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 321 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 322 SCALOC = ONE 323* 324 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 325 DA11 = ABS( A11 ) 326 IF( DA11.LE.SMIN ) THEN 327 A11 = SMIN 328 DA11 = SMIN 329 INFO = 1 330 END IF 331 DB = ABS( VEC( 1, 1 ) ) 332 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 333 IF( DB.GT.BIGNUM*DA11 ) 334 $ SCALOC = ONE / DB 335 END IF 336 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 337* 338 IF( SCALOC.NE.ONE ) THEN 339 DO 10 J = 1, N 340 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 341 10 CONTINUE 342 SCALE = SCALE*SCALOC 343 END IF 344 C( K1, L1 ) = X( 1, 1 ) 345* 346 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 347* 348 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 349 $ C( MIN( K2+1, M ), L1 ), 1 ) 350 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 351 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 352* 353 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 354 $ C( MIN( K2+1, M ), L1 ), 1 ) 355 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 356 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 357* 358 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ), 359 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 360 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 361 IF( IERR.NE.0 ) 362 $ INFO = 1 363* 364 IF( SCALOC.NE.ONE ) THEN 365 DO 20 J = 1, N 366 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 367 20 CONTINUE 368 SCALE = SCALE*SCALOC 369 END IF 370 C( K1, L1 ) = X( 1, 1 ) 371 C( K2, L1 ) = X( 2, 1 ) 372* 373 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 374* 375 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 376 $ C( MIN( K1+1, M ), L1 ), 1 ) 377 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 378 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 379* 380 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 381 $ C( MIN( K1+1, M ), L2 ), 1 ) 382 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 383 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 384* 385 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ), 386 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 387 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 388 IF( IERR.NE.0 ) 389 $ INFO = 1 390* 391 IF( SCALOC.NE.ONE ) THEN 392 DO 30 J = 1, N 393 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 394 30 CONTINUE 395 SCALE = SCALE*SCALOC 396 END IF 397 C( K1, L1 ) = X( 1, 1 ) 398 C( K1, L2 ) = X( 2, 1 ) 399* 400 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 401* 402 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 403 $ C( MIN( K2+1, M ), L1 ), 1 ) 404 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 405 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 406* 407 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 408 $ C( MIN( K2+1, M ), L2 ), 1 ) 409 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 410 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 411* 412 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 413 $ C( MIN( K2+1, M ), L1 ), 1 ) 414 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 415 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 416* 417 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 418 $ C( MIN( K2+1, M ), L2 ), 1 ) 419 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 420 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 421* 422 CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2, 423 $ A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC, 424 $ 2, SCALOC, X, 2, XNORM, IERR ) 425 IF( IERR.NE.0 ) 426 $ INFO = 1 427* 428 IF( SCALOC.NE.ONE ) THEN 429 DO 40 J = 1, N 430 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 431 40 CONTINUE 432 SCALE = SCALE*SCALOC 433 END IF 434 C( K1, L1 ) = X( 1, 1 ) 435 C( K1, L2 ) = X( 1, 2 ) 436 C( K2, L1 ) = X( 2, 1 ) 437 C( K2, L2 ) = X( 2, 2 ) 438 END IF 439* 440 50 CONTINUE 441* 442 60 CONTINUE 443* 444 ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN 445* 446* Solve A**T *X + ISGN*X*B = scale*C. 447* 448* The (K,L)th block of X is determined starting from 449* upper-left corner column by column by 450* 451* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 452* 453* Where 454* K-1 T L-1 455* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)] 456* I=1 J=1 457* 458* Start column loop (index = L) 459* L1 (L2): column index of the first (last) row of X(K,L) 460* 461 LNEXT = 1 462 DO 120 L = 1, N 463 IF( L.LT.LNEXT ) 464 $ GO TO 120 465 IF( L.EQ.N ) THEN 466 L1 = L 467 L2 = L 468 ELSE 469 IF( B( L+1, L ).NE.ZERO ) THEN 470 L1 = L 471 L2 = L + 1 472 LNEXT = L + 2 473 ELSE 474 L1 = L 475 L2 = L 476 LNEXT = L + 1 477 END IF 478 END IF 479* 480* Start row loop (index = K) 481* K1 (K2): row index of the first (last) row of X(K,L) 482* 483 KNEXT = 1 484 DO 110 K = 1, M 485 IF( K.LT.KNEXT ) 486 $ GO TO 110 487 IF( K.EQ.M ) THEN 488 K1 = K 489 K2 = K 490 ELSE 491 IF( A( K+1, K ).NE.ZERO ) THEN 492 K1 = K 493 K2 = K + 1 494 KNEXT = K + 2 495 ELSE 496 K1 = K 497 K2 = K 498 KNEXT = K + 1 499 END IF 500 END IF 501* 502 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 503 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 504 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 505 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 506 SCALOC = ONE 507* 508 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 509 DA11 = ABS( A11 ) 510 IF( DA11.LE.SMIN ) THEN 511 A11 = SMIN 512 DA11 = SMIN 513 INFO = 1 514 END IF 515 DB = ABS( VEC( 1, 1 ) ) 516 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 517 IF( DB.GT.BIGNUM*DA11 ) 518 $ SCALOC = ONE / DB 519 END IF 520 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 521* 522 IF( SCALOC.NE.ONE ) THEN 523 DO 70 J = 1, N 524 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 525 70 CONTINUE 526 SCALE = SCALE*SCALOC 527 END IF 528 C( K1, L1 ) = X( 1, 1 ) 529* 530 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 531* 532 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 533 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 534 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 535* 536 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 537 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 538 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 539* 540 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ), 541 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 542 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 543 IF( IERR.NE.0 ) 544 $ INFO = 1 545* 546 IF( SCALOC.NE.ONE ) THEN 547 DO 80 J = 1, N 548 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 549 80 CONTINUE 550 SCALE = SCALE*SCALOC 551 END IF 552 C( K1, L1 ) = X( 1, 1 ) 553 C( K2, L1 ) = X( 2, 1 ) 554* 555 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 556* 557 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 558 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 559 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 560* 561 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 562 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 563 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 564* 565 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ), 566 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 567 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 568 IF( IERR.NE.0 ) 569 $ INFO = 1 570* 571 IF( SCALOC.NE.ONE ) THEN 572 DO 90 J = 1, N 573 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 574 90 CONTINUE 575 SCALE = SCALE*SCALOC 576 END IF 577 C( K1, L1 ) = X( 1, 1 ) 578 C( K1, L2 ) = X( 2, 1 ) 579* 580 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 581* 582 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 583 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 584 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 585* 586 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 587 SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 588 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 589* 590 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 591 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 592 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 593* 594 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 595 SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 596 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 597* 598 CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ), 599 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 600 $ 2, XNORM, IERR ) 601 IF( IERR.NE.0 ) 602 $ INFO = 1 603* 604 IF( SCALOC.NE.ONE ) THEN 605 DO 100 J = 1, N 606 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 607 100 CONTINUE 608 SCALE = SCALE*SCALOC 609 END IF 610 C( K1, L1 ) = X( 1, 1 ) 611 C( K1, L2 ) = X( 1, 2 ) 612 C( K2, L1 ) = X( 2, 1 ) 613 C( K2, L2 ) = X( 2, 2 ) 614 END IF 615* 616 110 CONTINUE 617 120 CONTINUE 618* 619 ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN 620* 621* Solve A**T*X + ISGN*X*B**T = scale*C. 622* 623* The (K,L)th block of X is determined starting from 624* top-right corner column by column by 625* 626* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) 627* 628* Where 629* K-1 N 630* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. 631* I=1 J=L+1 632* 633* Start column loop (index = L) 634* L1 (L2): column index of the first (last) row of X(K,L) 635* 636 LNEXT = N 637 DO 180 L = N, 1, -1 638 IF( L.GT.LNEXT ) 639 $ GO TO 180 640 IF( L.EQ.1 ) THEN 641 L1 = L 642 L2 = L 643 ELSE 644 IF( B( L, L-1 ).NE.ZERO ) THEN 645 L1 = L - 1 646 L2 = L 647 LNEXT = L - 2 648 ELSE 649 L1 = L 650 L2 = L 651 LNEXT = L - 1 652 END IF 653 END IF 654* 655* Start row loop (index = K) 656* K1 (K2): row index of the first (last) row of X(K,L) 657* 658 KNEXT = 1 659 DO 170 K = 1, M 660 IF( K.LT.KNEXT ) 661 $ GO TO 170 662 IF( K.EQ.M ) THEN 663 K1 = K 664 K2 = K 665 ELSE 666 IF( A( K+1, K ).NE.ZERO ) THEN 667 K1 = K 668 K2 = K + 1 669 KNEXT = K + 2 670 ELSE 671 K1 = K 672 K2 = K 673 KNEXT = K + 1 674 END IF 675 END IF 676* 677 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 678 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 679 SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC, 680 $ B( L1, MIN( L1+1, N ) ), LDB ) 681 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 682 SCALOC = ONE 683* 684 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 685 DA11 = ABS( A11 ) 686 IF( DA11.LE.SMIN ) THEN 687 A11 = SMIN 688 DA11 = SMIN 689 INFO = 1 690 END IF 691 DB = ABS( VEC( 1, 1 ) ) 692 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 693 IF( DB.GT.BIGNUM*DA11 ) 694 $ SCALOC = ONE / DB 695 END IF 696 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 697* 698 IF( SCALOC.NE.ONE ) THEN 699 DO 130 J = 1, N 700 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 701 130 CONTINUE 702 SCALE = SCALE*SCALOC 703 END IF 704 C( K1, L1 ) = X( 1, 1 ) 705* 706 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 707* 708 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 709 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 710 $ B( L1, MIN( L2+1, N ) ), LDB ) 711 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 712* 713 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 714 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 715 $ B( L1, MIN( L2+1, N ) ), LDB ) 716 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 717* 718 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ), 719 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 720 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 721 IF( IERR.NE.0 ) 722 $ INFO = 1 723* 724 IF( SCALOC.NE.ONE ) THEN 725 DO 140 J = 1, N 726 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 727 140 CONTINUE 728 SCALE = SCALE*SCALOC 729 END IF 730 C( K1, L1 ) = X( 1, 1 ) 731 C( K2, L1 ) = X( 2, 1 ) 732* 733 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 734* 735 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 736 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 737 $ B( L1, MIN( L2+1, N ) ), LDB ) 738 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 739* 740 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 741 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 742 $ B( L2, MIN( L2+1, N ) ), LDB ) 743 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 744* 745 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ), 746 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 747 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 748 IF( IERR.NE.0 ) 749 $ INFO = 1 750* 751 IF( SCALOC.NE.ONE ) THEN 752 DO 150 J = 1, N 753 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 754 150 CONTINUE 755 SCALE = SCALE*SCALOC 756 END IF 757 C( K1, L1 ) = X( 1, 1 ) 758 C( K1, L2 ) = X( 2, 1 ) 759* 760 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 761* 762 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 763 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 764 $ B( L1, MIN( L2+1, N ) ), LDB ) 765 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 766* 767 SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 768 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 769 $ B( L2, MIN( L2+1, N ) ), LDB ) 770 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 771* 772 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 773 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 774 $ B( L1, MIN( L2+1, N ) ), LDB ) 775 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 776* 777 SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 778 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 779 $ B( L2, MIN( L2+1, N ) ), LDB ) 780 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 781* 782 CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ), 783 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 784 $ 2, XNORM, IERR ) 785 IF( IERR.NE.0 ) 786 $ INFO = 1 787* 788 IF( SCALOC.NE.ONE ) THEN 789 DO 160 J = 1, N 790 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 791 160 CONTINUE 792 SCALE = SCALE*SCALOC 793 END IF 794 C( K1, L1 ) = X( 1, 1 ) 795 C( K1, L2 ) = X( 1, 2 ) 796 C( K2, L1 ) = X( 2, 1 ) 797 C( K2, L2 ) = X( 2, 2 ) 798 END IF 799* 800 170 CONTINUE 801 180 CONTINUE 802* 803 ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN 804* 805* Solve A*X + ISGN*X*B**T = scale*C. 806* 807* The (K,L)th block of X is determined starting from 808* bottom-right corner column by column by 809* 810* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) 811* 812* Where 813* M N 814* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. 815* I=K+1 J=L+1 816* 817* Start column loop (index = L) 818* L1 (L2): column index of the first (last) row of X(K,L) 819* 820 LNEXT = N 821 DO 240 L = N, 1, -1 822 IF( L.GT.LNEXT ) 823 $ GO TO 240 824 IF( L.EQ.1 ) THEN 825 L1 = L 826 L2 = L 827 ELSE 828 IF( B( L, L-1 ).NE.ZERO ) THEN 829 L1 = L - 1 830 L2 = L 831 LNEXT = L - 2 832 ELSE 833 L1 = L 834 L2 = L 835 LNEXT = L - 1 836 END IF 837 END IF 838* 839* Start row loop (index = K) 840* K1 (K2): row index of the first (last) row of X(K,L) 841* 842 KNEXT = M 843 DO 230 K = M, 1, -1 844 IF( K.GT.KNEXT ) 845 $ GO TO 230 846 IF( K.EQ.1 ) THEN 847 K1 = K 848 K2 = K 849 ELSE 850 IF( A( K, K-1 ).NE.ZERO ) THEN 851 K1 = K - 1 852 K2 = K 853 KNEXT = K - 2 854 ELSE 855 K1 = K 856 K2 = K 857 KNEXT = K - 1 858 END IF 859 END IF 860* 861 IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 862 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 863 $ C( MIN( K1+1, M ), L1 ), 1 ) 864 SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC, 865 $ B( L1, MIN( L1+1, N ) ), LDB ) 866 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 867 SCALOC = ONE 868* 869 A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 870 DA11 = ABS( A11 ) 871 IF( DA11.LE.SMIN ) THEN 872 A11 = SMIN 873 DA11 = SMIN 874 INFO = 1 875 END IF 876 DB = ABS( VEC( 1, 1 ) ) 877 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 878 IF( DB.GT.BIGNUM*DA11 ) 879 $ SCALOC = ONE / DB 880 END IF 881 X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 882* 883 IF( SCALOC.NE.ONE ) THEN 884 DO 190 J = 1, N 885 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 886 190 CONTINUE 887 SCALE = SCALE*SCALOC 888 END IF 889 C( K1, L1 ) = X( 1, 1 ) 890* 891 ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 892* 893 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 894 $ C( MIN( K2+1, M ), L1 ), 1 ) 895 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 896 $ B( L1, MIN( L2+1, N ) ), LDB ) 897 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 898* 899 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 900 $ C( MIN( K2+1, M ), L1 ), 1 ) 901 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 902 $ B( L1, MIN( L2+1, N ) ), LDB ) 903 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 904* 905 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ), 906 $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 907 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 908 IF( IERR.NE.0 ) 909 $ INFO = 1 910* 911 IF( SCALOC.NE.ONE ) THEN 912 DO 200 J = 1, N 913 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 914 200 CONTINUE 915 SCALE = SCALE*SCALOC 916 END IF 917 C( K1, L1 ) = X( 1, 1 ) 918 C( K2, L1 ) = X( 2, 1 ) 919* 920 ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 921* 922 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 923 $ C( MIN( K1+1, M ), L1 ), 1 ) 924 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 925 $ B( L1, MIN( L2+1, N ) ), LDB ) 926 VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 927* 928 SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 929 $ C( MIN( K1+1, M ), L2 ), 1 ) 930 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 931 $ B( L2, MIN( L2+1, N ) ), LDB ) 932 VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 933* 934 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ), 935 $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 936 $ ZERO, X, 2, SCALOC, XNORM, IERR ) 937 IF( IERR.NE.0 ) 938 $ INFO = 1 939* 940 IF( SCALOC.NE.ONE ) THEN 941 DO 210 J = 1, N 942 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 943 210 CONTINUE 944 SCALE = SCALE*SCALOC 945 END IF 946 C( K1, L1 ) = X( 1, 1 ) 947 C( K1, L2 ) = X( 2, 1 ) 948* 949 ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 950* 951 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 952 $ C( MIN( K2+1, M ), L1 ), 1 ) 953 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 954 $ B( L1, MIN( L2+1, N ) ), LDB ) 955 VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 956* 957 SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 958 $ C( MIN( K2+1, M ), L2 ), 1 ) 959 SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 960 $ B( L2, MIN( L2+1, N ) ), LDB ) 961 VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 962* 963 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 964 $ C( MIN( K2+1, M ), L1 ), 1 ) 965 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 966 $ B( L1, MIN( L2+1, N ) ), LDB ) 967 VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 968* 969 SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 970 $ C( MIN( K2+1, M ), L2 ), 1 ) 971 SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 972 $ B( L2, MIN( L2+1, N ) ), LDB ) 973 VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 974* 975 CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ), 976 $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 977 $ 2, XNORM, IERR ) 978 IF( IERR.NE.0 ) 979 $ INFO = 1 980* 981 IF( SCALOC.NE.ONE ) THEN 982 DO 220 J = 1, N 983 CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 984 220 CONTINUE 985 SCALE = SCALE*SCALOC 986 END IF 987 C( K1, L1 ) = X( 1, 1 ) 988 C( K1, L2 ) = X( 1, 2 ) 989 C( K2, L1 ) = X( 2, 1 ) 990 C( K2, L2 ) = X( 2, 2 ) 991 END IF 992* 993 230 CONTINUE 994 240 CONTINUE 995* 996 END IF 997* 998 RETURN 999* 1000* End of DTRSYL 1001* 1002 END 1003