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