Annotation of rpl/lapack/lapack/dlasy2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
! 2: $ LDTR, B, LDB, SCALE, X, LDX, 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 LTRANL, LTRANR
! 11: INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
! 12: DOUBLE PRECISION SCALE, XNORM
! 13: * ..
! 14: * .. Array Arguments ..
! 15: DOUBLE PRECISION B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
! 16: $ X( LDX, * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
! 23: *
! 24: * op(TL)*X + ISGN*X*op(TR) = SCALE*B,
! 25: *
! 26: * where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
! 27: * -1. op(T) = T or T', where T' denotes the transpose of T.
! 28: *
! 29: * Arguments
! 30: * =========
! 31: *
! 32: * LTRANL (input) LOGICAL
! 33: * On entry, LTRANL specifies the op(TL):
! 34: * = .FALSE., op(TL) = TL,
! 35: * = .TRUE., op(TL) = TL'.
! 36: *
! 37: * LTRANR (input) LOGICAL
! 38: * On entry, LTRANR specifies the op(TR):
! 39: * = .FALSE., op(TR) = TR,
! 40: * = .TRUE., op(TR) = TR'.
! 41: *
! 42: * ISGN (input) INTEGER
! 43: * On entry, ISGN specifies the sign of the equation
! 44: * as described before. ISGN may only be 1 or -1.
! 45: *
! 46: * N1 (input) INTEGER
! 47: * On entry, N1 specifies the order of matrix TL.
! 48: * N1 may only be 0, 1 or 2.
! 49: *
! 50: * N2 (input) INTEGER
! 51: * On entry, N2 specifies the order of matrix TR.
! 52: * N2 may only be 0, 1 or 2.
! 53: *
! 54: * TL (input) DOUBLE PRECISION array, dimension (LDTL,2)
! 55: * On entry, TL contains an N1 by N1 matrix.
! 56: *
! 57: * LDTL (input) INTEGER
! 58: * The leading dimension of the matrix TL. LDTL >= max(1,N1).
! 59: *
! 60: * TR (input) DOUBLE PRECISION array, dimension (LDTR,2)
! 61: * On entry, TR contains an N2 by N2 matrix.
! 62: *
! 63: * LDTR (input) INTEGER
! 64: * The leading dimension of the matrix TR. LDTR >= max(1,N2).
! 65: *
! 66: * B (input) DOUBLE PRECISION array, dimension (LDB,2)
! 67: * On entry, the N1 by N2 matrix B contains the right-hand
! 68: * side of the equation.
! 69: *
! 70: * LDB (input) INTEGER
! 71: * The leading dimension of the matrix B. LDB >= max(1,N1).
! 72: *
! 73: * SCALE (output) DOUBLE PRECISION
! 74: * On exit, SCALE contains the scale factor. SCALE is chosen
! 75: * less than or equal to 1 to prevent the solution overflowing.
! 76: *
! 77: * X (output) DOUBLE PRECISION array, dimension (LDX,2)
! 78: * On exit, X contains the N1 by N2 solution.
! 79: *
! 80: * LDX (input) INTEGER
! 81: * The leading dimension of the matrix X. LDX >= max(1,N1).
! 82: *
! 83: * XNORM (output) DOUBLE PRECISION
! 84: * On exit, XNORM is the infinity-norm of the solution.
! 85: *
! 86: * INFO (output) INTEGER
! 87: * On exit, INFO is set to
! 88: * 0: successful exit.
! 89: * 1: TL and TR have too close eigenvalues, so TL or
! 90: * TR is perturbed to get a nonsingular equation.
! 91: * NOTE: In the interests of speed, this routine does not
! 92: * check the inputs for errors.
! 93: *
! 94: * =====================================================================
! 95: *
! 96: * .. Parameters ..
! 97: DOUBLE PRECISION ZERO, ONE
! 98: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 99: DOUBLE PRECISION TWO, HALF, EIGHT
! 100: PARAMETER ( TWO = 2.0D+0, HALF = 0.5D+0, EIGHT = 8.0D+0 )
! 101: * ..
! 102: * .. Local Scalars ..
! 103: LOGICAL BSWAP, XSWAP
! 104: INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K
! 105: DOUBLE PRECISION BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
! 106: $ TEMP, U11, U12, U22, XMAX
! 107: * ..
! 108: * .. Local Arrays ..
! 109: LOGICAL BSWPIV( 4 ), XSWPIV( 4 )
! 110: INTEGER JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
! 111: $ LOCU22( 4 )
! 112: DOUBLE PRECISION BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
! 113: * ..
! 114: * .. External Functions ..
! 115: INTEGER IDAMAX
! 116: DOUBLE PRECISION DLAMCH
! 117: EXTERNAL IDAMAX, DLAMCH
! 118: * ..
! 119: * .. External Subroutines ..
! 120: EXTERNAL DCOPY, DSWAP
! 121: * ..
! 122: * .. Intrinsic Functions ..
! 123: INTRINSIC ABS, MAX
! 124: * ..
! 125: * .. Data statements ..
! 126: DATA LOCU12 / 3, 4, 1, 2 / , LOCL21 / 2, 1, 4, 3 / ,
! 127: $ LOCU22 / 4, 3, 2, 1 /
! 128: DATA XSWPIV / .FALSE., .FALSE., .TRUE., .TRUE. /
! 129: DATA BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. /
! 130: * ..
! 131: * .. Executable Statements ..
! 132: *
! 133: * Do not check the input parameters for errors
! 134: *
! 135: INFO = 0
! 136: *
! 137: * Quick return if possible
! 138: *
! 139: IF( N1.EQ.0 .OR. N2.EQ.0 )
! 140: $ RETURN
! 141: *
! 142: * Set constants to control overflow
! 143: *
! 144: EPS = DLAMCH( 'P' )
! 145: SMLNUM = DLAMCH( 'S' ) / EPS
! 146: SGN = ISGN
! 147: *
! 148: K = N1 + N1 + N2 - 2
! 149: GO TO ( 10, 20, 30, 50 )K
! 150: *
! 151: * 1 by 1: TL11*X + SGN*X*TR11 = B11
! 152: *
! 153: 10 CONTINUE
! 154: TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 )
! 155: BET = ABS( TAU1 )
! 156: IF( BET.LE.SMLNUM ) THEN
! 157: TAU1 = SMLNUM
! 158: BET = SMLNUM
! 159: INFO = 1
! 160: END IF
! 161: *
! 162: SCALE = ONE
! 163: GAM = ABS( B( 1, 1 ) )
! 164: IF( SMLNUM*GAM.GT.BET )
! 165: $ SCALE = ONE / GAM
! 166: *
! 167: X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1
! 168: XNORM = ABS( X( 1, 1 ) )
! 169: RETURN
! 170: *
! 171: * 1 by 2:
! 172: * TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12]
! 173: * [TR21 TR22]
! 174: *
! 175: 20 CONTINUE
! 176: *
! 177: SMIN = MAX( EPS*MAX( ABS( TL( 1, 1 ) ), ABS( TR( 1, 1 ) ),
! 178: $ ABS( TR( 1, 2 ) ), ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ),
! 179: $ SMLNUM )
! 180: TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
! 181: TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
! 182: IF( LTRANR ) THEN
! 183: TMP( 2 ) = SGN*TR( 2, 1 )
! 184: TMP( 3 ) = SGN*TR( 1, 2 )
! 185: ELSE
! 186: TMP( 2 ) = SGN*TR( 1, 2 )
! 187: TMP( 3 ) = SGN*TR( 2, 1 )
! 188: END IF
! 189: BTMP( 1 ) = B( 1, 1 )
! 190: BTMP( 2 ) = B( 1, 2 )
! 191: GO TO 40
! 192: *
! 193: * 2 by 1:
! 194: * op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11]
! 195: * [TL21 TL22] [X21] [X21] [B21]
! 196: *
! 197: 30 CONTINUE
! 198: SMIN = MAX( EPS*MAX( ABS( TR( 1, 1 ) ), ABS( TL( 1, 1 ) ),
! 199: $ ABS( TL( 1, 2 ) ), ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ),
! 200: $ SMLNUM )
! 201: TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
! 202: TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
! 203: IF( LTRANL ) THEN
! 204: TMP( 2 ) = TL( 1, 2 )
! 205: TMP( 3 ) = TL( 2, 1 )
! 206: ELSE
! 207: TMP( 2 ) = TL( 2, 1 )
! 208: TMP( 3 ) = TL( 1, 2 )
! 209: END IF
! 210: BTMP( 1 ) = B( 1, 1 )
! 211: BTMP( 2 ) = B( 2, 1 )
! 212: 40 CONTINUE
! 213: *
! 214: * Solve 2 by 2 system using complete pivoting.
! 215: * Set pivots less than SMIN to SMIN.
! 216: *
! 217: IPIV = IDAMAX( 4, TMP, 1 )
! 218: U11 = TMP( IPIV )
! 219: IF( ABS( U11 ).LE.SMIN ) THEN
! 220: INFO = 1
! 221: U11 = SMIN
! 222: END IF
! 223: U12 = TMP( LOCU12( IPIV ) )
! 224: L21 = TMP( LOCL21( IPIV ) ) / U11
! 225: U22 = TMP( LOCU22( IPIV ) ) - U12*L21
! 226: XSWAP = XSWPIV( IPIV )
! 227: BSWAP = BSWPIV( IPIV )
! 228: IF( ABS( U22 ).LE.SMIN ) THEN
! 229: INFO = 1
! 230: U22 = SMIN
! 231: END IF
! 232: IF( BSWAP ) THEN
! 233: TEMP = BTMP( 2 )
! 234: BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
! 235: BTMP( 1 ) = TEMP
! 236: ELSE
! 237: BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
! 238: END IF
! 239: SCALE = ONE
! 240: IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR.
! 241: $ ( TWO*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) THEN
! 242: SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) )
! 243: BTMP( 1 ) = BTMP( 1 )*SCALE
! 244: BTMP( 2 ) = BTMP( 2 )*SCALE
! 245: END IF
! 246: X2( 2 ) = BTMP( 2 ) / U22
! 247: X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 )
! 248: IF( XSWAP ) THEN
! 249: TEMP = X2( 2 )
! 250: X2( 2 ) = X2( 1 )
! 251: X2( 1 ) = TEMP
! 252: END IF
! 253: X( 1, 1 ) = X2( 1 )
! 254: IF( N1.EQ.1 ) THEN
! 255: X( 1, 2 ) = X2( 2 )
! 256: XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
! 257: ELSE
! 258: X( 2, 1 ) = X2( 2 )
! 259: XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) )
! 260: END IF
! 261: RETURN
! 262: *
! 263: * 2 by 2:
! 264: * op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
! 265: * [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22]
! 266: *
! 267: * Solve equivalent 4 by 4 system using complete pivoting.
! 268: * Set pivots less than SMIN to SMIN.
! 269: *
! 270: 50 CONTINUE
! 271: SMIN = MAX( ABS( TR( 1, 1 ) ), ABS( TR( 1, 2 ) ),
! 272: $ ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) )
! 273: SMIN = MAX( SMIN, ABS( TL( 1, 1 ) ), ABS( TL( 1, 2 ) ),
! 274: $ ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) )
! 275: SMIN = MAX( EPS*SMIN, SMLNUM )
! 276: BTMP( 1 ) = ZERO
! 277: CALL DCOPY( 16, BTMP, 0, T16, 1 )
! 278: T16( 1, 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
! 279: T16( 2, 2 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
! 280: T16( 3, 3 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
! 281: T16( 4, 4 ) = TL( 2, 2 ) + SGN*TR( 2, 2 )
! 282: IF( LTRANL ) THEN
! 283: T16( 1, 2 ) = TL( 2, 1 )
! 284: T16( 2, 1 ) = TL( 1, 2 )
! 285: T16( 3, 4 ) = TL( 2, 1 )
! 286: T16( 4, 3 ) = TL( 1, 2 )
! 287: ELSE
! 288: T16( 1, 2 ) = TL( 1, 2 )
! 289: T16( 2, 1 ) = TL( 2, 1 )
! 290: T16( 3, 4 ) = TL( 1, 2 )
! 291: T16( 4, 3 ) = TL( 2, 1 )
! 292: END IF
! 293: IF( LTRANR ) THEN
! 294: T16( 1, 3 ) = SGN*TR( 1, 2 )
! 295: T16( 2, 4 ) = SGN*TR( 1, 2 )
! 296: T16( 3, 1 ) = SGN*TR( 2, 1 )
! 297: T16( 4, 2 ) = SGN*TR( 2, 1 )
! 298: ELSE
! 299: T16( 1, 3 ) = SGN*TR( 2, 1 )
! 300: T16( 2, 4 ) = SGN*TR( 2, 1 )
! 301: T16( 3, 1 ) = SGN*TR( 1, 2 )
! 302: T16( 4, 2 ) = SGN*TR( 1, 2 )
! 303: END IF
! 304: BTMP( 1 ) = B( 1, 1 )
! 305: BTMP( 2 ) = B( 2, 1 )
! 306: BTMP( 3 ) = B( 1, 2 )
! 307: BTMP( 4 ) = B( 2, 2 )
! 308: *
! 309: * Perform elimination
! 310: *
! 311: DO 100 I = 1, 3
! 312: XMAX = ZERO
! 313: DO 70 IP = I, 4
! 314: DO 60 JP = I, 4
! 315: IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN
! 316: XMAX = ABS( T16( IP, JP ) )
! 317: IPSV = IP
! 318: JPSV = JP
! 319: END IF
! 320: 60 CONTINUE
! 321: 70 CONTINUE
! 322: IF( IPSV.NE.I ) THEN
! 323: CALL DSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 )
! 324: TEMP = BTMP( I )
! 325: BTMP( I ) = BTMP( IPSV )
! 326: BTMP( IPSV ) = TEMP
! 327: END IF
! 328: IF( JPSV.NE.I )
! 329: $ CALL DSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
! 330: JPIV( I ) = JPSV
! 331: IF( ABS( T16( I, I ) ).LT.SMIN ) THEN
! 332: INFO = 1
! 333: T16( I, I ) = SMIN
! 334: END IF
! 335: DO 90 J = I + 1, 4
! 336: T16( J, I ) = T16( J, I ) / T16( I, I )
! 337: BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I )
! 338: DO 80 K = I + 1, 4
! 339: T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K )
! 340: 80 CONTINUE
! 341: 90 CONTINUE
! 342: 100 CONTINUE
! 343: IF( ABS( T16( 4, 4 ) ).LT.SMIN )
! 344: $ T16( 4, 4 ) = SMIN
! 345: SCALE = ONE
! 346: IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
! 347: $ ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
! 348: $ ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
! 349: $ ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
! 350: SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ),
! 351: $ ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
! 352: BTMP( 1 ) = BTMP( 1 )*SCALE
! 353: BTMP( 2 ) = BTMP( 2 )*SCALE
! 354: BTMP( 3 ) = BTMP( 3 )*SCALE
! 355: BTMP( 4 ) = BTMP( 4 )*SCALE
! 356: END IF
! 357: DO 120 I = 1, 4
! 358: K = 5 - I
! 359: TEMP = ONE / T16( K, K )
! 360: TMP( K ) = BTMP( K )*TEMP
! 361: DO 110 J = K + 1, 4
! 362: TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J )
! 363: 110 CONTINUE
! 364: 120 CONTINUE
! 365: DO 130 I = 1, 3
! 366: IF( JPIV( 4-I ).NE.4-I ) THEN
! 367: TEMP = TMP( 4-I )
! 368: TMP( 4-I ) = TMP( JPIV( 4-I ) )
! 369: TMP( JPIV( 4-I ) ) = TEMP
! 370: END IF
! 371: 130 CONTINUE
! 372: X( 1, 1 ) = TMP( 1 )
! 373: X( 2, 1 ) = TMP( 2 )
! 374: X( 1, 2 ) = TMP( 3 )
! 375: X( 2, 2 ) = TMP( 4 )
! 376: XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ),
! 377: $ ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) )
! 378: RETURN
! 379: *
! 380: * End of DLASY2
! 381: *
! 382: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasy2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack