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