Annotation of rpl/lapack/lapack/dlasd5.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
! 2: *
! 3: * -- LAPACK auxiliary 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: INTEGER I
! 10: DOUBLE PRECISION DSIGMA, RHO
! 11: * ..
! 12: * .. Array Arguments ..
! 13: DOUBLE PRECISION D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 )
! 14: * ..
! 15: *
! 16: * Purpose
! 17: * =======
! 18: *
! 19: * This subroutine computes the square root of the I-th eigenvalue
! 20: * of a positive symmetric rank-one modification of a 2-by-2 diagonal
! 21: * matrix
! 22: *
! 23: * diag( D ) * diag( D ) + RHO * Z * transpose(Z) .
! 24: *
! 25: * The diagonal entries in the array D are assumed to satisfy
! 26: *
! 27: * 0 <= D(i) < D(j) for i < j .
! 28: *
! 29: * We also assume RHO > 0 and that the Euclidean norm of the vector
! 30: * Z is one.
! 31: *
! 32: * Arguments
! 33: * =========
! 34: *
! 35: * I (input) INTEGER
! 36: * The index of the eigenvalue to be computed. I = 1 or I = 2.
! 37: *
! 38: * D (input) DOUBLE PRECISION array, dimension ( 2 )
! 39: * The original eigenvalues. We assume 0 <= D(1) < D(2).
! 40: *
! 41: * Z (input) DOUBLE PRECISION array, dimension ( 2 )
! 42: * The components of the updating vector.
! 43: *
! 44: * DELTA (output) DOUBLE PRECISION array, dimension ( 2 )
! 45: * Contains (D(j) - sigma_I) in its j-th component.
! 46: * The vector DELTA contains the information necessary
! 47: * to construct the eigenvectors.
! 48: *
! 49: * RHO (input) DOUBLE PRECISION
! 50: * The scalar in the symmetric updating formula.
! 51: *
! 52: * DSIGMA (output) DOUBLE PRECISION
! 53: * The computed sigma_I, the I-th updated eigenvalue.
! 54: *
! 55: * WORK (workspace) DOUBLE PRECISION array, dimension ( 2 )
! 56: * WORK contains (D(j) + sigma_I) in its j-th component.
! 57: *
! 58: * Further Details
! 59: * ===============
! 60: *
! 61: * Based on contributions by
! 62: * Ren-Cang Li, Computer Science Division, University of California
! 63: * at Berkeley, USA
! 64: *
! 65: * =====================================================================
! 66: *
! 67: * .. Parameters ..
! 68: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR
! 69: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
! 70: $ THREE = 3.0D+0, FOUR = 4.0D+0 )
! 71: * ..
! 72: * .. Local Scalars ..
! 73: DOUBLE PRECISION B, C, DEL, DELSQ, TAU, W
! 74: * ..
! 75: * .. Intrinsic Functions ..
! 76: INTRINSIC ABS, SQRT
! 77: * ..
! 78: * .. Executable Statements ..
! 79: *
! 80: DEL = D( 2 ) - D( 1 )
! 81: DELSQ = DEL*( D( 2 )+D( 1 ) )
! 82: IF( I.EQ.1 ) THEN
! 83: W = ONE + FOUR*RHO*( Z( 2 )*Z( 2 ) / ( D( 1 )+THREE*D( 2 ) )-
! 84: $ Z( 1 )*Z( 1 ) / ( THREE*D( 1 )+D( 2 ) ) ) / DEL
! 85: IF( W.GT.ZERO ) THEN
! 86: B = DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
! 87: C = RHO*Z( 1 )*Z( 1 )*DELSQ
! 88: *
! 89: * B > ZERO, always
! 90: *
! 91: * The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 )
! 92: *
! 93: TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) )
! 94: *
! 95: * The following TAU is DSIGMA - D( 1 )
! 96: *
! 97: TAU = TAU / ( D( 1 )+SQRT( D( 1 )*D( 1 )+TAU ) )
! 98: DSIGMA = D( 1 ) + TAU
! 99: DELTA( 1 ) = -TAU
! 100: DELTA( 2 ) = DEL - TAU
! 101: WORK( 1 ) = TWO*D( 1 ) + TAU
! 102: WORK( 2 ) = ( D( 1 )+TAU ) + D( 2 )
! 103: * DELTA( 1 ) = -Z( 1 ) / TAU
! 104: * DELTA( 2 ) = Z( 2 ) / ( DEL-TAU )
! 105: ELSE
! 106: B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
! 107: C = RHO*Z( 2 )*Z( 2 )*DELSQ
! 108: *
! 109: * The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
! 110: *
! 111: IF( B.GT.ZERO ) THEN
! 112: TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) )
! 113: ELSE
! 114: TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO
! 115: END IF
! 116: *
! 117: * The following TAU is DSIGMA - D( 2 )
! 118: *
! 119: TAU = TAU / ( D( 2 )+SQRT( ABS( D( 2 )*D( 2 )+TAU ) ) )
! 120: DSIGMA = D( 2 ) + TAU
! 121: DELTA( 1 ) = -( DEL+TAU )
! 122: DELTA( 2 ) = -TAU
! 123: WORK( 1 ) = D( 1 ) + TAU + D( 2 )
! 124: WORK( 2 ) = TWO*D( 2 ) + TAU
! 125: * DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
! 126: * DELTA( 2 ) = -Z( 2 ) / TAU
! 127: END IF
! 128: * TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
! 129: * DELTA( 1 ) = DELTA( 1 ) / TEMP
! 130: * DELTA( 2 ) = DELTA( 2 ) / TEMP
! 131: ELSE
! 132: *
! 133: * Now I=2
! 134: *
! 135: B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
! 136: C = RHO*Z( 2 )*Z( 2 )*DELSQ
! 137: *
! 138: * The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
! 139: *
! 140: IF( B.GT.ZERO ) THEN
! 141: TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO
! 142: ELSE
! 143: TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) )
! 144: END IF
! 145: *
! 146: * The following TAU is DSIGMA - D( 2 )
! 147: *
! 148: TAU = TAU / ( D( 2 )+SQRT( D( 2 )*D( 2 )+TAU ) )
! 149: DSIGMA = D( 2 ) + TAU
! 150: DELTA( 1 ) = -( DEL+TAU )
! 151: DELTA( 2 ) = -TAU
! 152: WORK( 1 ) = D( 1 ) + TAU + D( 2 )
! 153: WORK( 2 ) = TWO*D( 2 ) + TAU
! 154: * DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
! 155: * DELTA( 2 ) = -Z( 2 ) / TAU
! 156: * TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
! 157: * DELTA( 1 ) = DELTA( 1 ) / TEMP
! 158: * DELTA( 2 ) = DELTA( 2 ) / TEMP
! 159: END IF
! 160: RETURN
! 161: *
! 162: * End of DLASD5
! 163: *
! 164: END
CVSweb interface <joel.bertrand@systella.fr>