Annotation of rpl/lapack/lapack/dlanv2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
! 2: *
! 3: * -- LAPACK driver routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: DOUBLE PRECISION A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
! 10: * ..
! 11: *
! 12: * Purpose
! 13: * =======
! 14: *
! 15: * DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
! 16: * matrix in standard form:
! 17: *
! 18: * [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ]
! 19: * [ C D ] [ SN CS ] [ CC DD ] [-SN CS ]
! 20: *
! 21: * where either
! 22: * 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
! 23: * 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
! 24: * conjugate eigenvalues.
! 25: *
! 26: * Arguments
! 27: * =========
! 28: *
! 29: * A (input/output) DOUBLE PRECISION
! 30: * B (input/output) DOUBLE PRECISION
! 31: * C (input/output) DOUBLE PRECISION
! 32: * D (input/output) DOUBLE PRECISION
! 33: * On entry, the elements of the input matrix.
! 34: * On exit, they are overwritten by the elements of the
! 35: * standardised Schur form.
! 36: *
! 37: * RT1R (output) DOUBLE PRECISION
! 38: * RT1I (output) DOUBLE PRECISION
! 39: * RT2R (output) DOUBLE PRECISION
! 40: * RT2I (output) DOUBLE PRECISION
! 41: * The real and imaginary parts of the eigenvalues. If the
! 42: * eigenvalues are a complex conjugate pair, RT1I > 0.
! 43: *
! 44: * CS (output) DOUBLE PRECISION
! 45: * SN (output) DOUBLE PRECISION
! 46: * Parameters of the rotation matrix.
! 47: *
! 48: * Further Details
! 49: * ===============
! 50: *
! 51: * Modified by V. Sima, Research Institute for Informatics, Bucharest,
! 52: * Romania, to reduce the risk of cancellation errors,
! 53: * when computing real eigenvalues, and to ensure, if possible, that
! 54: * abs(RT1R) >= abs(RT2R).
! 55: *
! 56: * =====================================================================
! 57: *
! 58: * .. Parameters ..
! 59: DOUBLE PRECISION ZERO, HALF, ONE
! 60: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
! 61: DOUBLE PRECISION MULTPL
! 62: PARAMETER ( MULTPL = 4.0D+0 )
! 63: * ..
! 64: * .. Local Scalars ..
! 65: DOUBLE PRECISION AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB,
! 66: $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z
! 67: * ..
! 68: * .. External Functions ..
! 69: DOUBLE PRECISION DLAMCH, DLAPY2
! 70: EXTERNAL DLAMCH, DLAPY2
! 71: * ..
! 72: * .. Intrinsic Functions ..
! 73: INTRINSIC ABS, MAX, MIN, SIGN, SQRT
! 74: * ..
! 75: * .. Executable Statements ..
! 76: *
! 77: EPS = DLAMCH( 'P' )
! 78: IF( C.EQ.ZERO ) THEN
! 79: CS = ONE
! 80: SN = ZERO
! 81: GO TO 10
! 82: *
! 83: ELSE IF( B.EQ.ZERO ) THEN
! 84: *
! 85: * Swap rows and columns
! 86: *
! 87: CS = ZERO
! 88: SN = ONE
! 89: TEMP = D
! 90: D = A
! 91: A = TEMP
! 92: B = -C
! 93: C = ZERO
! 94: GO TO 10
! 95: ELSE IF( ( A-D ).EQ.ZERO .AND. SIGN( ONE, B ).NE.SIGN( ONE, C ) )
! 96: $ THEN
! 97: CS = ONE
! 98: SN = ZERO
! 99: GO TO 10
! 100: ELSE
! 101: *
! 102: TEMP = A - D
! 103: P = HALF*TEMP
! 104: BCMAX = MAX( ABS( B ), ABS( C ) )
! 105: BCMIS = MIN( ABS( B ), ABS( C ) )*SIGN( ONE, B )*SIGN( ONE, C )
! 106: SCALE = MAX( ABS( P ), BCMAX )
! 107: Z = ( P / SCALE )*P + ( BCMAX / SCALE )*BCMIS
! 108: *
! 109: * If Z is of the order of the machine accuracy, postpone the
! 110: * decision on the nature of eigenvalues
! 111: *
! 112: IF( Z.GE.MULTPL*EPS ) THEN
! 113: *
! 114: * Real eigenvalues. Compute A and D.
! 115: *
! 116: Z = P + SIGN( SQRT( SCALE )*SQRT( Z ), P )
! 117: A = D + Z
! 118: D = D - ( BCMAX / Z )*BCMIS
! 119: *
! 120: * Compute B and the rotation matrix
! 121: *
! 122: TAU = DLAPY2( C, Z )
! 123: CS = Z / TAU
! 124: SN = C / TAU
! 125: B = B - C
! 126: C = ZERO
! 127: ELSE
! 128: *
! 129: * Complex eigenvalues, or real (almost) equal eigenvalues.
! 130: * Make diagonal elements equal.
! 131: *
! 132: SIGMA = B + C
! 133: TAU = DLAPY2( SIGMA, TEMP )
! 134: CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) )
! 135: SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA )
! 136: *
! 137: * Compute [ AA BB ] = [ A B ] [ CS -SN ]
! 138: * [ CC DD ] [ C D ] [ SN CS ]
! 139: *
! 140: AA = A*CS + B*SN
! 141: BB = -A*SN + B*CS
! 142: CC = C*CS + D*SN
! 143: DD = -C*SN + D*CS
! 144: *
! 145: * Compute [ A B ] = [ CS SN ] [ AA BB ]
! 146: * [ C D ] [-SN CS ] [ CC DD ]
! 147: *
! 148: A = AA*CS + CC*SN
! 149: B = BB*CS + DD*SN
! 150: C = -AA*SN + CC*CS
! 151: D = -BB*SN + DD*CS
! 152: *
! 153: TEMP = HALF*( A+D )
! 154: A = TEMP
! 155: D = TEMP
! 156: *
! 157: IF( C.NE.ZERO ) THEN
! 158: IF( B.NE.ZERO ) THEN
! 159: IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN
! 160: *
! 161: * Real eigenvalues: reduce to upper triangular form
! 162: *
! 163: SAB = SQRT( ABS( B ) )
! 164: SAC = SQRT( ABS( C ) )
! 165: P = SIGN( SAB*SAC, C )
! 166: TAU = ONE / SQRT( ABS( B+C ) )
! 167: A = TEMP + P
! 168: D = TEMP - P
! 169: B = B - C
! 170: C = ZERO
! 171: CS1 = SAB*TAU
! 172: SN1 = SAC*TAU
! 173: TEMP = CS*CS1 - SN*SN1
! 174: SN = CS*SN1 + SN*CS1
! 175: CS = TEMP
! 176: END IF
! 177: ELSE
! 178: B = -C
! 179: C = ZERO
! 180: TEMP = CS
! 181: CS = -SN
! 182: SN = TEMP
! 183: END IF
! 184: END IF
! 185: END IF
! 186: *
! 187: END IF
! 188: *
! 189: 10 CONTINUE
! 190: *
! 191: * Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
! 192: *
! 193: RT1R = A
! 194: RT2R = D
! 195: IF( C.EQ.ZERO ) THEN
! 196: RT1I = ZERO
! 197: RT2I = ZERO
! 198: ELSE
! 199: RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) )
! 200: RT2I = -RT1I
! 201: END IF
! 202: RETURN
! 203: *
! 204: * End of DLANV2
! 205: *
! 206: END
CVSweb interface <joel.bertrand@systella.fr>