|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL LREAL, LTRAN 00011 INTEGER INFO, LDT, N 00012 DOUBLE PRECISION SCALE, W 00013 * .. 00014 * .. Array Arguments .. 00015 DOUBLE PRECISION B( * ), T( LDT, * ), WORK( * ), X( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLAQTR solves the real quasi-triangular system 00022 * 00023 * op(T)*p = scale*c, if LREAL = .TRUE. 00024 * 00025 * or the complex quasi-triangular systems 00026 * 00027 * op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE. 00028 * 00029 * in real arithmetic, where T is upper quasi-triangular. 00030 * If LREAL = .FALSE., then the first diagonal block of T must be 00031 * 1 by 1, B is the specially structured matrix 00032 * 00033 * B = [ b(1) b(2) ... b(n) ] 00034 * [ w ] 00035 * [ w ] 00036 * [ . ] 00037 * [ w ] 00038 * 00039 * op(A) = A or A', A' denotes the conjugate transpose of 00040 * matrix A. 00041 * 00042 * On input, X = [ c ]. On output, X = [ p ]. 00043 * [ d ] [ q ] 00044 * 00045 * This subroutine is designed for the condition number estimation 00046 * in routine DTRSNA. 00047 * 00048 * Arguments 00049 * ========= 00050 * 00051 * LTRAN (input) LOGICAL 00052 * On entry, LTRAN specifies the option of conjugate transpose: 00053 * = .FALSE., op(T+i*B) = T+i*B, 00054 * = .TRUE., op(T+i*B) = (T+i*B)'. 00055 * 00056 * LREAL (input) LOGICAL 00057 * On entry, LREAL specifies the input matrix structure: 00058 * = .FALSE., the input is complex 00059 * = .TRUE., the input is real 00060 * 00061 * N (input) INTEGER 00062 * On entry, N specifies the order of T+i*B. N >= 0. 00063 * 00064 * T (input) DOUBLE PRECISION array, dimension (LDT,N) 00065 * On entry, T contains a matrix in Schur canonical form. 00066 * If LREAL = .FALSE., then the first diagonal block of T mu 00067 * be 1 by 1. 00068 * 00069 * LDT (input) INTEGER 00070 * The leading dimension of the matrix T. LDT >= max(1,N). 00071 * 00072 * B (input) DOUBLE PRECISION array, dimension (N) 00073 * On entry, B contains the elements to form the matrix 00074 * B as described above. 00075 * If LREAL = .TRUE., B is not referenced. 00076 * 00077 * W (input) DOUBLE PRECISION 00078 * On entry, W is the diagonal element of the matrix B. 00079 * If LREAL = .TRUE., W is not referenced. 00080 * 00081 * SCALE (output) DOUBLE PRECISION 00082 * On exit, SCALE is the scale factor. 00083 * 00084 * X (input/output) DOUBLE PRECISION array, dimension (2*N) 00085 * On entry, X contains the right hand side of the system. 00086 * On exit, X is overwritten by the solution. 00087 * 00088 * WORK (workspace) DOUBLE PRECISION array, dimension (N) 00089 * 00090 * INFO (output) INTEGER 00091 * On exit, INFO is set to 00092 * 0: successful exit. 00093 * 1: the some diagonal 1 by 1 block has been perturbed by 00094 * a small number SMIN to keep nonsingularity. 00095 * 2: the some diagonal 2 by 2 block has been perturbed by 00096 * a small number in DLALN2 to keep nonsingularity. 00097 * NOTE: In the interests of speed, this routine does not 00098 * check the inputs for errors. 00099 * 00100 * ===================================================================== 00101 * 00102 * .. Parameters .. 00103 DOUBLE PRECISION ZERO, ONE 00104 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00105 * .. 00106 * .. Local Scalars .. 00107 LOGICAL NOTRAN 00108 INTEGER I, IERR, J, J1, J2, JNEXT, K, N1, N2 00109 DOUBLE PRECISION BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW, 00110 $ SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z 00111 * .. 00112 * .. Local Arrays .. 00113 DOUBLE PRECISION D( 2, 2 ), V( 2, 2 ) 00114 * .. 00115 * .. External Functions .. 00116 INTEGER IDAMAX 00117 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE 00118 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE 00119 * .. 00120 * .. External Subroutines .. 00121 EXTERNAL DAXPY, DLADIV, DLALN2, DSCAL 00122 * .. 00123 * .. Intrinsic Functions .. 00124 INTRINSIC ABS, MAX 00125 * .. 00126 * .. Executable Statements .. 00127 * 00128 * Do not test the input parameters for errors 00129 * 00130 NOTRAN = .NOT.LTRAN 00131 INFO = 0 00132 * 00133 * Quick return if possible 00134 * 00135 IF( N.EQ.0 ) 00136 $ RETURN 00137 * 00138 * Set constants to control overflow 00139 * 00140 EPS = DLAMCH( 'P' ) 00141 SMLNUM = DLAMCH( 'S' ) / EPS 00142 BIGNUM = ONE / SMLNUM 00143 * 00144 XNORM = DLANGE( 'M', N, N, T, LDT, D ) 00145 IF( .NOT.LREAL ) 00146 $ XNORM = MAX( XNORM, ABS( W ), DLANGE( 'M', N, 1, B, N, D ) ) 00147 SMIN = MAX( SMLNUM, EPS*XNORM ) 00148 * 00149 * Compute 1-norm of each column of strictly upper triangular 00150 * part of T to control overflow in triangular solver. 00151 * 00152 WORK( 1 ) = ZERO 00153 DO 10 J = 2, N 00154 WORK( J ) = DASUM( J-1, T( 1, J ), 1 ) 00155 10 CONTINUE 00156 * 00157 IF( .NOT.LREAL ) THEN 00158 DO 20 I = 2, N 00159 WORK( I ) = WORK( I ) + ABS( B( I ) ) 00160 20 CONTINUE 00161 END IF 00162 * 00163 N2 = 2*N 00164 N1 = N 00165 IF( .NOT.LREAL ) 00166 $ N1 = N2 00167 K = IDAMAX( N1, X, 1 ) 00168 XMAX = ABS( X( K ) ) 00169 SCALE = ONE 00170 * 00171 IF( XMAX.GT.BIGNUM ) THEN 00172 SCALE = BIGNUM / XMAX 00173 CALL DSCAL( N1, SCALE, X, 1 ) 00174 XMAX = BIGNUM 00175 END IF 00176 * 00177 IF( LREAL ) THEN 00178 * 00179 IF( NOTRAN ) THEN 00180 * 00181 * Solve T*p = scale*c 00182 * 00183 JNEXT = N 00184 DO 30 J = N, 1, -1 00185 IF( J.GT.JNEXT ) 00186 $ GO TO 30 00187 J1 = J 00188 J2 = J 00189 JNEXT = J - 1 00190 IF( J.GT.1 ) THEN 00191 IF( T( J, J-1 ).NE.ZERO ) THEN 00192 J1 = J - 1 00193 JNEXT = J - 2 00194 END IF 00195 END IF 00196 * 00197 IF( J1.EQ.J2 ) THEN 00198 * 00199 * Meet 1 by 1 diagonal block 00200 * 00201 * Scale to avoid overflow when computing 00202 * x(j) = b(j)/T(j,j) 00203 * 00204 XJ = ABS( X( J1 ) ) 00205 TJJ = ABS( T( J1, J1 ) ) 00206 TMP = T( J1, J1 ) 00207 IF( TJJ.LT.SMIN ) THEN 00208 TMP = SMIN 00209 TJJ = SMIN 00210 INFO = 1 00211 END IF 00212 * 00213 IF( XJ.EQ.ZERO ) 00214 $ GO TO 30 00215 * 00216 IF( TJJ.LT.ONE ) THEN 00217 IF( XJ.GT.BIGNUM*TJJ ) THEN 00218 REC = ONE / XJ 00219 CALL DSCAL( N, REC, X, 1 ) 00220 SCALE = SCALE*REC 00221 XMAX = XMAX*REC 00222 END IF 00223 END IF 00224 X( J1 ) = X( J1 ) / TMP 00225 XJ = ABS( X( J1 ) ) 00226 * 00227 * Scale x if necessary to avoid overflow when adding a 00228 * multiple of column j1 of T. 00229 * 00230 IF( XJ.GT.ONE ) THEN 00231 REC = ONE / XJ 00232 IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN 00233 CALL DSCAL( N, REC, X, 1 ) 00234 SCALE = SCALE*REC 00235 END IF 00236 END IF 00237 IF( J1.GT.1 ) THEN 00238 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 ) 00239 K = IDAMAX( J1-1, X, 1 ) 00240 XMAX = ABS( X( K ) ) 00241 END IF 00242 * 00243 ELSE 00244 * 00245 * Meet 2 by 2 diagonal block 00246 * 00247 * Call 2 by 2 linear system solve, to take 00248 * care of possible overflow by scaling factor. 00249 * 00250 D( 1, 1 ) = X( J1 ) 00251 D( 2, 1 ) = X( J2 ) 00252 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ), 00253 $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2, 00254 $ SCALOC, XNORM, IERR ) 00255 IF( IERR.NE.0 ) 00256 $ INFO = 2 00257 * 00258 IF( SCALOC.NE.ONE ) THEN 00259 CALL DSCAL( N, SCALOC, X, 1 ) 00260 SCALE = SCALE*SCALOC 00261 END IF 00262 X( J1 ) = V( 1, 1 ) 00263 X( J2 ) = V( 2, 1 ) 00264 * 00265 * Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2)) 00266 * to avoid overflow in updating right-hand side. 00267 * 00268 XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) ) 00269 IF( XJ.GT.ONE ) THEN 00270 REC = ONE / XJ 00271 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT. 00272 $ ( BIGNUM-XMAX )*REC ) THEN 00273 CALL DSCAL( N, REC, X, 1 ) 00274 SCALE = SCALE*REC 00275 END IF 00276 END IF 00277 * 00278 * Update right-hand side 00279 * 00280 IF( J1.GT.1 ) THEN 00281 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 ) 00282 CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 ) 00283 K = IDAMAX( J1-1, X, 1 ) 00284 XMAX = ABS( X( K ) ) 00285 END IF 00286 * 00287 END IF 00288 * 00289 30 CONTINUE 00290 * 00291 ELSE 00292 * 00293 * Solve T'*p = scale*c 00294 * 00295 JNEXT = 1 00296 DO 40 J = 1, N 00297 IF( J.LT.JNEXT ) 00298 $ GO TO 40 00299 J1 = J 00300 J2 = J 00301 JNEXT = J + 1 00302 IF( J.LT.N ) THEN 00303 IF( T( J+1, J ).NE.ZERO ) THEN 00304 J2 = J + 1 00305 JNEXT = J + 2 00306 END IF 00307 END IF 00308 * 00309 IF( J1.EQ.J2 ) THEN 00310 * 00311 * 1 by 1 diagonal block 00312 * 00313 * Scale if necessary to avoid overflow in forming the 00314 * right-hand side element by inner product. 00315 * 00316 XJ = ABS( X( J1 ) ) 00317 IF( XMAX.GT.ONE ) THEN 00318 REC = ONE / XMAX 00319 IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN 00320 CALL DSCAL( N, REC, X, 1 ) 00321 SCALE = SCALE*REC 00322 XMAX = XMAX*REC 00323 END IF 00324 END IF 00325 * 00326 X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 ) 00327 * 00328 XJ = ABS( X( J1 ) ) 00329 TJJ = ABS( T( J1, J1 ) ) 00330 TMP = T( J1, J1 ) 00331 IF( TJJ.LT.SMIN ) THEN 00332 TMP = SMIN 00333 TJJ = SMIN 00334 INFO = 1 00335 END IF 00336 * 00337 IF( TJJ.LT.ONE ) THEN 00338 IF( XJ.GT.BIGNUM*TJJ ) THEN 00339 REC = ONE / XJ 00340 CALL DSCAL( N, REC, X, 1 ) 00341 SCALE = SCALE*REC 00342 XMAX = XMAX*REC 00343 END IF 00344 END IF 00345 X( J1 ) = X( J1 ) / TMP 00346 XMAX = MAX( XMAX, ABS( X( J1 ) ) ) 00347 * 00348 ELSE 00349 * 00350 * 2 by 2 diagonal block 00351 * 00352 * Scale if necessary to avoid overflow in forming the 00353 * right-hand side elements by inner product. 00354 * 00355 XJ = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ) ) 00356 IF( XMAX.GT.ONE ) THEN 00357 REC = ONE / XMAX 00358 IF( MAX( WORK( J2 ), WORK( J1 ) ).GT.( BIGNUM-XJ )* 00359 $ REC ) THEN 00360 CALL DSCAL( N, REC, X, 1 ) 00361 SCALE = SCALE*REC 00362 XMAX = XMAX*REC 00363 END IF 00364 END IF 00365 * 00366 D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 00367 $ 1 ) 00368 D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X, 00369 $ 1 ) 00370 * 00371 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J1, J1 ), 00372 $ LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2, 00373 $ SCALOC, XNORM, IERR ) 00374 IF( IERR.NE.0 ) 00375 $ INFO = 2 00376 * 00377 IF( SCALOC.NE.ONE ) THEN 00378 CALL DSCAL( N, SCALOC, X, 1 ) 00379 SCALE = SCALE*SCALOC 00380 END IF 00381 X( J1 ) = V( 1, 1 ) 00382 X( J2 ) = V( 2, 1 ) 00383 XMAX = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ), XMAX ) 00384 * 00385 END IF 00386 40 CONTINUE 00387 END IF 00388 * 00389 ELSE 00390 * 00391 SMINW = MAX( EPS*ABS( W ), SMIN ) 00392 IF( NOTRAN ) THEN 00393 * 00394 * Solve (T + iB)*(p+iq) = c+id 00395 * 00396 JNEXT = N 00397 DO 70 J = N, 1, -1 00398 IF( J.GT.JNEXT ) 00399 $ GO TO 70 00400 J1 = J 00401 J2 = J 00402 JNEXT = J - 1 00403 IF( J.GT.1 ) THEN 00404 IF( T( J, J-1 ).NE.ZERO ) THEN 00405 J1 = J - 1 00406 JNEXT = J - 2 00407 END IF 00408 END IF 00409 * 00410 IF( J1.EQ.J2 ) THEN 00411 * 00412 * 1 by 1 diagonal block 00413 * 00414 * Scale if necessary to avoid overflow in division 00415 * 00416 Z = W 00417 IF( J1.EQ.1 ) 00418 $ Z = B( 1 ) 00419 XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) ) 00420 TJJ = ABS( T( J1, J1 ) ) + ABS( Z ) 00421 TMP = T( J1, J1 ) 00422 IF( TJJ.LT.SMINW ) THEN 00423 TMP = SMINW 00424 TJJ = SMINW 00425 INFO = 1 00426 END IF 00427 * 00428 IF( XJ.EQ.ZERO ) 00429 $ GO TO 70 00430 * 00431 IF( TJJ.LT.ONE ) THEN 00432 IF( XJ.GT.BIGNUM*TJJ ) THEN 00433 REC = ONE / XJ 00434 CALL DSCAL( N2, REC, X, 1 ) 00435 SCALE = SCALE*REC 00436 XMAX = XMAX*REC 00437 END IF 00438 END IF 00439 CALL DLADIV( X( J1 ), X( N+J1 ), TMP, Z, SR, SI ) 00440 X( J1 ) = SR 00441 X( N+J1 ) = SI 00442 XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) ) 00443 * 00444 * Scale x if necessary to avoid overflow when adding a 00445 * multiple of column j1 of T. 00446 * 00447 IF( XJ.GT.ONE ) THEN 00448 REC = ONE / XJ 00449 IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN 00450 CALL DSCAL( N2, REC, X, 1 ) 00451 SCALE = SCALE*REC 00452 END IF 00453 END IF 00454 * 00455 IF( J1.GT.1 ) THEN 00456 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 ) 00457 CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1, 00458 $ X( N+1 ), 1 ) 00459 * 00460 X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) 00461 X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) 00462 * 00463 XMAX = ZERO 00464 DO 50 K = 1, J1 - 1 00465 XMAX = MAX( XMAX, ABS( X( K ) )+ 00466 $ ABS( X( K+N ) ) ) 00467 50 CONTINUE 00468 END IF 00469 * 00470 ELSE 00471 * 00472 * Meet 2 by 2 diagonal block 00473 * 00474 D( 1, 1 ) = X( J1 ) 00475 D( 2, 1 ) = X( J2 ) 00476 D( 1, 2 ) = X( N+J1 ) 00477 D( 2, 2 ) = X( N+J2 ) 00478 CALL DLALN2( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ), 00479 $ LDT, ONE, ONE, D, 2, ZERO, -W, V, 2, 00480 $ SCALOC, XNORM, IERR ) 00481 IF( IERR.NE.0 ) 00482 $ INFO = 2 00483 * 00484 IF( SCALOC.NE.ONE ) THEN 00485 CALL DSCAL( 2*N, SCALOC, X, 1 ) 00486 SCALE = SCALOC*SCALE 00487 END IF 00488 X( J1 ) = V( 1, 1 ) 00489 X( J2 ) = V( 2, 1 ) 00490 X( N+J1 ) = V( 1, 2 ) 00491 X( N+J2 ) = V( 2, 2 ) 00492 * 00493 * Scale X(J1), .... to avoid overflow in 00494 * updating right hand side. 00495 * 00496 XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ), 00497 $ ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) ) 00498 IF( XJ.GT.ONE ) THEN 00499 REC = ONE / XJ 00500 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT. 00501 $ ( BIGNUM-XMAX )*REC ) THEN 00502 CALL DSCAL( N2, REC, X, 1 ) 00503 SCALE = SCALE*REC 00504 END IF 00505 END IF 00506 * 00507 * Update the right-hand side. 00508 * 00509 IF( J1.GT.1 ) THEN 00510 CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 ) 00511 CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 ) 00512 * 00513 CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1, 00514 $ X( N+1 ), 1 ) 00515 CALL DAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1, 00516 $ X( N+1 ), 1 ) 00517 * 00518 X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) + 00519 $ B( J2 )*X( N+J2 ) 00520 X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) - 00521 $ B( J2 )*X( J2 ) 00522 * 00523 XMAX = ZERO 00524 DO 60 K = 1, J1 - 1 00525 XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ), 00526 $ XMAX ) 00527 60 CONTINUE 00528 END IF 00529 * 00530 END IF 00531 70 CONTINUE 00532 * 00533 ELSE 00534 * 00535 * Solve (T + iB)'*(p+iq) = c+id 00536 * 00537 JNEXT = 1 00538 DO 80 J = 1, N 00539 IF( J.LT.JNEXT ) 00540 $ GO TO 80 00541 J1 = J 00542 J2 = J 00543 JNEXT = J + 1 00544 IF( J.LT.N ) THEN 00545 IF( T( J+1, J ).NE.ZERO ) THEN 00546 J2 = J + 1 00547 JNEXT = J + 2 00548 END IF 00549 END IF 00550 * 00551 IF( J1.EQ.J2 ) THEN 00552 * 00553 * 1 by 1 diagonal block 00554 * 00555 * Scale if necessary to avoid overflow in forming the 00556 * right-hand side element by inner product. 00557 * 00558 XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) ) 00559 IF( XMAX.GT.ONE ) THEN 00560 REC = ONE / XMAX 00561 IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN 00562 CALL DSCAL( N2, REC, X, 1 ) 00563 SCALE = SCALE*REC 00564 XMAX = XMAX*REC 00565 END IF 00566 END IF 00567 * 00568 X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 ) 00569 X( N+J1 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1, 00570 $ X( N+1 ), 1 ) 00571 IF( J1.GT.1 ) THEN 00572 X( J1 ) = X( J1 ) - B( J1 )*X( N+1 ) 00573 X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 ) 00574 END IF 00575 XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) ) 00576 * 00577 Z = W 00578 IF( J1.EQ.1 ) 00579 $ Z = B( 1 ) 00580 * 00581 * Scale if necessary to avoid overflow in 00582 * complex division 00583 * 00584 TJJ = ABS( T( J1, J1 ) ) + ABS( Z ) 00585 TMP = T( J1, J1 ) 00586 IF( TJJ.LT.SMINW ) THEN 00587 TMP = SMINW 00588 TJJ = SMINW 00589 INFO = 1 00590 END IF 00591 * 00592 IF( TJJ.LT.ONE ) THEN 00593 IF( XJ.GT.BIGNUM*TJJ ) THEN 00594 REC = ONE / XJ 00595 CALL DSCAL( N2, REC, X, 1 ) 00596 SCALE = SCALE*REC 00597 XMAX = XMAX*REC 00598 END IF 00599 END IF 00600 CALL DLADIV( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI ) 00601 X( J1 ) = SR 00602 X( J1+N ) = SI 00603 XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX ) 00604 * 00605 ELSE 00606 * 00607 * 2 by 2 diagonal block 00608 * 00609 * Scale if necessary to avoid overflow in forming the 00610 * right-hand side element by inner product. 00611 * 00612 XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ), 00613 $ ABS( X( J2 ) )+ABS( X( N+J2 ) ) ) 00614 IF( XMAX.GT.ONE ) THEN 00615 REC = ONE / XMAX 00616 IF( MAX( WORK( J1 ), WORK( J2 ) ).GT. 00617 $ ( BIGNUM-XJ ) / XMAX ) THEN 00618 CALL DSCAL( N2, REC, X, 1 ) 00619 SCALE = SCALE*REC 00620 XMAX = XMAX*REC 00621 END IF 00622 END IF 00623 * 00624 D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 00625 $ 1 ) 00626 D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X, 00627 $ 1 ) 00628 D( 1, 2 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1, 00629 $ X( N+1 ), 1 ) 00630 D( 2, 2 ) = X( N+J2 ) - DDOT( J1-1, T( 1, J2 ), 1, 00631 $ X( N+1 ), 1 ) 00632 D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 ) 00633 D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 ) 00634 D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 ) 00635 D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 ) 00636 * 00637 CALL DLALN2( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ), 00638 $ LDT, ONE, ONE, D, 2, ZERO, W, V, 2, 00639 $ SCALOC, XNORM, IERR ) 00640 IF( IERR.NE.0 ) 00641 $ INFO = 2 00642 * 00643 IF( SCALOC.NE.ONE ) THEN 00644 CALL DSCAL( N2, SCALOC, X, 1 ) 00645 SCALE = SCALOC*SCALE 00646 END IF 00647 X( J1 ) = V( 1, 1 ) 00648 X( J2 ) = V( 2, 1 ) 00649 X( N+J1 ) = V( 1, 2 ) 00650 X( N+J2 ) = V( 2, 2 ) 00651 XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ), 00652 $ ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX ) 00653 * 00654 END IF 00655 * 00656 80 CONTINUE 00657 * 00658 END IF 00659 * 00660 END IF 00661 * 00662 RETURN 00663 * 00664 * End of DLAQTR 00665 * 00666 END
1.7.3 
