![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLALN2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLALN2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaln2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaln2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaln2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, 00022 * LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * LOGICAL LTRANS 00026 * INTEGER INFO, LDA, LDB, LDX, NA, NW 00027 * DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM 00028 * .. 00029 * .. Array Arguments .. 00030 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> DLALN2 solves a system of the form (ca A - w D ) X = s B 00040 *> or (ca A**T - w D) X = s B with possible scaling ("s") and 00041 *> perturbation of A. (A**T means A-transpose.) 00042 *> 00043 *> A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA 00044 *> real diagonal matrix, w is a real or complex value, and X and B are 00045 *> NA x 1 matrices -- real if w is real, complex if w is complex. NA 00046 *> may be 1 or 2. 00047 *> 00048 *> If w is complex, X and B are represented as NA x 2 matrices, 00049 *> the first column of each being the real part and the second 00050 *> being the imaginary part. 00051 *> 00052 *> "s" is a scaling factor (.LE. 1), computed by DLALN2, which is 00053 *> so chosen that X can be computed without overflow. X is further 00054 *> scaled if necessary to assure that norm(ca A - w D)*norm(X) is less 00055 *> than overflow. 00056 *> 00057 *> If both singular values of (ca A - w D) are less than SMIN, 00058 *> SMIN*identity will be used instead of (ca A - w D). If only one 00059 *> singular value is less than SMIN, one element of (ca A - w D) will be 00060 *> perturbed enough to make the smallest singular value roughly SMIN. 00061 *> If both singular values are at least SMIN, (ca A - w D) will not be 00062 *> perturbed. In any case, the perturbation will be at most some small 00063 *> multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values 00064 *> are computed by infinity-norm approximations, and thus will only be 00065 *> correct to a factor of 2 or so. 00066 *> 00067 *> Note: all input quantities are assumed to be smaller than overflow 00068 *> by a reasonable factor. (See BIGNUM.) 00069 *> \endverbatim 00070 * 00071 * Arguments: 00072 * ========== 00073 * 00074 *> \param[in] LTRANS 00075 *> \verbatim 00076 *> LTRANS is LOGICAL 00077 *> =.TRUE.: A-transpose will be used. 00078 *> =.FALSE.: A will be used (not transposed.) 00079 *> \endverbatim 00080 *> 00081 *> \param[in] NA 00082 *> \verbatim 00083 *> NA is INTEGER 00084 *> The size of the matrix A. It may (only) be 1 or 2. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] NW 00088 *> \verbatim 00089 *> NW is INTEGER 00090 *> 1 if "w" is real, 2 if "w" is complex. It may only be 1 00091 *> or 2. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] SMIN 00095 *> \verbatim 00096 *> SMIN is DOUBLE PRECISION 00097 *> The desired lower bound on the singular values of A. This 00098 *> should be a safe distance away from underflow or overflow, 00099 *> say, between (underflow/machine precision) and (machine 00100 *> precision * overflow ). (See BIGNUM and ULP.) 00101 *> \endverbatim 00102 *> 00103 *> \param[in] CA 00104 *> \verbatim 00105 *> CA is DOUBLE PRECISION 00106 *> The coefficient c, which A is multiplied by. 00107 *> \endverbatim 00108 *> 00109 *> \param[in] A 00110 *> \verbatim 00111 *> A is DOUBLE PRECISION array, dimension (LDA,NA) 00112 *> The NA x NA matrix A. 00113 *> \endverbatim 00114 *> 00115 *> \param[in] LDA 00116 *> \verbatim 00117 *> LDA is INTEGER 00118 *> The leading dimension of A. It must be at least NA. 00119 *> \endverbatim 00120 *> 00121 *> \param[in] D1 00122 *> \verbatim 00123 *> D1 is DOUBLE PRECISION 00124 *> The 1,1 element in the diagonal matrix D. 00125 *> \endverbatim 00126 *> 00127 *> \param[in] D2 00128 *> \verbatim 00129 *> D2 is DOUBLE PRECISION 00130 *> The 2,2 element in the diagonal matrix D. Not used if NW=1. 00131 *> \endverbatim 00132 *> 00133 *> \param[in] B 00134 *> \verbatim 00135 *> B is DOUBLE PRECISION array, dimension (LDB,NW) 00136 *> The NA x NW matrix B (right-hand side). If NW=2 ("w" is 00137 *> complex), column 1 contains the real part of B and column 2 00138 *> contains the imaginary part. 00139 *> \endverbatim 00140 *> 00141 *> \param[in] LDB 00142 *> \verbatim 00143 *> LDB is INTEGER 00144 *> The leading dimension of B. It must be at least NA. 00145 *> \endverbatim 00146 *> 00147 *> \param[in] WR 00148 *> \verbatim 00149 *> WR is DOUBLE PRECISION 00150 *> The real part of the scalar "w". 00151 *> \endverbatim 00152 *> 00153 *> \param[in] WI 00154 *> \verbatim 00155 *> WI is DOUBLE PRECISION 00156 *> The imaginary part of the scalar "w". Not used if NW=1. 00157 *> \endverbatim 00158 *> 00159 *> \param[out] X 00160 *> \verbatim 00161 *> X is DOUBLE PRECISION array, dimension (LDX,NW) 00162 *> The NA x NW matrix X (unknowns), as computed by DLALN2. 00163 *> If NW=2 ("w" is complex), on exit, column 1 will contain 00164 *> the real part of X and column 2 will contain the imaginary 00165 *> part. 00166 *> \endverbatim 00167 *> 00168 *> \param[in] LDX 00169 *> \verbatim 00170 *> LDX is INTEGER 00171 *> The leading dimension of X. It must be at least NA. 00172 *> \endverbatim 00173 *> 00174 *> \param[out] SCALE 00175 *> \verbatim 00176 *> SCALE is DOUBLE PRECISION 00177 *> The scale factor that B must be multiplied by to insure 00178 *> that overflow does not occur when computing X. Thus, 00179 *> (ca A - w D) X will be SCALE*B, not B (ignoring 00180 *> perturbations of A.) It will be at most 1. 00181 *> \endverbatim 00182 *> 00183 *> \param[out] XNORM 00184 *> \verbatim 00185 *> XNORM is DOUBLE PRECISION 00186 *> The infinity-norm of X, when X is regarded as an NA x NW 00187 *> real matrix. 00188 *> \endverbatim 00189 *> 00190 *> \param[out] INFO 00191 *> \verbatim 00192 *> INFO is INTEGER 00193 *> An error flag. It will be set to zero if no error occurs, 00194 *> a negative number if an argument is in error, or a positive 00195 *> number if ca A - w D had to be perturbed. 00196 *> The possible values are: 00197 *> = 0: No error occurred, and (ca A - w D) did not have to be 00198 *> perturbed. 00199 *> = 1: (ca A - w D) had to be perturbed to make its smallest 00200 *> (or only) singular value greater than SMIN. 00201 *> NOTE: In the interests of speed, this routine does not 00202 *> check the inputs for errors. 00203 *> \endverbatim 00204 * 00205 * Authors: 00206 * ======== 00207 * 00208 *> \author Univ. of Tennessee 00209 *> \author Univ. of California Berkeley 00210 *> \author Univ. of Colorado Denver 00211 *> \author NAG Ltd. 00212 * 00213 *> \date November 2011 00214 * 00215 *> \ingroup doubleOTHERauxiliary 00216 * 00217 * ===================================================================== 00218 SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, 00219 $ LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) 00220 * 00221 * -- LAPACK auxiliary routine (version 3.4.0) -- 00222 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00223 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00224 * November 2011 00225 * 00226 * .. Scalar Arguments .. 00227 LOGICAL LTRANS 00228 INTEGER INFO, LDA, LDB, LDX, NA, NW 00229 DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM 00230 * .. 00231 * .. Array Arguments .. 00232 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ) 00233 * .. 00234 * 00235 * ===================================================================== 00236 * 00237 * .. Parameters .. 00238 DOUBLE PRECISION ZERO, ONE 00239 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00240 DOUBLE PRECISION TWO 00241 PARAMETER ( TWO = 2.0D0 ) 00242 * .. 00243 * .. Local Scalars .. 00244 INTEGER ICMAX, J 00245 DOUBLE PRECISION BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21, 00246 $ CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21, 00247 $ LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R, 00248 $ UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S, 00249 $ UR22, XI1, XI2, XR1, XR2 00250 * .. 00251 * .. Local Arrays .. 00252 LOGICAL RSWAP( 4 ), ZSWAP( 4 ) 00253 INTEGER IPIVOT( 4, 4 ) 00254 DOUBLE PRECISION CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 ) 00255 * .. 00256 * .. External Functions .. 00257 DOUBLE PRECISION DLAMCH 00258 EXTERNAL DLAMCH 00259 * .. 00260 * .. External Subroutines .. 00261 EXTERNAL DLADIV 00262 * .. 00263 * .. Intrinsic Functions .. 00264 INTRINSIC ABS, MAX 00265 * .. 00266 * .. Equivalences .. 00267 EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ), 00268 $ ( CR( 1, 1 ), CRV( 1 ) ) 00269 * .. 00270 * .. Data statements .. 00271 DATA ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. / 00272 DATA RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. / 00273 DATA IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4, 00274 $ 3, 2, 1 / 00275 * .. 00276 * .. Executable Statements .. 00277 * 00278 * Compute BIGNUM 00279 * 00280 SMLNUM = TWO*DLAMCH( 'Safe minimum' ) 00281 BIGNUM = ONE / SMLNUM 00282 SMINI = MAX( SMIN, SMLNUM ) 00283 * 00284 * Don't check for input errors 00285 * 00286 INFO = 0 00287 * 00288 * Standard Initializations 00289 * 00290 SCALE = ONE 00291 * 00292 IF( NA.EQ.1 ) THEN 00293 * 00294 * 1 x 1 (i.e., scalar) system C X = B 00295 * 00296 IF( NW.EQ.1 ) THEN 00297 * 00298 * Real 1x1 system. 00299 * 00300 * C = ca A - w D 00301 * 00302 CSR = CA*A( 1, 1 ) - WR*D1 00303 CNORM = ABS( CSR ) 00304 * 00305 * If | C | < SMINI, use C = SMINI 00306 * 00307 IF( CNORM.LT.SMINI ) THEN 00308 CSR = SMINI 00309 CNORM = SMINI 00310 INFO = 1 00311 END IF 00312 * 00313 * Check scaling for X = B / C 00314 * 00315 BNORM = ABS( B( 1, 1 ) ) 00316 IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN 00317 IF( BNORM.GT.BIGNUM*CNORM ) 00318 $ SCALE = ONE / BNORM 00319 END IF 00320 * 00321 * Compute X 00322 * 00323 X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR 00324 XNORM = ABS( X( 1, 1 ) ) 00325 ELSE 00326 * 00327 * Complex 1x1 system (w is complex) 00328 * 00329 * C = ca A - w D 00330 * 00331 CSR = CA*A( 1, 1 ) - WR*D1 00332 CSI = -WI*D1 00333 CNORM = ABS( CSR ) + ABS( CSI ) 00334 * 00335 * If | C | < SMINI, use C = SMINI 00336 * 00337 IF( CNORM.LT.SMINI ) THEN 00338 CSR = SMINI 00339 CSI = ZERO 00340 CNORM = SMINI 00341 INFO = 1 00342 END IF 00343 * 00344 * Check scaling for X = B / C 00345 * 00346 BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) ) 00347 IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN 00348 IF( BNORM.GT.BIGNUM*CNORM ) 00349 $ SCALE = ONE / BNORM 00350 END IF 00351 * 00352 * Compute X 00353 * 00354 CALL DLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI, 00355 $ X( 1, 1 ), X( 1, 2 ) ) 00356 XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) ) 00357 END IF 00358 * 00359 ELSE 00360 * 00361 * 2x2 System 00362 * 00363 * Compute the real part of C = ca A - w D (or ca A**T - w D ) 00364 * 00365 CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1 00366 CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2 00367 IF( LTRANS ) THEN 00368 CR( 1, 2 ) = CA*A( 2, 1 ) 00369 CR( 2, 1 ) = CA*A( 1, 2 ) 00370 ELSE 00371 CR( 2, 1 ) = CA*A( 2, 1 ) 00372 CR( 1, 2 ) = CA*A( 1, 2 ) 00373 END IF 00374 * 00375 IF( NW.EQ.1 ) THEN 00376 * 00377 * Real 2x2 system (w is real) 00378 * 00379 * Find the largest element in C 00380 * 00381 CMAX = ZERO 00382 ICMAX = 0 00383 * 00384 DO 10 J = 1, 4 00385 IF( ABS( CRV( J ) ).GT.CMAX ) THEN 00386 CMAX = ABS( CRV( J ) ) 00387 ICMAX = J 00388 END IF 00389 10 CONTINUE 00390 * 00391 * If norm(C) < SMINI, use SMINI*identity. 00392 * 00393 IF( CMAX.LT.SMINI ) THEN 00394 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) ) 00395 IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN 00396 IF( BNORM.GT.BIGNUM*SMINI ) 00397 $ SCALE = ONE / BNORM 00398 END IF 00399 TEMP = SCALE / SMINI 00400 X( 1, 1 ) = TEMP*B( 1, 1 ) 00401 X( 2, 1 ) = TEMP*B( 2, 1 ) 00402 XNORM = TEMP*BNORM 00403 INFO = 1 00404 RETURN 00405 END IF 00406 * 00407 * Gaussian elimination with complete pivoting. 00408 * 00409 UR11 = CRV( ICMAX ) 00410 CR21 = CRV( IPIVOT( 2, ICMAX ) ) 00411 UR12 = CRV( IPIVOT( 3, ICMAX ) ) 00412 CR22 = CRV( IPIVOT( 4, ICMAX ) ) 00413 UR11R = ONE / UR11 00414 LR21 = UR11R*CR21 00415 UR22 = CR22 - UR12*LR21 00416 * 00417 * If smaller pivot < SMINI, use SMINI 00418 * 00419 IF( ABS( UR22 ).LT.SMINI ) THEN 00420 UR22 = SMINI 00421 INFO = 1 00422 END IF 00423 IF( RSWAP( ICMAX ) ) THEN 00424 BR1 = B( 2, 1 ) 00425 BR2 = B( 1, 1 ) 00426 ELSE 00427 BR1 = B( 1, 1 ) 00428 BR2 = B( 2, 1 ) 00429 END IF 00430 BR2 = BR2 - LR21*BR1 00431 BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) ) 00432 IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN 00433 IF( BBND.GE.BIGNUM*ABS( UR22 ) ) 00434 $ SCALE = ONE / BBND 00435 END IF 00436 * 00437 XR2 = ( BR2*SCALE ) / UR22 00438 XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 ) 00439 IF( ZSWAP( ICMAX ) ) THEN 00440 X( 1, 1 ) = XR2 00441 X( 2, 1 ) = XR1 00442 ELSE 00443 X( 1, 1 ) = XR1 00444 X( 2, 1 ) = XR2 00445 END IF 00446 XNORM = MAX( ABS( XR1 ), ABS( XR2 ) ) 00447 * 00448 * Further scaling if norm(A) norm(X) > overflow 00449 * 00450 IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN 00451 IF( XNORM.GT.BIGNUM / CMAX ) THEN 00452 TEMP = CMAX / BIGNUM 00453 X( 1, 1 ) = TEMP*X( 1, 1 ) 00454 X( 2, 1 ) = TEMP*X( 2, 1 ) 00455 XNORM = TEMP*XNORM 00456 SCALE = TEMP*SCALE 00457 END IF 00458 END IF 00459 ELSE 00460 * 00461 * Complex 2x2 system (w is complex) 00462 * 00463 * Find the largest element in C 00464 * 00465 CI( 1, 1 ) = -WI*D1 00466 CI( 2, 1 ) = ZERO 00467 CI( 1, 2 ) = ZERO 00468 CI( 2, 2 ) = -WI*D2 00469 CMAX = ZERO 00470 ICMAX = 0 00471 * 00472 DO 20 J = 1, 4 00473 IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN 00474 CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) ) 00475 ICMAX = J 00476 END IF 00477 20 CONTINUE 00478 * 00479 * If norm(C) < SMINI, use SMINI*identity. 00480 * 00481 IF( CMAX.LT.SMINI ) THEN 00482 BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ), 00483 $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) ) 00484 IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN 00485 IF( BNORM.GT.BIGNUM*SMINI ) 00486 $ SCALE = ONE / BNORM 00487 END IF 00488 TEMP = SCALE / SMINI 00489 X( 1, 1 ) = TEMP*B( 1, 1 ) 00490 X( 2, 1 ) = TEMP*B( 2, 1 ) 00491 X( 1, 2 ) = TEMP*B( 1, 2 ) 00492 X( 2, 2 ) = TEMP*B( 2, 2 ) 00493 XNORM = TEMP*BNORM 00494 INFO = 1 00495 RETURN 00496 END IF 00497 * 00498 * Gaussian elimination with complete pivoting. 00499 * 00500 UR11 = CRV( ICMAX ) 00501 UI11 = CIV( ICMAX ) 00502 CR21 = CRV( IPIVOT( 2, ICMAX ) ) 00503 CI21 = CIV( IPIVOT( 2, ICMAX ) ) 00504 UR12 = CRV( IPIVOT( 3, ICMAX ) ) 00505 UI12 = CIV( IPIVOT( 3, ICMAX ) ) 00506 CR22 = CRV( IPIVOT( 4, ICMAX ) ) 00507 CI22 = CIV( IPIVOT( 4, ICMAX ) ) 00508 IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN 00509 * 00510 * Code when off-diagonals of pivoted C are real 00511 * 00512 IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN 00513 TEMP = UI11 / UR11 00514 UR11R = ONE / ( UR11*( ONE+TEMP**2 ) ) 00515 UI11R = -TEMP*UR11R 00516 ELSE 00517 TEMP = UR11 / UI11 00518 UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) ) 00519 UR11R = -TEMP*UI11R 00520 END IF 00521 LR21 = CR21*UR11R 00522 LI21 = CR21*UI11R 00523 UR12S = UR12*UR11R 00524 UI12S = UR12*UI11R 00525 UR22 = CR22 - UR12*LR21 00526 UI22 = CI22 - UR12*LI21 00527 ELSE 00528 * 00529 * Code when diagonals of pivoted C are real 00530 * 00531 UR11R = ONE / UR11 00532 UI11R = ZERO 00533 LR21 = CR21*UR11R 00534 LI21 = CI21*UR11R 00535 UR12S = UR12*UR11R 00536 UI12S = UI12*UR11R 00537 UR22 = CR22 - UR12*LR21 + UI12*LI21 00538 UI22 = -UR12*LI21 - UI12*LR21 00539 END IF 00540 U22ABS = ABS( UR22 ) + ABS( UI22 ) 00541 * 00542 * If smaller pivot < SMINI, use SMINI 00543 * 00544 IF( U22ABS.LT.SMINI ) THEN 00545 UR22 = SMINI 00546 UI22 = ZERO 00547 INFO = 1 00548 END IF 00549 IF( RSWAP( ICMAX ) ) THEN 00550 BR2 = B( 1, 1 ) 00551 BR1 = B( 2, 1 ) 00552 BI2 = B( 1, 2 ) 00553 BI1 = B( 2, 2 ) 00554 ELSE 00555 BR1 = B( 1, 1 ) 00556 BR2 = B( 2, 1 ) 00557 BI1 = B( 1, 2 ) 00558 BI2 = B( 2, 2 ) 00559 END IF 00560 BR2 = BR2 - LR21*BR1 + LI21*BI1 00561 BI2 = BI2 - LI21*BR1 - LR21*BI1 00562 BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )* 00563 $ ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ), 00564 $ ABS( BR2 )+ABS( BI2 ) ) 00565 IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN 00566 IF( BBND.GE.BIGNUM*U22ABS ) THEN 00567 SCALE = ONE / BBND 00568 BR1 = SCALE*BR1 00569 BI1 = SCALE*BI1 00570 BR2 = SCALE*BR2 00571 BI2 = SCALE*BI2 00572 END IF 00573 END IF 00574 * 00575 CALL DLADIV( BR2, BI2, UR22, UI22, XR2, XI2 ) 00576 XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2 00577 XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2 00578 IF( ZSWAP( ICMAX ) ) THEN 00579 X( 1, 1 ) = XR2 00580 X( 2, 1 ) = XR1 00581 X( 1, 2 ) = XI2 00582 X( 2, 2 ) = XI1 00583 ELSE 00584 X( 1, 1 ) = XR1 00585 X( 2, 1 ) = XR2 00586 X( 1, 2 ) = XI1 00587 X( 2, 2 ) = XI2 00588 END IF 00589 XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) ) 00590 * 00591 * Further scaling if norm(A) norm(X) > overflow 00592 * 00593 IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN 00594 IF( XNORM.GT.BIGNUM / CMAX ) THEN 00595 TEMP = CMAX / BIGNUM 00596 X( 1, 1 ) = TEMP*X( 1, 1 ) 00597 X( 2, 1 ) = TEMP*X( 2, 1 ) 00598 X( 1, 2 ) = TEMP*X( 1, 2 ) 00599 X( 2, 2 ) = TEMP*X( 2, 2 ) 00600 XNORM = TEMP*XNORM 00601 SCALE = TEMP*SCALE 00602 END IF 00603 END IF 00604 END IF 00605 END IF 00606 * 00607 RETURN 00608 * 00609 * End of DLALN2 00610 * 00611 END