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