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