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