LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaic1.f
Go to the documentation of this file.
1*> \brief \b DLAIC1 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
9*> Download DLAIC1 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaic1.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaic1.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaic1.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
22*
23* .. Scalar Arguments ..
24* INTEGER J, JOB
25* DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION W( J ), X( J )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DLAIC1 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 DLAIC1 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 DOUBLE PRECISION array, dimension (J)
83*> The j-vector x.
84*> \endverbatim
85*>
86*> \param[in] SEST
87*> \verbatim
88*> SEST is DOUBLE PRECISION
89*> Estimated singular value of j by j matrix L
90*> \endverbatim
91*>
92*> \param[in] W
93*> \verbatim
94*> W is DOUBLE PRECISION array, dimension (J)
95*> The j-vector w.
96*> \endverbatim
97*>
98*> \param[in] GAMMA
99*> \verbatim
100*> GAMMA is DOUBLE PRECISION
101*> The diagonal element gamma.
102*> \endverbatim
103*>
104*> \param[out] SESTPR
105*> \verbatim
106*> SESTPR is DOUBLE PRECISION
107*> Estimated singular value of (j+1) by (j+1) matrix Lhat.
108*> \endverbatim
109*>
110*> \param[out] S
111*> \verbatim
112*> S is DOUBLE PRECISION
113*> Sine needed in forming xhat.
114*> \endverbatim
115*>
116*> \param[out] C
117*> \verbatim
118*> C is DOUBLE PRECISION
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*> \ingroup laic1
131*
132* =====================================================================
133 SUBROUTINE dlaic1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
134*
135* -- LAPACK auxiliary routine --
136* -- LAPACK is a software package provided by Univ. of Tennessee, --
137* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*
139* .. Scalar Arguments ..
140 INTEGER J, JOB
141 DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR
142* ..
143* .. Array Arguments ..
144 DOUBLE PRECISION 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, ALPHA, B, COSINE, EPS,
157 $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
158* ..
159* .. Intrinsic Functions ..
160 INTRINSIC abs, max, sign, sqrt
161* ..
162* .. External Functions ..
163 DOUBLE PRECISION DDOT, DLAMCH
164 EXTERNAL ddot, dlamch
165* ..
166* .. Executable Statements ..
167*
168 eps = dlamch( 'Epsilon' )
169 alpha = ddot( j, x, 1, w, 1 )
170*
171 absalp = abs( alpha )
172 absgam = abs( gamma )
173 absest = abs( sest )
174*
175 IF( job.EQ.1 ) THEN
176*
177* Estimating largest singular value
178*
179* special cases
180*
181 IF( sest.EQ.zero ) THEN
182 s1 = max( absgam, absalp )
183 IF( s1.EQ.zero ) THEN
184 s = zero
185 c = one
186 sestpr = zero
187 ELSE
188 s = alpha / s1
189 c = gamma / s1
190 tmp = sqrt( s*s+c*c )
191 s = s / tmp
192 c = c / tmp
193 sestpr = s1*tmp
194 END IF
195 RETURN
196 ELSE IF( absgam.LE.eps*absest ) THEN
197 s = one
198 c = zero
199 tmp = max( absest, absalp )
200 s1 = absest / tmp
201 s2 = absalp / tmp
202 sestpr = tmp*sqrt( s1*s1+s2*s2 )
203 RETURN
204 ELSE IF( absalp.LE.eps*absest ) THEN
205 s1 = absgam
206 s2 = absest
207 IF( s1.LE.s2 ) THEN
208 s = one
209 c = zero
210 sestpr = s2
211 ELSE
212 s = zero
213 c = one
214 sestpr = s1
215 END IF
216 RETURN
217 ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam ) THEN
218 s1 = absgam
219 s2 = absalp
220 IF( s1.LE.s2 ) THEN
221 tmp = s1 / s2
222 s = sqrt( one+tmp*tmp )
223 sestpr = s2*s
224 c = ( gamma / s2 ) / s
225 s = sign( one, alpha ) / s
226 ELSE
227 tmp = s2 / s1
228 c = sqrt( one+tmp*tmp )
229 sestpr = s1*c
230 s = ( alpha / s1 ) / c
231 c = sign( one, gamma ) / c
232 END IF
233 RETURN
234 ELSE
235*
236* normal case
237*
238 zeta1 = alpha / absest
239 zeta2 = gamma / absest
240*
241 b = ( one-zeta1*zeta1-zeta2*zeta2 )*half
242 c = zeta1*zeta1
243 IF( b.GT.zero ) THEN
244 t = c / ( b+sqrt( b*b+c ) )
245 ELSE
246 t = sqrt( b*b+c ) - b
247 END IF
248*
249 sine = -zeta1 / t
250 cosine = -zeta2 / ( one+t )
251 tmp = sqrt( sine*sine+cosine*cosine )
252 s = sine / tmp
253 c = cosine / tmp
254 sestpr = sqrt( t+one )*absest
255 RETURN
256 END IF
257*
258 ELSE IF( job.EQ.2 ) THEN
259*
260* Estimating smallest singular value
261*
262* special cases
263*
264 IF( sest.EQ.zero ) THEN
265 sestpr = zero
266 IF( max( absgam, absalp ).EQ.zero ) THEN
267 sine = one
268 cosine = zero
269 ELSE
270 sine = -gamma
271 cosine = alpha
272 END IF
273 s1 = max( abs( sine ), abs( cosine ) )
274 s = sine / s1
275 c = cosine / s1
276 tmp = sqrt( s*s+c*c )
277 s = s / tmp
278 c = c / tmp
279 RETURN
280 ELSE IF( absgam.LE.eps*absest ) THEN
281 s = zero
282 c = one
283 sestpr = absgam
284 RETURN
285 ELSE IF( absalp.LE.eps*absest ) THEN
286 s1 = absgam
287 s2 = absest
288 IF( s1.LE.s2 ) THEN
289 s = zero
290 c = one
291 sestpr = s1
292 ELSE
293 s = one
294 c = zero
295 sestpr = s2
296 END IF
297 RETURN
298 ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam ) THEN
299 s1 = absgam
300 s2 = absalp
301 IF( s1.LE.s2 ) THEN
302 tmp = s1 / s2
303 c = sqrt( one+tmp*tmp )
304 sestpr = absest*( tmp / c )
305 s = -( gamma / s2 ) / c
306 c = sign( one, alpha ) / c
307 ELSE
308 tmp = s2 / s1
309 s = sqrt( one+tmp*tmp )
310 sestpr = absest / s
311 c = ( alpha / s1 ) / s
312 s = -sign( one, gamma ) / s
313 END IF
314 RETURN
315 ELSE
316*
317* normal case
318*
319 zeta1 = alpha / absest
320 zeta2 = gamma / absest
321*
322 norma = max( one+zeta1*zeta1+abs( zeta1*zeta2 ),
323 $ abs( zeta1*zeta2 )+zeta2*zeta2 )
324*
325* See if root is closer to zero or to ONE
326*
327 test = one + two*( zeta1-zeta2 )*( zeta1+zeta2 )
328 IF( test.GE.zero ) THEN
329*
330* root is close to zero, compute directly
331*
332 b = ( zeta1*zeta1+zeta2*zeta2+one )*half
333 c = zeta2*zeta2
334 t = c / ( b+sqrt( abs( b*b-c ) ) )
335 sine = zeta1 / ( one-t )
336 cosine = -zeta2 / t
337 sestpr = sqrt( t+four*eps*eps*norma )*absest
338 ELSE
339*
340* root is closer to ONE, shift by that amount
341*
342 b = ( zeta2*zeta2+zeta1*zeta1-one )*half
343 c = zeta1*zeta1
344 IF( b.GE.zero ) THEN
345 t = -c / ( b+sqrt( b*b+c ) )
346 ELSE
347 t = b - sqrt( b*b+c )
348 END IF
349 sine = -zeta1 / t
350 cosine = -zeta2 / ( one+t )
351 sestpr = sqrt( one+t+four*eps*eps*norma )*absest
352 END IF
353 tmp = sqrt( sine*sine+cosine*cosine )
354 s = sine / tmp
355 c = cosine / tmp
356 RETURN
357*
358 END IF
359 END IF
360 RETURN
361*
362* End of DLAIC1
363*
364 END
subroutine dlaic1(job, j, x, sest, w, gamma, sestpr, s, c)
DLAIC1 applies one step of incremental condition estimation.
Definition dlaic1.f:134