![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, 2: $ LDB, WR, WI, X, LDX, SCALE, XNORM, 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 LTRANS 11: INTEGER INFO, LDA, LDB, LDX, NA, NW 12: DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM 13: * .. 14: * .. Array Arguments .. 15: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * DLALN2 solves a system of the form (ca A - w D ) X = s B 22: * or (ca A' - w D) X = s B with possible scaling ("s") and 23: * perturbation of A. (A' means A-transpose.) 24: * 25: * A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA 26: * real diagonal matrix, w is a real or complex value, and X and B are 27: * NA x 1 matrices -- real if w is real, complex if w is complex. NA 28: * may be 1 or 2. 29: * 30: * If w is complex, X and B are represented as NA x 2 matrices, 31: * the first column of each being the real part and the second 32: * being the imaginary part. 33: * 34: * "s" is a scaling factor (.LE. 1), computed by DLALN2, which is 35: * so chosen that X can be computed without overflow. X is further 36: * scaled if necessary to assure that norm(ca A - w D)*norm(X) is less 37: * than overflow. 38: * 39: * If both singular values of (ca A - w D) are less than SMIN, 40: * SMIN*identity will be used instead of (ca A - w D). If only one 41: * singular value is less than SMIN, one element of (ca A - w D) will be 42: * perturbed enough to make the smallest singular value roughly SMIN. 43: * If both singular values are at least SMIN, (ca A - w D) will not be 44: * perturbed. In any case, the perturbation will be at most some small 45: * multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values 46: * are computed by infinity-norm approximations, and thus will only be 47: * correct to a factor of 2 or so. 48: * 49: * Note: all input quantities are assumed to be smaller than overflow 50: * by a reasonable factor. (See BIGNUM.) 51: * 52: * Arguments 53: * ========== 54: * 55: * LTRANS (input) LOGICAL 56: * =.TRUE.: A-transpose will be used. 57: * =.FALSE.: A will be used (not transposed.) 58: * 59: * NA (input) INTEGER 60: * The size of the matrix A. It may (only) be 1 or 2. 61: * 62: * NW (input) INTEGER 63: * 1 if "w" is real, 2 if "w" is complex. It may only be 1 64: * or 2. 65: * 66: * SMIN (input) DOUBLE PRECISION 67: * The desired lower bound on the singular values of A. This 68: * should be a safe distance away from underflow or overflow, 69: * say, between (underflow/machine precision) and (machine 70: * precision * overflow ). (See BIGNUM and ULP.) 71: * 72: * CA (input) DOUBLE PRECISION 73: * The coefficient c, which A is multiplied by. 74: * 75: * A (input) DOUBLE PRECISION array, dimension (LDA,NA) 76: * The NA x NA matrix A. 77: * 78: * LDA (input) INTEGER 79: * The leading dimension of A. It must be at least NA. 80: * 81: * D1 (input) DOUBLE PRECISION 82: * The 1,1 element in the diagonal matrix D. 83: * 84: * D2 (input) DOUBLE PRECISION 85: * The 2,2 element in the diagonal matrix D. Not used if NW=1. 86: * 87: * B (input) DOUBLE PRECISION array, dimension (LDB,NW) 88: * The NA x NW matrix B (right-hand side). If NW=2 ("w" is 89: * complex), column 1 contains the real part of B and column 2 90: * contains the imaginary part. 91: * 92: * LDB (input) INTEGER 93: * The leading dimension of B. It must be at least NA. 94: * 95: * WR (input) DOUBLE PRECISION 96: * The real part of the scalar "w". 97: * 98: * WI (input) DOUBLE PRECISION 99: * The imaginary part of the scalar "w". Not used if NW=1. 100: * 101: * X (output) DOUBLE PRECISION array, dimension (LDX,NW) 102: * The NA x NW matrix X (unknowns), as computed by DLALN2. 103: * If NW=2 ("w" is complex), on exit, column 1 will contain 104: * the real part of X and column 2 will contain the imaginary 105: * part. 106: * 107: * LDX (input) INTEGER 108: * The leading dimension of X. It must be at least NA. 109: * 110: * SCALE (output) DOUBLE PRECISION 111: * The scale factor that B must be multiplied by to insure 112: * that overflow does not occur when computing X. Thus, 113: * (ca A - w D) X will be SCALE*B, not B (ignoring 114: * perturbations of A.) It will be at most 1. 115: * 116: * XNORM (output) DOUBLE PRECISION 117: * The infinity-norm of X, when X is regarded as an NA x NW 118: * real matrix. 119: * 120: * INFO (output) INTEGER 121: * An error flag. It will be set to zero if no error occurs, 122: * a negative number if an argument is in error, or a positive 123: * number if ca A - w D had to be perturbed. 124: * The possible values are: 125: * = 0: No error occurred, and (ca A - w D) did not have to be 126: * perturbed. 127: * = 1: (ca A - w D) had to be perturbed to make its smallest 128: * (or only) singular value greater than SMIN. 129: * NOTE: In the interests of speed, this routine does not 130: * check the inputs for errors. 131: * 132: * ===================================================================== 133: * 134: * .. Parameters .. 135: DOUBLE PRECISION ZERO, ONE 136: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 137: DOUBLE PRECISION TWO 138: PARAMETER ( TWO = 2.0D0 ) 139: * .. 140: * .. Local Scalars .. 141: INTEGER ICMAX, J 142: DOUBLE PRECISION BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21, 143: $ CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21, 144: $ LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R, 145: $ UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S, 146: $ UR22, XI1, XI2, XR1, XR2 147: * .. 148: * .. Local Arrays .. 149: LOGICAL RSWAP( 4 ), ZSWAP( 4 ) 150: INTEGER IPIVOT( 4, 4 ) 151: DOUBLE PRECISION CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 ) 152: * .. 153: * .. External Functions .. 154: DOUBLE PRECISION DLAMCH 155: EXTERNAL DLAMCH 156: * .. 157: * .. External Subroutines .. 158: EXTERNAL DLADIV 159: * .. 160: * .. Intrinsic Functions .. 161: INTRINSIC ABS, MAX 162: * .. 163: * .. Equivalences .. 164: EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ), 165: $ ( CR( 1, 1 ), CRV( 1 ) ) 166: * .. 167: * .. Data statements .. 168: DATA ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. / 169: DATA RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. / 170: DATA IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4, 171: $ 3, 2, 1 / 172: * .. 173: * .. Executable Statements .. 174: * 175: * Compute BIGNUM 176: * 177: SMLNUM = TWO*DLAMCH( 'Safe minimum' ) 178: BIGNUM = ONE / SMLNUM 179: SMINI = MAX( SMIN, SMLNUM ) 180: * 181: * Don't check for input errors 182: * 183: INFO = 0 184: * 185: * Standard Initializations 186: * 187: SCALE = ONE 188: * 189: IF( NA.EQ.1 ) THEN 190: * 191: * 1 x 1 (i.e., scalar) system C X = B 192: * 193: IF( NW.EQ.1 ) THEN 194: * 195: * Real 1x1 system. 196: * 197: * C = ca A - w D 198: * 199: CSR = CA*A( 1, 1 ) - WR*D1 200: CNORM = ABS( CSR ) 201: * 202: * If | C | < SMINI, use C = SMINI 203: * 204: IF( CNORM.LT.SMINI ) THEN 205: CSR = SMINI 206: CNORM = SMINI 207: INFO = 1 208: END IF 209: * 210: * Check scaling for X = B / C 211: * 212: BNORM = ABS( B( 1, 1 ) ) 213: IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN 214: IF( BNORM.GT.BIGNUM*CNORM ) 215: $ SCALE = ONE / BNORM 216: END IF 217: * 218: * Compute X 219: * 220: X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR 221: XNORM = ABS( X( 1, 1 ) ) 222: ELSE 223: * 224: * Complex 1x1 system (w is complex) 225: * 226: * C = ca A - w D 227: * 228: CSR = CA*A( 1, 1 ) - WR*D1 229: CSI = -WI*D1 230: CNORM = ABS( CSR ) + ABS( CSI ) 231: * 232: * If | C | < SMINI, use C = SMINI 233: * 234: IF( CNORM.LT.SMINI ) THEN 235: CSR = SMINI 236: CSI = ZERO 237: CNORM = SMINI 238: INFO = 1 239: END IF 240: * 241: * Check scaling for X = B / C 242: * 243: BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) ) 244: IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN 245: IF( BNORM.GT.BIGNUM*CNORM ) 246: $ SCALE = ONE / BNORM 247: END IF 248: * 249: * Compute X 250: * 251: CALL DLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI, 252: $ X( 1, 1 ), X( 1, 2 ) ) 253: XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) ) 254: END IF 255: * 256: ELSE 257: * 258: * 2x2 System 259: * 260: * Compute the real part of C = ca A - w D (or ca A' - w D ) 261: * 262: CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1 263: CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2 264: IF( LTRANS ) THEN 265: CR( 1, 2 ) = CA*A( 2, 1 ) 266: CR( 2, 1 ) = CA*A( 1, 2 ) 267: ELSE 268: CR( 2, 1 ) = CA*A( 2, 1 ) 269: CR( 1, 2 ) = CA*A( 1, 2 ) 270: END IF 271: * 272: IF( NW.EQ.1 ) THEN 273: * 274: * Real 2x2 system (w is real) 275: * 276: * Find the largest element in C 277: * 278: CMAX = ZERO 279: ICMAX = 0 280: * 281: DO 10 J = 1, 4 282: IF( ABS( CRV( J ) ).GT.CMAX ) THEN 283: CMAX = ABS( CRV( J ) ) 284: ICMAX = J 285: END IF 286: 10 CONTINUE 287: * 288: * If norm(C) < SMINI, use SMINI*identity. 289: * 290: IF( CMAX.LT.SMINI ) THEN 291: BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) ) 292: IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN 293: IF( BNORM.GT.BIGNUM*SMINI ) 294: $ SCALE = ONE / BNORM 295: END IF 296: TEMP = SCALE / SMINI 297: X( 1, 1 ) = TEMP*B( 1, 1 ) 298: X( 2, 1 ) = TEMP*B( 2, 1 ) 299: XNORM = TEMP*BNORM 300: INFO = 1 301: RETURN 302: END IF 303: * 304: * Gaussian elimination with complete pivoting. 305: * 306: UR11 = CRV( ICMAX ) 307: CR21 = CRV( IPIVOT( 2, ICMAX ) ) 308: UR12 = CRV( IPIVOT( 3, ICMAX ) ) 309: CR22 = CRV( IPIVOT( 4, ICMAX ) ) 310: UR11R = ONE / UR11 311: LR21 = UR11R*CR21 312: UR22 = CR22 - UR12*LR21 313: * 314: * If smaller pivot < SMINI, use SMINI 315: * 316: IF( ABS( UR22 ).LT.SMINI ) THEN 317: UR22 = SMINI 318: INFO = 1 319: END IF 320: IF( RSWAP( ICMAX ) ) THEN 321: BR1 = B( 2, 1 ) 322: BR2 = B( 1, 1 ) 323: ELSE 324: BR1 = B( 1, 1 ) 325: BR2 = B( 2, 1 ) 326: END IF 327: BR2 = BR2 - LR21*BR1 328: BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) ) 329: IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN 330: IF( BBND.GE.BIGNUM*ABS( UR22 ) ) 331: $ SCALE = ONE / BBND 332: END IF 333: * 334: XR2 = ( BR2*SCALE ) / UR22 335: XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 ) 336: IF( ZSWAP( ICMAX ) ) THEN 337: X( 1, 1 ) = XR2 338: X( 2, 1 ) = XR1 339: ELSE 340: X( 1, 1 ) = XR1 341: X( 2, 1 ) = XR2 342: END IF 343: XNORM = MAX( ABS( XR1 ), ABS( XR2 ) ) 344: * 345: * Further scaling if norm(A) norm(X) > overflow 346: * 347: IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN 348: IF( XNORM.GT.BIGNUM / CMAX ) THEN 349: TEMP = CMAX / BIGNUM 350: X( 1, 1 ) = TEMP*X( 1, 1 ) 351: X( 2, 1 ) = TEMP*X( 2, 1 ) 352: XNORM = TEMP*XNORM 353: SCALE = TEMP*SCALE 354: END IF 355: END IF 356: ELSE 357: * 358: * Complex 2x2 system (w is complex) 359: * 360: * Find the largest element in C 361: * 362: CI( 1, 1 ) = -WI*D1 363: CI( 2, 1 ) = ZERO 364: CI( 1, 2 ) = ZERO 365: CI( 2, 2 ) = -WI*D2 366: CMAX = ZERO 367: ICMAX = 0 368: * 369: DO 20 J = 1, 4 370: IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN 371: CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) ) 372: ICMAX = J 373: END IF 374: 20 CONTINUE 375: * 376: * If norm(C) < SMINI, use SMINI*identity. 377: * 378: IF( CMAX.LT.SMINI ) THEN 379: BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ), 380: $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) ) 381: IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN 382: IF( BNORM.GT.BIGNUM*SMINI ) 383: $ SCALE = ONE / BNORM 384: END IF 385: TEMP = SCALE / SMINI 386: X( 1, 1 ) = TEMP*B( 1, 1 ) 387: X( 2, 1 ) = TEMP*B( 2, 1 ) 388: X( 1, 2 ) = TEMP*B( 1, 2 ) 389: X( 2, 2 ) = TEMP*B( 2, 2 ) 390: XNORM = TEMP*BNORM 391: INFO = 1 392: RETURN 393: END IF 394: * 395: * Gaussian elimination with complete pivoting. 396: * 397: UR11 = CRV( ICMAX ) 398: UI11 = CIV( ICMAX ) 399: CR21 = CRV( IPIVOT( 2, ICMAX ) ) 400: CI21 = CIV( IPIVOT( 2, ICMAX ) ) 401: UR12 = CRV( IPIVOT( 3, ICMAX ) ) 402: UI12 = CIV( IPIVOT( 3, ICMAX ) ) 403: CR22 = CRV( IPIVOT( 4, ICMAX ) ) 404: CI22 = CIV( IPIVOT( 4, ICMAX ) ) 405: IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN 406: * 407: * Code when off-diagonals of pivoted C are real 408: * 409: IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN 410: TEMP = UI11 / UR11 411: UR11R = ONE / ( UR11*( ONE+TEMP**2 ) ) 412: UI11R = -TEMP*UR11R 413: ELSE 414: TEMP = UR11 / UI11 415: UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) ) 416: UR11R = -TEMP*UI11R 417: END IF 418: LR21 = CR21*UR11R 419: LI21 = CR21*UI11R 420: UR12S = UR12*UR11R 421: UI12S = UR12*UI11R 422: UR22 = CR22 - UR12*LR21 423: UI22 = CI22 - UR12*LI21 424: ELSE 425: * 426: * Code when diagonals of pivoted C are real 427: * 428: UR11R = ONE / UR11 429: UI11R = ZERO 430: LR21 = CR21*UR11R 431: LI21 = CI21*UR11R 432: UR12S = UR12*UR11R 433: UI12S = UI12*UR11R 434: UR22 = CR22 - UR12*LR21 + UI12*LI21 435: UI22 = -UR12*LI21 - UI12*LR21 436: END IF 437: U22ABS = ABS( UR22 ) + ABS( UI22 ) 438: * 439: * If smaller pivot < SMINI, use SMINI 440: * 441: IF( U22ABS.LT.SMINI ) THEN 442: UR22 = SMINI 443: UI22 = ZERO 444: INFO = 1 445: END IF 446: IF( RSWAP( ICMAX ) ) THEN 447: BR2 = B( 1, 1 ) 448: BR1 = B( 2, 1 ) 449: BI2 = B( 1, 2 ) 450: BI1 = B( 2, 2 ) 451: ELSE 452: BR1 = B( 1, 1 ) 453: BR2 = B( 2, 1 ) 454: BI1 = B( 1, 2 ) 455: BI2 = B( 2, 2 ) 456: END IF 457: BR2 = BR2 - LR21*BR1 + LI21*BI1 458: BI2 = BI2 - LI21*BR1 - LR21*BI1 459: BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )* 460: $ ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ), 461: $ ABS( BR2 )+ABS( BI2 ) ) 462: IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN 463: IF( BBND.GE.BIGNUM*U22ABS ) THEN 464: SCALE = ONE / BBND 465: BR1 = SCALE*BR1 466: BI1 = SCALE*BI1 467: BR2 = SCALE*BR2 468: BI2 = SCALE*BI2 469: END IF 470: END IF 471: * 472: CALL DLADIV( BR2, BI2, UR22, UI22, XR2, XI2 ) 473: XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2 474: XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2 475: IF( ZSWAP( ICMAX ) ) THEN 476: X( 1, 1 ) = XR2 477: X( 2, 1 ) = XR1 478: X( 1, 2 ) = XI2 479: X( 2, 2 ) = XI1 480: ELSE 481: X( 1, 1 ) = XR1 482: X( 2, 1 ) = XR2 483: X( 1, 2 ) = XI1 484: X( 2, 2 ) = XI2 485: END IF 486: XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) ) 487: * 488: * Further scaling if norm(A) norm(X) > overflow 489: * 490: IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN 491: IF( XNORM.GT.BIGNUM / CMAX ) THEN 492: TEMP = CMAX / BIGNUM 493: X( 1, 1 ) = TEMP*X( 1, 1 ) 494: X( 2, 1 ) = TEMP*X( 2, 1 ) 495: X( 1, 2 ) = TEMP*X( 1, 2 ) 496: X( 2, 2 ) = TEMP*X( 2, 2 ) 497: XNORM = TEMP*XNORM 498: SCALE = TEMP*SCALE 499: END IF 500: END IF 501: END IF 502: END IF 503: * 504: RETURN 505: * 506: * End of DLALN2 507: * 508: END