Annotation of rpl/lapack/lapack/dlasd5.f, revision 1.3
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>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd5.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack