1 SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER J, JOB
10 DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION W( J ), X( J )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DLAIC1 applies one step of incremental condition estimation in
20 * its simplest version:
21 *
22 * Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
23 * lower triangular matrix L, such that
24 * twonorm(L*x) = sest
25 * Then DLAIC1 computes sestpr, s, c such that
26 * the vector
27 * [ s*x ]
28 * xhat = [ c ]
29 * is an approximate singular vector of
30 * [ L 0 ]
31 * Lhat = [ w**T gamma ]
32 * in the sense that
33 * twonorm(Lhat*xhat) = sestpr.
34 *
35 * Depending on JOB, an estimate for the largest or smallest singular
36 * value is computed.
37 *
38 * Note that [s c]**T and sestpr**2 is an eigenpair of the system
39 *
40 * diag(sest*sest, 0) + [alpha gamma] * [ alpha ]
41 * [ gamma ]
42 *
43 * where alpha = x**T*w.
44 *
45 * Arguments
46 * =========
47 *
48 * JOB (input) INTEGER
49 * = 1: an estimate for the largest singular value is computed.
50 * = 2: an estimate for the smallest singular value is computed.
51 *
52 * J (input) INTEGER
53 * Length of X and W
54 *
55 * X (input) DOUBLE PRECISION array, dimension (J)
56 * The j-vector x.
57 *
58 * SEST (input) DOUBLE PRECISION
59 * Estimated singular value of j by j matrix L
60 *
61 * W (input) DOUBLE PRECISION array, dimension (J)
62 * The j-vector w.
63 *
64 * GAMMA (input) DOUBLE PRECISION
65 * The diagonal element gamma.
66 *
67 * SESTPR (output) DOUBLE PRECISION
68 * Estimated singular value of (j+1) by (j+1) matrix Lhat.
69 *
70 * S (output) DOUBLE PRECISION
71 * Sine needed in forming xhat.
72 *
73 * C (output) DOUBLE PRECISION
74 * Cosine needed in forming xhat.
75 *
76 * =====================================================================
77 *
78 * .. Parameters ..
79 DOUBLE PRECISION ZERO, ONE, TWO
80 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
81 DOUBLE PRECISION HALF, FOUR
82 PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 )
83 * ..
84 * .. Local Scalars ..
85 DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
86 $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC ABS, MAX, SIGN, SQRT
90 * ..
91 * .. External Functions ..
92 DOUBLE PRECISION DDOT, DLAMCH
93 EXTERNAL DDOT, DLAMCH
94 * ..
95 * .. Executable Statements ..
96 *
97 EPS = DLAMCH( 'Epsilon' )
98 ALPHA = DDOT( J, X, 1, W, 1 )
99 *
100 ABSALP = ABS( ALPHA )
101 ABSGAM = ABS( GAMMA )
102 ABSEST = ABS( SEST )
103 *
104 IF( JOB.EQ.1 ) THEN
105 *
106 * Estimating largest singular value
107 *
108 * special cases
109 *
110 IF( SEST.EQ.ZERO ) THEN
111 S1 = MAX( ABSGAM, ABSALP )
112 IF( S1.EQ.ZERO ) THEN
113 S = ZERO
114 C = ONE
115 SESTPR = ZERO
116 ELSE
117 S = ALPHA / S1
118 C = GAMMA / S1
119 TMP = SQRT( S*S+C*C )
120 S = S / TMP
121 C = C / TMP
122 SESTPR = S1*TMP
123 END IF
124 RETURN
125 ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
126 S = ONE
127 C = ZERO
128 TMP = MAX( ABSEST, ABSALP )
129 S1 = ABSEST / TMP
130 S2 = ABSALP / TMP
131 SESTPR = TMP*SQRT( S1*S1+S2*S2 )
132 RETURN
133 ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
134 S1 = ABSGAM
135 S2 = ABSEST
136 IF( S1.LE.S2 ) THEN
137 S = ONE
138 C = ZERO
139 SESTPR = S2
140 ELSE
141 S = ZERO
142 C = ONE
143 SESTPR = S1
144 END IF
145 RETURN
146 ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
147 S1 = ABSGAM
148 S2 = ABSALP
149 IF( S1.LE.S2 ) THEN
150 TMP = S1 / S2
151 S = SQRT( ONE+TMP*TMP )
152 SESTPR = S2*S
153 C = ( GAMMA / S2 ) / S
154 S = SIGN( ONE, ALPHA ) / S
155 ELSE
156 TMP = S2 / S1
157 C = SQRT( ONE+TMP*TMP )
158 SESTPR = S1*C
159 S = ( ALPHA / S1 ) / C
160 C = SIGN( ONE, GAMMA ) / C
161 END IF
162 RETURN
163 ELSE
164 *
165 * normal case
166 *
167 ZETA1 = ALPHA / ABSEST
168 ZETA2 = GAMMA / ABSEST
169 *
170 B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
171 C = ZETA1*ZETA1
172 IF( B.GT.ZERO ) THEN
173 T = C / ( B+SQRT( B*B+C ) )
174 ELSE
175 T = SQRT( B*B+C ) - B
176 END IF
177 *
178 SINE = -ZETA1 / T
179 COSINE = -ZETA2 / ( ONE+T )
180 TMP = SQRT( SINE*SINE+COSINE*COSINE )
181 S = SINE / TMP
182 C = COSINE / TMP
183 SESTPR = SQRT( T+ONE )*ABSEST
184 RETURN
185 END IF
186 *
187 ELSE IF( JOB.EQ.2 ) THEN
188 *
189 * Estimating smallest singular value
190 *
191 * special cases
192 *
193 IF( SEST.EQ.ZERO ) THEN
194 SESTPR = ZERO
195 IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
196 SINE = ONE
197 COSINE = ZERO
198 ELSE
199 SINE = -GAMMA
200 COSINE = ALPHA
201 END IF
202 S1 = MAX( ABS( SINE ), ABS( COSINE ) )
203 S = SINE / S1
204 C = COSINE / S1
205 TMP = SQRT( S*S+C*C )
206 S = S / TMP
207 C = C / TMP
208 RETURN
209 ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
210 S = ZERO
211 C = ONE
212 SESTPR = ABSGAM
213 RETURN
214 ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
215 S1 = ABSGAM
216 S2 = ABSEST
217 IF( S1.LE.S2 ) THEN
218 S = ZERO
219 C = ONE
220 SESTPR = S1
221 ELSE
222 S = ONE
223 C = ZERO
224 SESTPR = S2
225 END IF
226 RETURN
227 ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
228 S1 = ABSGAM
229 S2 = ABSALP
230 IF( S1.LE.S2 ) THEN
231 TMP = S1 / S2
232 C = SQRT( ONE+TMP*TMP )
233 SESTPR = ABSEST*( TMP / C )
234 S = -( GAMMA / S2 ) / C
235 C = SIGN( ONE, ALPHA ) / C
236 ELSE
237 TMP = S2 / S1
238 S = SQRT( ONE+TMP*TMP )
239 SESTPR = ABSEST / S
240 C = ( ALPHA / S1 ) / S
241 S = -SIGN( ONE, GAMMA ) / S
242 END IF
243 RETURN
244 ELSE
245 *
246 * normal case
247 *
248 ZETA1 = ALPHA / ABSEST
249 ZETA2 = GAMMA / ABSEST
250 *
251 NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
252 $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
253 *
254 * See if root is closer to zero or to ONE
255 *
256 TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
257 IF( TEST.GE.ZERO ) THEN
258 *
259 * root is close to zero, compute directly
260 *
261 B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
262 C = ZETA2*ZETA2
263 T = C / ( B+SQRT( ABS( B*B-C ) ) )
264 SINE = ZETA1 / ( ONE-T )
265 COSINE = -ZETA2 / T
266 SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
267 ELSE
268 *
269 * root is closer to ONE, shift by that amount
270 *
271 B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
272 C = ZETA1*ZETA1
273 IF( B.GE.ZERO ) THEN
274 T = -C / ( B+SQRT( B*B+C ) )
275 ELSE
276 T = B - SQRT( B*B+C )
277 END IF
278 SINE = -ZETA1 / T
279 COSINE = -ZETA2 / ( ONE+T )
280 SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
281 END IF
282 TMP = SQRT( SINE*SINE+COSINE*COSINE )
283 S = SINE / TMP
284 C = COSINE / TMP
285 RETURN
286 *
287 END IF
288 END IF
289 RETURN
290 *
291 * End of DLAIC1
292 *
293 END
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER J, JOB
10 DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION W( J ), X( J )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DLAIC1 applies one step of incremental condition estimation in
20 * its simplest version:
21 *
22 * Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
23 * lower triangular matrix L, such that
24 * twonorm(L*x) = sest
25 * Then DLAIC1 computes sestpr, s, c such that
26 * the vector
27 * [ s*x ]
28 * xhat = [ c ]
29 * is an approximate singular vector of
30 * [ L 0 ]
31 * Lhat = [ w**T gamma ]
32 * in the sense that
33 * twonorm(Lhat*xhat) = sestpr.
34 *
35 * Depending on JOB, an estimate for the largest or smallest singular
36 * value is computed.
37 *
38 * Note that [s c]**T and sestpr**2 is an eigenpair of the system
39 *
40 * diag(sest*sest, 0) + [alpha gamma] * [ alpha ]
41 * [ gamma ]
42 *
43 * where alpha = x**T*w.
44 *
45 * Arguments
46 * =========
47 *
48 * JOB (input) INTEGER
49 * = 1: an estimate for the largest singular value is computed.
50 * = 2: an estimate for the smallest singular value is computed.
51 *
52 * J (input) INTEGER
53 * Length of X and W
54 *
55 * X (input) DOUBLE PRECISION array, dimension (J)
56 * The j-vector x.
57 *
58 * SEST (input) DOUBLE PRECISION
59 * Estimated singular value of j by j matrix L
60 *
61 * W (input) DOUBLE PRECISION array, dimension (J)
62 * The j-vector w.
63 *
64 * GAMMA (input) DOUBLE PRECISION
65 * The diagonal element gamma.
66 *
67 * SESTPR (output) DOUBLE PRECISION
68 * Estimated singular value of (j+1) by (j+1) matrix Lhat.
69 *
70 * S (output) DOUBLE PRECISION
71 * Sine needed in forming xhat.
72 *
73 * C (output) DOUBLE PRECISION
74 * Cosine needed in forming xhat.
75 *
76 * =====================================================================
77 *
78 * .. Parameters ..
79 DOUBLE PRECISION ZERO, ONE, TWO
80 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
81 DOUBLE PRECISION HALF, FOUR
82 PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 )
83 * ..
84 * .. Local Scalars ..
85 DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
86 $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC ABS, MAX, SIGN, SQRT
90 * ..
91 * .. External Functions ..
92 DOUBLE PRECISION DDOT, DLAMCH
93 EXTERNAL DDOT, DLAMCH
94 * ..
95 * .. Executable Statements ..
96 *
97 EPS = DLAMCH( 'Epsilon' )
98 ALPHA = DDOT( J, X, 1, W, 1 )
99 *
100 ABSALP = ABS( ALPHA )
101 ABSGAM = ABS( GAMMA )
102 ABSEST = ABS( SEST )
103 *
104 IF( JOB.EQ.1 ) THEN
105 *
106 * Estimating largest singular value
107 *
108 * special cases
109 *
110 IF( SEST.EQ.ZERO ) THEN
111 S1 = MAX( ABSGAM, ABSALP )
112 IF( S1.EQ.ZERO ) THEN
113 S = ZERO
114 C = ONE
115 SESTPR = ZERO
116 ELSE
117 S = ALPHA / S1
118 C = GAMMA / S1
119 TMP = SQRT( S*S+C*C )
120 S = S / TMP
121 C = C / TMP
122 SESTPR = S1*TMP
123 END IF
124 RETURN
125 ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
126 S = ONE
127 C = ZERO
128 TMP = MAX( ABSEST, ABSALP )
129 S1 = ABSEST / TMP
130 S2 = ABSALP / TMP
131 SESTPR = TMP*SQRT( S1*S1+S2*S2 )
132 RETURN
133 ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
134 S1 = ABSGAM
135 S2 = ABSEST
136 IF( S1.LE.S2 ) THEN
137 S = ONE
138 C = ZERO
139 SESTPR = S2
140 ELSE
141 S = ZERO
142 C = ONE
143 SESTPR = S1
144 END IF
145 RETURN
146 ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
147 S1 = ABSGAM
148 S2 = ABSALP
149 IF( S1.LE.S2 ) THEN
150 TMP = S1 / S2
151 S = SQRT( ONE+TMP*TMP )
152 SESTPR = S2*S
153 C = ( GAMMA / S2 ) / S
154 S = SIGN( ONE, ALPHA ) / S
155 ELSE
156 TMP = S2 / S1
157 C = SQRT( ONE+TMP*TMP )
158 SESTPR = S1*C
159 S = ( ALPHA / S1 ) / C
160 C = SIGN( ONE, GAMMA ) / C
161 END IF
162 RETURN
163 ELSE
164 *
165 * normal case
166 *
167 ZETA1 = ALPHA / ABSEST
168 ZETA2 = GAMMA / ABSEST
169 *
170 B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
171 C = ZETA1*ZETA1
172 IF( B.GT.ZERO ) THEN
173 T = C / ( B+SQRT( B*B+C ) )
174 ELSE
175 T = SQRT( B*B+C ) - B
176 END IF
177 *
178 SINE = -ZETA1 / T
179 COSINE = -ZETA2 / ( ONE+T )
180 TMP = SQRT( SINE*SINE+COSINE*COSINE )
181 S = SINE / TMP
182 C = COSINE / TMP
183 SESTPR = SQRT( T+ONE )*ABSEST
184 RETURN
185 END IF
186 *
187 ELSE IF( JOB.EQ.2 ) THEN
188 *
189 * Estimating smallest singular value
190 *
191 * special cases
192 *
193 IF( SEST.EQ.ZERO ) THEN
194 SESTPR = ZERO
195 IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
196 SINE = ONE
197 COSINE = ZERO
198 ELSE
199 SINE = -GAMMA
200 COSINE = ALPHA
201 END IF
202 S1 = MAX( ABS( SINE ), ABS( COSINE ) )
203 S = SINE / S1
204 C = COSINE / S1
205 TMP = SQRT( S*S+C*C )
206 S = S / TMP
207 C = C / TMP
208 RETURN
209 ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
210 S = ZERO
211 C = ONE
212 SESTPR = ABSGAM
213 RETURN
214 ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
215 S1 = ABSGAM
216 S2 = ABSEST
217 IF( S1.LE.S2 ) THEN
218 S = ZERO
219 C = ONE
220 SESTPR = S1
221 ELSE
222 S = ONE
223 C = ZERO
224 SESTPR = S2
225 END IF
226 RETURN
227 ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
228 S1 = ABSGAM
229 S2 = ABSALP
230 IF( S1.LE.S2 ) THEN
231 TMP = S1 / S2
232 C = SQRT( ONE+TMP*TMP )
233 SESTPR = ABSEST*( TMP / C )
234 S = -( GAMMA / S2 ) / C
235 C = SIGN( ONE, ALPHA ) / C
236 ELSE
237 TMP = S2 / S1
238 S = SQRT( ONE+TMP*TMP )
239 SESTPR = ABSEST / S
240 C = ( ALPHA / S1 ) / S
241 S = -SIGN( ONE, GAMMA ) / S
242 END IF
243 RETURN
244 ELSE
245 *
246 * normal case
247 *
248 ZETA1 = ALPHA / ABSEST
249 ZETA2 = GAMMA / ABSEST
250 *
251 NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
252 $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
253 *
254 * See if root is closer to zero or to ONE
255 *
256 TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
257 IF( TEST.GE.ZERO ) THEN
258 *
259 * root is close to zero, compute directly
260 *
261 B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
262 C = ZETA2*ZETA2
263 T = C / ( B+SQRT( ABS( B*B-C ) ) )
264 SINE = ZETA1 / ( ONE-T )
265 COSINE = -ZETA2 / T
266 SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
267 ELSE
268 *
269 * root is closer to ONE, shift by that amount
270 *
271 B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
272 C = ZETA1*ZETA1
273 IF( B.GE.ZERO ) THEN
274 T = -C / ( B+SQRT( B*B+C ) )
275 ELSE
276 T = B - SQRT( B*B+C )
277 END IF
278 SINE = -ZETA1 / T
279 COSINE = -ZETA2 / ( ONE+T )
280 SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
281 END IF
282 TMP = SQRT( SINE*SINE+COSINE*COSINE )
283 S = SINE / TMP
284 C = COSINE / TMP
285 RETURN
286 *
287 END IF
288 END IF
289 RETURN
290 *
291 * End of DLAIC1
292 *
293 END