Annotation of rpl/lapack/lapack/dlaqtr.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaqtr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack