LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ claic1()

subroutine claic1 ( integer job,
integer j,
complex, dimension( j ) x,
real sest,
complex, dimension( j ) w,
complex gamma,
real sestpr,
complex s,
complex c )

CLAIC1 applies one step of incremental condition estimation.

Download CLAIC1 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> CLAIC1 applies one step of incremental condition estimation in
!> its simplest version:
!>
!> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
!> lower triangular matrix L, such that
!>          twonorm(L*x) = sest
!> Then CLAIC1 computes sestpr, s, c such that
!> the vector
!>                 [ s*x ]
!>          xhat = [  c  ]
!> is an approximate singular vector of
!>                 [ L      0  ]
!>          Lhat = [ w**H gamma ]
!> in the sense that
!>          twonorm(Lhat*xhat) = sestpr.
!>
!> Depending on JOB, an estimate for the largest or smallest singular
!> value is computed.
!>
!> Note that [s c]**H and sestpr**2 is an eigenpair of the system
!>
!>     diag(sest*sest, 0) + [alpha  gamma] * [ conjg(alpha) ]
!>                                           [ conjg(gamma) ]
!>
!> where  alpha =  x**H*w.
!>
Parameters
[in]JOB
!>          JOB is INTEGER
!>          = 1: an estimate for the largest singular value is computed.
!>          = 2: an estimate for the smallest singular value is computed.
!>
[in]J
!>          J is INTEGER
!>          Length of X and W
!>
[in]X
!>          X is COMPLEX array, dimension (J)
!>          The j-vector x.
!>
[in]SEST
!>          SEST is REAL
!>          Estimated singular value of j by j matrix L
!>
[in]W
!>          W is COMPLEX array, dimension (J)
!>          The j-vector w.
!>
[in]GAMMA
!>          GAMMA is COMPLEX
!>          The diagonal element gamma.
!>
[out]SESTPR
!>          SESTPR is REAL
!>          Estimated singular value of (j+1) by (j+1) matrix Lhat.
!>
[out]S
!>          S is COMPLEX
!>          Sine needed in forming xhat.
!>
[out]C
!>          C is COMPLEX
!>          Cosine needed in forming xhat.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 132 of file claic1.f.

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 REAL SEST, SESTPR
141 COMPLEX C, GAMMA, S
142* ..
143* .. Array Arguments ..
144 COMPLEX W( J ), X( J )
145* ..
146*
147* =====================================================================
148*
149* .. Parameters ..
150 REAL ZERO, ONE, TWO
151 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
152 REAL HALF, FOUR
153 parameter( half = 0.5e0, four = 4.0e0 )
154* ..
155* .. Local Scalars ..
156 REAL ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
157 $ SCL, T, TEST, TMP, ZETA1, ZETA2
158 COMPLEX ALPHA, COSINE, SINE
159* ..
160* .. Intrinsic Functions ..
161 INTRINSIC abs, conjg, max, sqrt
162* ..
163* .. External Functions ..
164 REAL SLAMCH
165 COMPLEX CDOTC
166 EXTERNAL slamch, cdotc
167* ..
168* .. Executable Statements ..
169*
170 eps = slamch( 'Epsilon' )
171 alpha = cdotc( 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 = real( sqrt( s*conjg( s )+c*conjg( 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 = real( c / ( b+sqrt( b*b+c ) ) )
247 ELSE
248 t = real( sqrt( b*b+c ) - b )
249 END IF
250*
251 sine = -( alpha / absest ) / t
252 cosine = -( gamma / absest ) / ( one+t )
253 tmp = real( sqrt( sine * conjg( sine )
254 $ + cosine * conjg( 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 = -conjg( gamma )
274 cosine = conjg( alpha )
275 END IF
276 s1 = max( abs( sine ), abs( cosine ) )
277 s = sine / s1
278 c = cosine / s1
279 tmp = real( sqrt( s*conjg( s )+c*conjg( 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 scl = sqrt( one+tmp*tmp )
307 sestpr = absest*( tmp / scl )
308 s = -( conjg( gamma ) / s2 ) / scl
309 c = ( conjg( alpha ) / s2 ) / scl
310 ELSE
311 tmp = s2 / s1
312 scl = sqrt( one+tmp*tmp )
313 sestpr = absest / scl
314 s = -( conjg( gamma ) / s1 ) / scl
315 c = ( conjg( alpha ) / s1 ) / scl
316 END IF
317 RETURN
318 ELSE
319*
320* normal case
321*
322 zeta1 = absalp / absest
323 zeta2 = absgam / absest
324*
325 norma = max( one+zeta1*zeta1+zeta1*zeta2,
326 $ 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 = real( c / ( b+sqrt( abs( b*b-c ) ) ) )
338 sine = ( alpha / absest ) / ( one-t )
339 cosine = -( gamma / absest ) / 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 = real( -c / ( b+sqrt( b*b+c ) ) )
349 ELSE
350 t = real( b - sqrt( b*b+c ) )
351 END IF
352 sine = -( alpha / absest ) / t
353 cosine = -( gamma / absest ) / ( one+t )
354 sestpr = sqrt( one+t+four*eps*eps*norma )*absest
355 END IF
356 tmp = real( sqrt( sine * conjg( sine )
357 $ + cosine * conjg( cosine ) ) )
358 s = sine / tmp
359 c = cosine / tmp
360 RETURN
361*
362 END IF
363 END IF
364 RETURN
365*
366* End of CLAIC1
367*
cdotc
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
Here is the caller graph for this function: