Annotation of rpl/lapack/lapack/dlaln2.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaln2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack