1 | *> \brief \b SLAIC1 |
---|
2 | * |
---|
3 | * =========== DOCUMENTATION =========== |
---|
4 | * |
---|
5 | * Online html documentation available at |
---|
6 | * http://www.netlib.org/lapack/explore-html/ |
---|
7 | * |
---|
8 | *> \htmlonly |
---|
9 | *> Download SLAIC1 + dependencies |
---|
10 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaic1.f"> |
---|
11 | *> [TGZ]</a> |
---|
12 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaic1.f"> |
---|
13 | *> [ZIP]</a> |
---|
14 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaic1.f"> |
---|
15 | *> [TXT]</a> |
---|
16 | *> \endhtmlonly |
---|
17 | * |
---|
18 | * Definition: |
---|
19 | * =========== |
---|
20 | * |
---|
21 | * SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) |
---|
22 | * |
---|
23 | * .. Scalar Arguments .. |
---|
24 | * INTEGER J, JOB |
---|
25 | * REAL C, GAMMA, S, SEST, SESTPR |
---|
26 | * .. |
---|
27 | * .. Array Arguments .. |
---|
28 | * REAL W( J ), X( J ) |
---|
29 | * .. |
---|
30 | * |
---|
31 | * |
---|
32 | *> \par Purpose: |
---|
33 | * ============= |
---|
34 | *> |
---|
35 | *> \verbatim |
---|
36 | *> |
---|
37 | *> SLAIC1 applies one step of incremental condition estimation in |
---|
38 | *> its simplest version: |
---|
39 | *> |
---|
40 | *> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j |
---|
41 | *> lower triangular matrix L, such that |
---|
42 | *> twonorm(L*x) = sest |
---|
43 | *> Then SLAIC1 computes sestpr, s, c such that |
---|
44 | *> the vector |
---|
45 | *> [ s*x ] |
---|
46 | *> xhat = [ c ] |
---|
47 | *> is an approximate singular vector of |
---|
48 | *> [ L 0 ] |
---|
49 | *> Lhat = [ w**T gamma ] |
---|
50 | *> in the sense that |
---|
51 | *> twonorm(Lhat*xhat) = sestpr. |
---|
52 | *> |
---|
53 | *> Depending on JOB, an estimate for the largest or smallest singular |
---|
54 | *> value is computed. |
---|
55 | *> |
---|
56 | *> Note that [s c]**T and sestpr**2 is an eigenpair of the system |
---|
57 | *> |
---|
58 | *> diag(sest*sest, 0) + [alpha gamma] * [ alpha ] |
---|
59 | *> [ gamma ] |
---|
60 | *> |
---|
61 | *> where alpha = x**T*w. |
---|
62 | *> \endverbatim |
---|
63 | * |
---|
64 | * Arguments: |
---|
65 | * ========== |
---|
66 | * |
---|
67 | *> \param[in] JOB |
---|
68 | *> \verbatim |
---|
69 | *> JOB is INTEGER |
---|
70 | *> = 1: an estimate for the largest singular value is computed. |
---|
71 | *> = 2: an estimate for the smallest singular value is computed. |
---|
72 | *> \endverbatim |
---|
73 | *> |
---|
74 | *> \param[in] J |
---|
75 | *> \verbatim |
---|
76 | *> J is INTEGER |
---|
77 | *> Length of X and W |
---|
78 | *> \endverbatim |
---|
79 | *> |
---|
80 | *> \param[in] X |
---|
81 | *> \verbatim |
---|
82 | *> X is REAL array, dimension (J) |
---|
83 | *> The j-vector x. |
---|
84 | *> \endverbatim |
---|
85 | *> |
---|
86 | *> \param[in] SEST |
---|
87 | *> \verbatim |
---|
88 | *> SEST is REAL |
---|
89 | *> Estimated singular value of j by j matrix L |
---|
90 | *> \endverbatim |
---|
91 | *> |
---|
92 | *> \param[in] W |
---|
93 | *> \verbatim |
---|
94 | *> W is REAL array, dimension (J) |
---|
95 | *> The j-vector w. |
---|
96 | *> \endverbatim |
---|
97 | *> |
---|
98 | *> \param[in] GAMMA |
---|
99 | *> \verbatim |
---|
100 | *> GAMMA is REAL |
---|
101 | *> The diagonal element gamma. |
---|
102 | *> \endverbatim |
---|
103 | *> |
---|
104 | *> \param[out] SESTPR |
---|
105 | *> \verbatim |
---|
106 | *> SESTPR is REAL |
---|
107 | *> Estimated singular value of (j+1) by (j+1) matrix Lhat. |
---|
108 | *> \endverbatim |
---|
109 | *> |
---|
110 | *> \param[out] S |
---|
111 | *> \verbatim |
---|
112 | *> S is REAL |
---|
113 | *> Sine needed in forming xhat. |
---|
114 | *> \endverbatim |
---|
115 | *> |
---|
116 | *> \param[out] C |
---|
117 | *> \verbatim |
---|
118 | *> C is REAL |
---|
119 | *> Cosine needed in forming xhat. |
---|
120 | *> \endverbatim |
---|
121 | * |
---|
122 | * Authors: |
---|
123 | * ======== |
---|
124 | * |
---|
125 | *> \author Univ. of Tennessee |
---|
126 | *> \author Univ. of California Berkeley |
---|
127 | *> \author Univ. of Colorado Denver |
---|
128 | *> \author NAG Ltd. |
---|
129 | * |
---|
130 | *> \date November 2011 |
---|
131 | * |
---|
132 | *> \ingroup realOTHERauxiliary |
---|
133 | * |
---|
134 | * ===================================================================== |
---|
135 | SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) |
---|
136 | * |
---|
137 | * -- LAPACK auxiliary routine (version 3.4.0) -- |
---|
138 | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
---|
139 | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
---|
140 | * November 2011 |
---|
141 | * |
---|
142 | * .. Scalar Arguments .. |
---|
143 | INTEGER J, JOB |
---|
144 | REAL C, GAMMA, S, SEST, SESTPR |
---|
145 | * .. |
---|
146 | * .. Array Arguments .. |
---|
147 | REAL W( J ), X( J ) |
---|
148 | * .. |
---|
149 | * |
---|
150 | * ===================================================================== |
---|
151 | * |
---|
152 | * .. Parameters .. |
---|
153 | REAL ZERO, ONE, TWO |
---|
154 | PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) |
---|
155 | REAL HALF, FOUR |
---|
156 | PARAMETER ( HALF = 0.5E0, FOUR = 4.0E0 ) |
---|
157 | * .. |
---|
158 | * .. Local Scalars .. |
---|
159 | REAL ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS, |
---|
160 | $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2 |
---|
161 | * .. |
---|
162 | * .. Intrinsic Functions .. |
---|
163 | INTRINSIC ABS, MAX, SIGN, SQRT |
---|
164 | * .. |
---|
165 | * .. External Functions .. |
---|
166 | REAL SDOT, SLAMCH |
---|
167 | EXTERNAL SDOT, SLAMCH |
---|
168 | * .. |
---|
169 | * .. Executable Statements .. |
---|
170 | * |
---|
171 | EPS = SLAMCH( 'Epsilon' ) |
---|
172 | ALPHA = SDOT( J, X, 1, W, 1 ) |
---|
173 | * |
---|
174 | ABSALP = ABS( ALPHA ) |
---|
175 | ABSGAM = ABS( GAMMA ) |
---|
176 | ABSEST = ABS( SEST ) |
---|
177 | * |
---|
178 | IF( JOB.EQ.1 ) THEN |
---|
179 | * |
---|
180 | * Estimating largest singular value |
---|
181 | * |
---|
182 | * special cases |
---|
183 | * |
---|
184 | IF( SEST.EQ.ZERO ) THEN |
---|
185 | S1 = MAX( ABSGAM, ABSALP ) |
---|
186 | IF( S1.EQ.ZERO ) THEN |
---|
187 | S = ZERO |
---|
188 | C = ONE |
---|
189 | SESTPR = ZERO |
---|
190 | ELSE |
---|
191 | S = ALPHA / S1 |
---|
192 | C = GAMMA / S1 |
---|
193 | TMP = SQRT( S*S+C*C ) |
---|
194 | S = S / TMP |
---|
195 | C = C / TMP |
---|
196 | SESTPR = S1*TMP |
---|
197 | END IF |
---|
198 | RETURN |
---|
199 | ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN |
---|
200 | S = ONE |
---|
201 | C = ZERO |
---|
202 | TMP = MAX( ABSEST, ABSALP ) |
---|
203 | S1 = ABSEST / TMP |
---|
204 | S2 = ABSALP / TMP |
---|
205 | SESTPR = TMP*SQRT( S1*S1+S2*S2 ) |
---|
206 | RETURN |
---|
207 | ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN |
---|
208 | S1 = ABSGAM |
---|
209 | S2 = ABSEST |
---|
210 | IF( S1.LE.S2 ) THEN |
---|
211 | S = ONE |
---|
212 | C = ZERO |
---|
213 | SESTPR = S2 |
---|
214 | ELSE |
---|
215 | S = ZERO |
---|
216 | C = ONE |
---|
217 | SESTPR = S1 |
---|
218 | END IF |
---|
219 | RETURN |
---|
220 | ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN |
---|
221 | S1 = ABSGAM |
---|
222 | S2 = ABSALP |
---|
223 | IF( S1.LE.S2 ) THEN |
---|
224 | TMP = S1 / S2 |
---|
225 | S = SQRT( ONE+TMP*TMP ) |
---|
226 | SESTPR = S2*S |
---|
227 | C = ( GAMMA / S2 ) / S |
---|
228 | S = SIGN( ONE, ALPHA ) / S |
---|
229 | ELSE |
---|
230 | TMP = S2 / S1 |
---|
231 | C = SQRT( ONE+TMP*TMP ) |
---|
232 | SESTPR = S1*C |
---|
233 | S = ( ALPHA / S1 ) / C |
---|
234 | C = SIGN( ONE, GAMMA ) / C |
---|
235 | END IF |
---|
236 | RETURN |
---|
237 | ELSE |
---|
238 | * |
---|
239 | * normal case |
---|
240 | * |
---|
241 | ZETA1 = ALPHA / ABSEST |
---|
242 | ZETA2 = GAMMA / ABSEST |
---|
243 | * |
---|
244 | B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF |
---|
245 | C = ZETA1*ZETA1 |
---|
246 | IF( B.GT.ZERO ) THEN |
---|
247 | T = C / ( B+SQRT( B*B+C ) ) |
---|
248 | ELSE |
---|
249 | T = SQRT( B*B+C ) - B |
---|
250 | END IF |
---|
251 | * |
---|
252 | SINE = -ZETA1 / T |
---|
253 | COSINE = -ZETA2 / ( ONE+T ) |
---|
254 | TMP = SQRT( SINE*SINE+COSINE*COSINE ) |
---|
255 | S = SINE / TMP |
---|
256 | C = COSINE / TMP |
---|
257 | SESTPR = SQRT( T+ONE )*ABSEST |
---|
258 | RETURN |
---|
259 | END IF |
---|
260 | * |
---|
261 | ELSE IF( JOB.EQ.2 ) THEN |
---|
262 | * |
---|
263 | * Estimating smallest singular value |
---|
264 | * |
---|
265 | * special cases |
---|
266 | * |
---|
267 | IF( SEST.EQ.ZERO ) THEN |
---|
268 | SESTPR = ZERO |
---|
269 | IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN |
---|
270 | SINE = ONE |
---|
271 | COSINE = ZERO |
---|
272 | ELSE |
---|
273 | SINE = -GAMMA |
---|
274 | COSINE = ALPHA |
---|
275 | END IF |
---|
276 | S1 = MAX( ABS( SINE ), ABS( COSINE ) ) |
---|
277 | S = SINE / S1 |
---|
278 | C = COSINE / S1 |
---|
279 | TMP = SQRT( S*S+C*C ) |
---|
280 | S = S / TMP |
---|
281 | C = C / TMP |
---|
282 | RETURN |
---|
283 | ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN |
---|
284 | S = ZERO |
---|
285 | C = ONE |
---|
286 | SESTPR = ABSGAM |
---|
287 | RETURN |
---|
288 | ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN |
---|
289 | S1 = ABSGAM |
---|
290 | S2 = ABSEST |
---|
291 | IF( S1.LE.S2 ) THEN |
---|
292 | S = ZERO |
---|
293 | C = ONE |
---|
294 | SESTPR = S1 |
---|
295 | ELSE |
---|
296 | S = ONE |
---|
297 | C = ZERO |
---|
298 | SESTPR = S2 |
---|
299 | END IF |
---|
300 | RETURN |
---|
301 | ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN |
---|
302 | S1 = ABSGAM |
---|
303 | S2 = ABSALP |
---|
304 | IF( S1.LE.S2 ) THEN |
---|
305 | TMP = S1 / S2 |
---|
306 | C = SQRT( ONE+TMP*TMP ) |
---|
307 | SESTPR = ABSEST*( TMP / C ) |
---|
308 | S = -( GAMMA / S2 ) / C |
---|
309 | C = SIGN( ONE, ALPHA ) / C |
---|
310 | ELSE |
---|
311 | TMP = S2 / S1 |
---|
312 | S = SQRT( ONE+TMP*TMP ) |
---|
313 | SESTPR = ABSEST / S |
---|
314 | C = ( ALPHA / S1 ) / S |
---|
315 | S = -SIGN( ONE, GAMMA ) / S |
---|
316 | END IF |
---|
317 | RETURN |
---|
318 | ELSE |
---|
319 | * |
---|
320 | * normal case |
---|
321 | * |
---|
322 | ZETA1 = ALPHA / ABSEST |
---|
323 | ZETA2 = GAMMA / ABSEST |
---|
324 | * |
---|
325 | NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ), |
---|
326 | $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 ) |
---|
327 | * |
---|
328 | * See if root is closer to zero or to ONE |
---|
329 | * |
---|
330 | TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 ) |
---|
331 | IF( TEST.GE.ZERO ) THEN |
---|
332 | * |
---|
333 | * root is close to zero, compute directly |
---|
334 | * |
---|
335 | B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF |
---|
336 | C = ZETA2*ZETA2 |
---|
337 | T = C / ( B+SQRT( ABS( B*B-C ) ) ) |
---|
338 | SINE = ZETA1 / ( ONE-T ) |
---|
339 | COSINE = -ZETA2 / T |
---|
340 | SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST |
---|
341 | ELSE |
---|
342 | * |
---|
343 | * root is closer to ONE, shift by that amount |
---|
344 | * |
---|
345 | B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF |
---|
346 | C = ZETA1*ZETA1 |
---|
347 | IF( B.GE.ZERO ) THEN |
---|
348 | T = -C / ( B+SQRT( B*B+C ) ) |
---|
349 | ELSE |
---|
350 | T = B - SQRT( B*B+C ) |
---|
351 | END IF |
---|
352 | SINE = -ZETA1 / T |
---|
353 | COSINE = -ZETA2 / ( ONE+T ) |
---|
354 | SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST |
---|
355 | END IF |
---|
356 | TMP = SQRT( SINE*SINE+COSINE*COSINE ) |
---|
357 | S = SINE / TMP |
---|
358 | C = COSINE / TMP |
---|
359 | RETURN |
---|
360 | * |
---|
361 | END IF |
---|
362 | END IF |
---|
363 | RETURN |
---|
364 | * |
---|
365 | * End of SLAIC1 |
---|
366 | * |
---|
367 | END |
---|