1*> \brief \b DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of special form, in real arithmetic. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLAQTR + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqtr.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqtr.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqtr.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, 22* INFO ) 23* 24* .. Scalar Arguments .. 25* LOGICAL LREAL, LTRAN 26* INTEGER INFO, LDT, N 27* DOUBLE PRECISION SCALE, W 28* .. 29* .. Array Arguments .. 30* DOUBLE PRECISION B( * ), T( LDT, * ), WORK( * ), X( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DLAQTR solves the real quasi-triangular system 40*> 41*> op(T)*p = scale*c, if LREAL = .TRUE. 42*> 43*> or the complex quasi-triangular systems 44*> 45*> op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE. 46*> 47*> in real arithmetic, where T is upper quasi-triangular. 48*> If LREAL = .FALSE., then the first diagonal block of T must be 49*> 1 by 1, B is the specially structured matrix 50*> 51*> B = [ b(1) b(2) ... b(n) ] 52*> [ w ] 53*> [ w ] 54*> [ . ] 55*> [ w ] 56*> 57*> op(A) = A or A**T, A**T denotes the transpose of 58*> matrix A. 59*> 60*> On input, X = [ c ]. On output, X = [ p ]. 61*> [ d ] [ q ] 62*> 63*> This subroutine is designed for the condition number estimation 64*> in routine DTRSNA. 65*> \endverbatim 66* 67* Arguments: 68* ========== 69* 70*> \param[in] LTRAN 71*> \verbatim 72*> LTRAN is LOGICAL 73*> On entry, LTRAN specifies the option of conjugate transpose: 74*> = .FALSE., op(T+i*B) = T+i*B, 75*> = .TRUE., op(T+i*B) = (T+i*B)**T. 76*> \endverbatim 77*> 78*> \param[in] LREAL 79*> \verbatim 80*> LREAL is LOGICAL 81*> On entry, LREAL specifies the input matrix structure: 82*> = .FALSE., the input is complex 83*> = .TRUE., the input is real 84*> \endverbatim 85*> 86*> \param[in] N 87*> \verbatim 88*> N is INTEGER 89*> On entry, N specifies the order of T+i*B. N >= 0. 90*> \endverbatim 91*> 92*> \param[in] T 93*> \verbatim 94*> T is DOUBLE PRECISION array, dimension (LDT,N) 95*> On entry, T contains a matrix in Schur canonical form. 96*> If LREAL = .FALSE., then the first diagonal block of T mu 97*> be 1 by 1. 98*> \endverbatim 99*> 100*> \param[in] LDT 101*> \verbatim 102*> LDT is INTEGER 103*> The leading dimension of the matrix T. LDT >= max(1,N). 104*> \endverbatim 105*> 106*> \param[in] B 107*> \verbatim 108*> B is DOUBLE PRECISION array, dimension (N) 109*> On entry, B contains the elements to form the matrix 110*> B as described above. 111*> If LREAL = .TRUE., B is not referenced. 112*> \endverbatim 113*> 114*> \param[in] W 115*> \verbatim 116*> W is DOUBLE PRECISION 117*> On entry, W is the diagonal element of the matrix B. 118*> If LREAL = .TRUE., W is not referenced. 119*> \endverbatim 120*> 121*> \param[out] SCALE 122*> \verbatim 123*> SCALE is DOUBLE PRECISION 124*> On exit, SCALE is the scale factor. 125*> \endverbatim 126*> 127*> \param[in,out] X 128*> \verbatim 129*> X is DOUBLE PRECISION array, dimension (2*N) 130*> On entry, X contains the right hand side of the system. 131*> On exit, X is overwritten by the solution. 132*> \endverbatim 133*> 134*> \param[out] WORK 135*> \verbatim 136*> WORK is DOUBLE PRECISION array, dimension (N) 137*> \endverbatim 138*> 139*> \param[out] INFO 140*> \verbatim 141*> INFO is INTEGER 142*> On exit, INFO is set to 143*> 0: successful exit. 144*> 1: the some diagonal 1 by 1 block has been perturbed by 145*> a small number SMIN to keep nonsingularity. 146*> 2: the some diagonal 2 by 2 block has been perturbed by 147*> a small number in DLALN2 to keep nonsingularity. 148*> NOTE: In the interests of speed, this routine does not 149*> check the inputs for errors. 150*> \endverbatim 151* 152* Authors: 153* ======== 154* 155*> \author Univ. of Tennessee 156*> \author Univ. of California Berkeley 157*> \author Univ. of Colorado Denver 158*> \author NAG Ltd. 159* 160*> \ingroup doubleOTHERauxiliary 161* 162* ===================================================================== 163 SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, 164 $ INFO ) 165* 166* -- LAPACK auxiliary routine -- 167* -- LAPACK is a software package provided by Univ. of Tennessee, -- 168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 169* 170* .. Scalar Arguments .. 171 LOGICAL LREAL, LTRAN 172 INTEGER INFO, LDT, N 173 DOUBLE PRECISION SCALE, W 174* .. 175* .. Array Arguments .. 176 DOUBLE PRECISION B( * ), T( LDT, * ), WORK( * ), X( * ) 177* .. 178* 179* ===================================================================== 180* 181* .. Parameters .. 182 DOUBLE PRECISION ZERO, ONE 183 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 184* .. 185* .. Local Scalars .. 186 LOGICAL NOTRAN 187 INTEGER I, IERR, J, J1, J2, JNEXT, K, N1, N2 188 DOUBLE PRECISION BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW, 189 $ SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z 190* .. 191* .. Local Arrays .. 192 DOUBLE PRECISION D( 2, 2 ), V( 2, 2 ) 193* .. 194* .. External Functions .. 195 INTEGER IDAMAX 196 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE 197 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE 198* .. 199* .. External Subroutines .. 200 EXTERNAL DAXPY, DLADIV, DLALN2, DSCAL 201* .. 202* .. Intrinsic Functions .. 203 INTRINSIC ABS, MAX 204* .. 205* .. Executable Statements .. 206* 207* Do not test the input parameters for errors 208* 209 NOTRAN = .NOT.LTRAN 210 INFO = 0 211* 212* Quick return if possible 213* 214 IF( N.EQ.0 ) 215 $ RETURN 216* 217* Set constants to control overflow 218* 219 EPS = DLAMCH( 'P' ) 220 SMLNUM = DLAMCH( 'S' ) / EPS 221 BIGNUM = ONE / SMLNUM 222* 223 XNORM = DLANGE( 'M', N, N, T, LDT, D ) 224 IF( .NOT.LREAL ) 225 $ XNORM = MAX( XNORM, ABS( W ), DLANGE( 'M', N, 1, B, N, D ) ) 226 SMIN = MAX( SMLNUM, EPS*XNORM ) 227* 228* Compute 1-norm of each column of strictly upper triangular 229* part of T to control overflow in triangular solver. 230* 231 WORK( 1 ) = ZERO 232 DO 10 J = 2, N 233 WORK( J ) = DASUM( J-1, T( 1, J ), 1 ) 234 10 CONTINUE 235* 236 IF( .NOT.LREAL ) THEN 237 DO 20 I = 2, N 238 WORK( I ) = WORK( I ) + ABS( B( I ) ) 239 20 CONTINUE 240 END IF 241* 242 N2 = 2*N 243 N1 = N 244 IF( .NOT.LREAL ) 245 $ N1 = N2 246 K = IDAMAX( N1, X, 1 ) 247 XMAX = ABS( X( K ) ) 248 SCALE = ONE 249* 250 IF( XMAX.GT.BIGNUM ) THEN 251 SCALE = BIGNUM / XMAX 252 CALL DSCAL( N1, SCALE, X, 1 ) 253 XMAX = BIGNUM 254 END IF 255* 256 IF( LREAL ) THEN 257* 258 IF( NOTRAN ) THEN 259* 260* Solve T*p = scale*c 261* 262 JNEXT = N 263 DO 30 J = N, 1, -1 264 IF( J.GT.JNEXT ) 265 $ GO TO 30 266 J1 = J 267 J2 = J 268 JNEXT = J - 1 269 IF( J.GT.1 ) THEN 270 IF( T( J, J-1 ).NE.ZERO ) THEN 271 J1 = J - 1 272 JNEXT = J - 2 273 END IF 274 END IF 275* 276 IF( J1.EQ.J2 ) THEN 277* 278* Meet 1 by 1 diagonal block 279* 280* Scale to avoid overflow when computing 281* x(j) = b(j)/T(j,j) 282* 283 XJ = ABS( X( J1 ) ) 284 TJJ = ABS( T( J1, J1 ) ) 285 TMP = T( J1, J1 ) 286 IF( TJJ.LT.SMIN ) THEN 287 TMP = SMIN 288 TJJ = SMIN 289 INFO = 1 290 END IF 291* 292 IF( XJ.EQ.ZERO ) 293 $ GO TO 30 294* 295 IF( TJJ.LT.ONE ) THEN 296 IF( XJ.GT.BIGNUM*TJJ ) THEN 297 REC = ONE / XJ 298 CALL DSCAL( N, REC, X, 1 ) 299 SCALE = SCALE*REC 300 XMAX = XMAX*REC 301 END IF 302 END IF 303 X( J1 ) = X( J1 ) / TMP 304 XJ = ABS( X( J1 ) ) 305* 306* Scale x if necessary to avoid overflow when adding a 307* multiple of column j1 of T. 308* 309 IF( XJ.GT.ONE ) THEN 310 REC = ONE / XJ 311 IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN 312 CALL DSCAL( N, REC, X, 1 ) 313 SCALE = SCALE*REC 314 END IF 315 END IF 316 IF( J1.GT.1 ) THEN 317 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 ) 318 K = IDAMAX( J1-1, X, 1 ) 319 XMAX = ABS( X( K ) ) 320 END IF 321* 322 ELSE 323* 324* Meet 2 by 2 diagonal block 325* 326* Call 2 by 2 linear system solve, to take 327* care of possible overflow by scaling factor. 328* 329 D( 1, 1 ) = X( J1 ) 330 D( 2, 1 ) = X( J2 ) 331 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ), 332 $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2, 333 $ SCALOC, XNORM, IERR ) 334 IF( IERR.NE.0 ) 335 $ INFO = 2 336* 337 IF( SCALOC.NE.ONE ) THEN 338 CALL DSCAL( N, SCALOC, X, 1 ) 339 SCALE = SCALE*SCALOC 340 END IF 341 X( J1 ) = V( 1, 1 ) 342 X( J2 ) = V( 2, 1 ) 343* 344* Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2)) 345* to avoid overflow in updating right-hand side. 346* 347 XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) ) 348 IF( XJ.GT.ONE ) THEN 349 REC = ONE / XJ 350 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT. 351 $ ( BIGNUM-XMAX )*REC ) THEN 352 CALL DSCAL( N, REC, X, 1 ) 353 SCALE = SCALE*REC 354 END IF 355 END IF 356* 357* Update right-hand side 358* 359 IF( J1.GT.1 ) THEN 360 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 ) 361 CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 ) 362 K = IDAMAX( J1-1, X, 1 ) 363 XMAX = ABS( X( K ) ) 364 END IF 365* 366 END IF 367* 368 30 CONTINUE 369* 370 ELSE 371* 372* Solve T**T*p = scale*c 373* 374 JNEXT = 1 375 DO 40 J = 1, N 376 IF( J.LT.JNEXT ) 377 $ GO TO 40 378 J1 = J 379 J2 = J 380 JNEXT = J + 1 381 IF( J.LT.N ) THEN 382 IF( T( J+1, J ).NE.ZERO ) THEN 383 J2 = J + 1 384 JNEXT = J + 2 385 END IF 386 END IF 387* 388 IF( J1.EQ.J2 ) THEN 389* 390* 1 by 1 diagonal block 391* 392* Scale if necessary to avoid overflow in forming the 393* right-hand side element by inner product. 394* 395 XJ = ABS( X( J1 ) ) 396 IF( XMAX.GT.ONE ) THEN 397 REC = ONE / XMAX 398 IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN 399 CALL DSCAL( N, REC, X, 1 ) 400 SCALE = SCALE*REC 401 XMAX = XMAX*REC 402 END IF 403 END IF 404* 405 X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 ) 406* 407 XJ = ABS( X( J1 ) ) 408 TJJ = ABS( T( J1, J1 ) ) 409 TMP = T( J1, J1 ) 410 IF( TJJ.LT.SMIN ) THEN 411 TMP = SMIN 412 TJJ = SMIN 413 INFO = 1 414 END IF 415* 416 IF( TJJ.LT.ONE ) THEN 417 IF( XJ.GT.BIGNUM*TJJ ) THEN 418 REC = ONE / XJ 419 CALL DSCAL( N, REC, X, 1 ) 420 SCALE = SCALE*REC 421 XMAX = XMAX*REC 422 END IF 423 END IF 424 X( J1 ) = X( J1 ) / TMP 425 XMAX = MAX( XMAX, ABS( X( J1 ) ) ) 426* 427 ELSE 428* 429* 2 by 2 diagonal block 430* 431* Scale if necessary to avoid overflow in forming the 432* right-hand side elements by inner product. 433* 434 XJ = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ) ) 435 IF( XMAX.GT.ONE ) THEN 436 REC = ONE / XMAX 437 IF( MAX( WORK( J2 ), WORK( J1 ) ).GT.( BIGNUM-XJ )* 438 $ REC ) THEN 439 CALL DSCAL( N, REC, X, 1 ) 440 SCALE = SCALE*REC 441 XMAX = XMAX*REC 442 END IF 443 END IF 444* 445 D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 446 $ 1 ) 447 D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X, 448 $ 1 ) 449* 450 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J1, J1 ), 451 $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2, 452 $ SCALOC, XNORM, IERR ) 453 IF( IERR.NE.0 ) 454 $ INFO = 2 455* 456 IF( SCALOC.NE.ONE ) THEN 457 CALL DSCAL( N, SCALOC, X, 1 ) 458 SCALE = SCALE*SCALOC 459 END IF 460 X( J1 ) = V( 1, 1 ) 461 X( J2 ) = V( 2, 1 ) 462 XMAX = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ), XMAX ) 463* 464 END IF 465 40 CONTINUE 466 END IF 467* 468 ELSE 469* 470 SMINW = MAX( EPS*ABS( W ), SMIN ) 471 IF( NOTRAN ) THEN 472* 473* Solve (T + iB)*(p+iq) = c+id 474* 475 JNEXT = N 476 DO 70 J = N, 1, -1 477 IF( J.GT.JNEXT ) 478 $ GO TO 70 479 J1 = J 480 J2 = J 481 JNEXT = J - 1 482 IF( J.GT.1 ) THEN 483 IF( T( J, J-1 ).NE.ZERO ) THEN 484 J1 = J - 1 485 JNEXT = J - 2 486 END IF 487 END IF 488* 489 IF( J1.EQ.J2 ) THEN 490* 491* 1 by 1 diagonal block 492* 493* Scale if necessary to avoid overflow in division 494* 495 Z = W 496 IF( J1.EQ.1 ) 497 $ Z = B( 1 ) 498 XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) ) 499 TJJ = ABS( T( J1, J1 ) ) + ABS( Z ) 500 TMP = T( J1, J1 ) 501 IF( TJJ.LT.SMINW ) THEN 502 TMP = SMINW 503 TJJ = SMINW 504 INFO = 1 505 END IF 506* 507 IF( XJ.EQ.ZERO ) 508 $ GO TO 70 509* 510 IF( TJJ.LT.ONE ) THEN 511 IF( XJ.GT.BIGNUM*TJJ ) THEN 512 REC = ONE / XJ 513 CALL DSCAL( N2, REC, X, 1 ) 514 SCALE = SCALE*REC 515 XMAX = XMAX*REC 516 END IF 517 END IF 518 CALL DLADIV( X( J1 ), X( N+J1 ), TMP, Z, SR, SI ) 519 X( J1 ) = SR 520 X( N+J1 ) = SI 521 XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) ) 522* 523* Scale x if necessary to avoid overflow when adding a 524* multiple of column j1 of T. 525* 526 IF( XJ.GT.ONE ) THEN 527 REC = ONE / XJ 528 IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN 529 CALL DSCAL( N2, REC, X, 1 ) 530 SCALE = SCALE*REC 531 END IF 532 END IF 533* 534 IF( J1.GT.1 ) THEN 535 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 ) 536 CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1, 537 $ X( N+1 ), 1 ) 538* 539 X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) 540 X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) 541* 542 XMAX = ZERO 543 DO 50 K = 1, J1 - 1 544 XMAX = MAX( XMAX, ABS( X( K ) )+ 545 $ ABS( X( K+N ) ) ) 546 50 CONTINUE 547 END IF 548* 549 ELSE 550* 551* Meet 2 by 2 diagonal block 552* 553 D( 1, 1 ) = X( J1 ) 554 D( 2, 1 ) = X( J2 ) 555 D( 1, 2 ) = X( N+J1 ) 556 D( 2, 2 ) = X( N+J2 ) 557 CALL DLALN2( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ), 558 $ LDT, ONE, ONE, D, 2, ZERO, -W, V, 2, 559 $ SCALOC, XNORM, IERR ) 560 IF( IERR.NE.0 ) 561 $ INFO = 2 562* 563 IF( SCALOC.NE.ONE ) THEN 564 CALL DSCAL( 2*N, SCALOC, X, 1 ) 565 SCALE = SCALOC*SCALE 566 END IF 567 X( J1 ) = V( 1, 1 ) 568 X( J2 ) = V( 2, 1 ) 569 X( N+J1 ) = V( 1, 2 ) 570 X( N+J2 ) = V( 2, 2 ) 571* 572* Scale X(J1), .... to avoid overflow in 573* updating right hand side. 574* 575 XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ), 576 $ ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) ) 577 IF( XJ.GT.ONE ) THEN 578 REC = ONE / XJ 579 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT. 580 $ ( BIGNUM-XMAX )*REC ) THEN 581 CALL DSCAL( N2, REC, X, 1 ) 582 SCALE = SCALE*REC 583 END IF 584 END IF 585* 586* Update the right-hand side. 587* 588 IF( J1.GT.1 ) THEN 589 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 ) 590 CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 ) 591* 592 CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1, 593 $ X( N+1 ), 1 ) 594 CALL DAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1, 595 $ X( N+1 ), 1 ) 596* 597 X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) + 598 $ B( J2 )*X( N+J2 ) 599 X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) - 600 $ B( J2 )*X( J2 ) 601* 602 XMAX = ZERO 603 DO 60 K = 1, J1 - 1 604 XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ), 605 $ XMAX ) 606 60 CONTINUE 607 END IF 608* 609 END IF 610 70 CONTINUE 611* 612 ELSE 613* 614* Solve (T + iB)**T*(p+iq) = c+id 615* 616 JNEXT = 1 617 DO 80 J = 1, N 618 IF( J.LT.JNEXT ) 619 $ GO TO 80 620 J1 = J 621 J2 = J 622 JNEXT = J + 1 623 IF( J.LT.N ) THEN 624 IF( T( J+1, J ).NE.ZERO ) THEN 625 J2 = J + 1 626 JNEXT = J + 2 627 END IF 628 END IF 629* 630 IF( J1.EQ.J2 ) THEN 631* 632* 1 by 1 diagonal block 633* 634* Scale if necessary to avoid overflow in forming the 635* right-hand side element by inner product. 636* 637 XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) ) 638 IF( XMAX.GT.ONE ) THEN 639 REC = ONE / XMAX 640 IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN 641 CALL DSCAL( N2, REC, X, 1 ) 642 SCALE = SCALE*REC 643 XMAX = XMAX*REC 644 END IF 645 END IF 646* 647 X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 ) 648 X( N+J1 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1, 649 $ X( N+1 ), 1 ) 650 IF( J1.GT.1 ) THEN 651 X( J1 ) = X( J1 ) - B( J1 )*X( N+1 ) 652 X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 ) 653 END IF 654 XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) ) 655* 656 Z = W 657 IF( J1.EQ.1 ) 658 $ Z = B( 1 ) 659* 660* Scale if necessary to avoid overflow in 661* complex division 662* 663 TJJ = ABS( T( J1, J1 ) ) + ABS( Z ) 664 TMP = T( J1, J1 ) 665 IF( TJJ.LT.SMINW ) THEN 666 TMP = SMINW 667 TJJ = SMINW 668 INFO = 1 669 END IF 670* 671 IF( TJJ.LT.ONE ) THEN 672 IF( XJ.GT.BIGNUM*TJJ ) THEN 673 REC = ONE / XJ 674 CALL DSCAL( N2, REC, X, 1 ) 675 SCALE = SCALE*REC 676 XMAX = XMAX*REC 677 END IF 678 END IF 679 CALL DLADIV( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI ) 680 X( J1 ) = SR 681 X( J1+N ) = SI 682 XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX ) 683* 684 ELSE 685* 686* 2 by 2 diagonal block 687* 688* Scale if necessary to avoid overflow in forming the 689* right-hand side element by inner product. 690* 691 XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ), 692 $ ABS( X( J2 ) )+ABS( X( N+J2 ) ) ) 693 IF( XMAX.GT.ONE ) THEN 694 REC = ONE / XMAX 695 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT. 696 $ ( BIGNUM-XJ ) / XMAX ) THEN 697 CALL DSCAL( N2, REC, X, 1 ) 698 SCALE = SCALE*REC 699 XMAX = XMAX*REC 700 END IF 701 END IF 702* 703 D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 704 $ 1 ) 705 D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X, 706 $ 1 ) 707 D( 1, 2 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1, 708 $ X( N+1 ), 1 ) 709 D( 2, 2 ) = X( N+J2 ) - DDOT( J1-1, T( 1, J2 ), 1, 710 $ X( N+1 ), 1 ) 711 D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 ) 712 D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 ) 713 D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 ) 714 D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 ) 715* 716 CALL DLALN2( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ), 717 $ LDT, ONE, ONE, D, 2, ZERO, W, V, 2, 718 $ SCALOC, XNORM, IERR ) 719 IF( IERR.NE.0 ) 720 $ INFO = 2 721* 722 IF( SCALOC.NE.ONE ) THEN 723 CALL DSCAL( N2, SCALOC, X, 1 ) 724 SCALE = SCALOC*SCALE 725 END IF 726 X( J1 ) = V( 1, 1 ) 727 X( J2 ) = V( 2, 1 ) 728 X( N+J1 ) = V( 1, 2 ) 729 X( N+J2 ) = V( 2, 2 ) 730 XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ), 731 $ ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX ) 732* 733 END IF 734* 735 80 CONTINUE 736* 737 END IF 738* 739 END IF 740* 741 RETURN 742* 743* End of DLAQTR 744* 745 END 746