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
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