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

◆ slaic1()

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

SLAIC1 applies one step of incremental condition estimation.

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

Purpose:
!>
!> SLAIC1 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 SLAIC1 computes sestpr, s, c such that
!> the vector
!>                 [ s*x ]
!>          xhat = [  c  ]
!> is an approximate singular vector of
!>                 [ L      0  ]
!>          Lhat = [ w**T 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]**T and sestpr**2 is an eigenpair of the system
!>
!>     diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
!>                                           [ gamma ]
!>
!> where  alpha =  x**T*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 REAL 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 REAL array, dimension (J)
!>          The j-vector w.
!>
[in]GAMMA
!>          GAMMA is REAL
!>          The diagonal element gamma.
!>
[out]SESTPR
!>          SESTPR is REAL
!>          Estimated singular value of (j+1) by (j+1) matrix Lhat.
!>
[out]S
!>          S is REAL
!>          Sine needed in forming xhat.
!>
[out]C
!>          C is REAL
!>          Cosine needed in forming xhat.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 131 of file slaic1.f.

132*
133* -- LAPACK auxiliary routine --
134* -- LAPACK is a software package provided by Univ. of Tennessee, --
135* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136*
137* .. Scalar Arguments ..
138 INTEGER J, JOB
139 REAL C, GAMMA, S, SEST, SESTPR
140* ..
141* .. Array Arguments ..
142 REAL W( J ), X( J )
143* ..
144*
145* =====================================================================
146*
147* .. Parameters ..
148 REAL ZERO, ONE, TWO
149 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
150 REAL HALF, FOUR
151 parameter( half = 0.5e0, four = 4.0e0 )
152* ..
153* .. Local Scalars ..
154 REAL ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
155 $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
156* ..
157* .. Intrinsic Functions ..
158 INTRINSIC abs, max, sign, sqrt
159* ..
160* .. External Functions ..
161 REAL SDOT, SLAMCH
162 EXTERNAL sdot, slamch
163* ..
164* .. Executable Statements ..
165*
166 eps = slamch( 'Epsilon' )
167 alpha = sdot( j, x, 1, w, 1 )
168*
169 absalp = abs( alpha )
170 absgam = abs( gamma )
171 absest = abs( sest )
172*
173 IF( job.EQ.1 ) THEN
174*
175* Estimating largest singular value
176*
177* special cases
178*
179 IF( sest.EQ.zero ) THEN
180 s1 = max( absgam, absalp )
181 IF( s1.EQ.zero ) THEN
182 s = zero
183 c = one
184 sestpr = zero
185 ELSE
186 s = alpha / s1
187 c = gamma / s1
188 tmp = sqrt( s*s+c*c )
189 s = s / tmp
190 c = c / tmp
191 sestpr = s1*tmp
192 END IF
193 RETURN
194 ELSE IF( absgam.LE.eps*absest ) THEN
195 s = one
196 c = zero
197 tmp = max( absest, absalp )
198 s1 = absest / tmp
199 s2 = absalp / tmp
200 sestpr = tmp*sqrt( s1*s1+s2*s2 )
201 RETURN
202 ELSE IF( absalp.LE.eps*absest ) THEN
203 s1 = absgam
204 s2 = absest
205 IF( s1.LE.s2 ) THEN
206 s = one
207 c = zero
208 sestpr = s2
209 ELSE
210 s = zero
211 c = one
212 sestpr = s1
213 END IF
214 RETURN
215 ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam ) THEN
216 s1 = absgam
217 s2 = absalp
218 IF( s1.LE.s2 ) THEN
219 tmp = s1 / s2
220 s = sqrt( one+tmp*tmp )
221 sestpr = s2*s
222 c = ( gamma / s2 ) / s
223 s = sign( one, alpha ) / s
224 ELSE
225 tmp = s2 / s1
226 c = sqrt( one+tmp*tmp )
227 sestpr = s1*c
228 s = ( alpha / s1 ) / c
229 c = sign( one, gamma ) / c
230 END IF
231 RETURN
232 ELSE
233*
234* normal case
235*
236 zeta1 = alpha / absest
237 zeta2 = gamma / absest
238*
239 b = ( one-zeta1*zeta1-zeta2*zeta2 )*half
240 c = zeta1*zeta1
241 IF( b.GT.zero ) THEN
242 t = c / ( b+sqrt( b*b+c ) )
243 ELSE
244 t = sqrt( b*b+c ) - b
245 END IF
246*
247 sine = -zeta1 / t
248 cosine = -zeta2 / ( one+t )
249 tmp = sqrt( sine*sine+cosine*cosine )
250 s = sine / tmp
251 c = cosine / tmp
252 sestpr = sqrt( t+one )*absest
253 RETURN
254 END IF
255*
256 ELSE IF( job.EQ.2 ) THEN
257*
258* Estimating smallest singular value
259*
260* special cases
261*
262 IF( sest.EQ.zero ) THEN
263 sestpr = zero
264 IF( max( absgam, absalp ).EQ.zero ) THEN
265 sine = one
266 cosine = zero
267 ELSE
268 sine = -gamma
269 cosine = alpha
270 END IF
271 s1 = max( abs( sine ), abs( cosine ) )
272 s = sine / s1
273 c = cosine / s1
274 tmp = sqrt( s*s+c*c )
275 s = s / tmp
276 c = c / tmp
277 RETURN
278 ELSE IF( absgam.LE.eps*absest ) THEN
279 s = zero
280 c = one
281 sestpr = absgam
282 RETURN
283 ELSE IF( absalp.LE.eps*absest ) THEN
284 s1 = absgam
285 s2 = absest
286 IF( s1.LE.s2 ) THEN
287 s = zero
288 c = one
289 sestpr = s1
290 ELSE
291 s = one
292 c = zero
293 sestpr = s2
294 END IF
295 RETURN
296 ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam ) THEN
297 s1 = absgam
298 s2 = absalp
299 IF( s1.LE.s2 ) THEN
300 tmp = s1 / s2
301 c = sqrt( one+tmp*tmp )
302 sestpr = absest*( tmp / c )
303 s = -( gamma / s2 ) / c
304 c = sign( one, alpha ) / c
305 ELSE
306 tmp = s2 / s1
307 s = sqrt( one+tmp*tmp )
308 sestpr = absest / s
309 c = ( alpha / s1 ) / s
310 s = -sign( one, gamma ) / s
311 END IF
312 RETURN
313 ELSE
314*
315* normal case
316*
317 zeta1 = alpha / absest
318 zeta2 = gamma / absest
319*
320 norma = max( one+zeta1*zeta1+abs( zeta1*zeta2 ),
321 $ abs( zeta1*zeta2 )+zeta2*zeta2 )
322*
323* See if root is closer to zero or to ONE
324*
325 test = one + two*( zeta1-zeta2 )*( zeta1+zeta2 )
326 IF( test.GE.zero ) THEN
327*
328* root is close to zero, compute directly
329*
330 b = ( zeta1*zeta1+zeta2*zeta2+one )*half
331 c = zeta2*zeta2
332 t = c / ( b+sqrt( abs( b*b-c ) ) )
333 sine = zeta1 / ( one-t )
334 cosine = -zeta2 / t
335 sestpr = sqrt( t+four*eps*eps*norma )*absest
336 ELSE
337*
338* root is closer to ONE, shift by that amount
339*
340 b = ( zeta2*zeta2+zeta1*zeta1-one )*half
341 c = zeta1*zeta1
342 IF( b.GE.zero ) THEN
343 t = -c / ( b+sqrt( b*b+c ) )
344 ELSE
345 t = b - sqrt( b*b+c )
346 END IF
347 sine = -zeta1 / t
348 cosine = -zeta2 / ( one+t )
349 sestpr = sqrt( one+t+four*eps*eps*norma )*absest
350 END IF
351 tmp = sqrt( sine*sine+cosine*cosine )
352 s = sine / tmp
353 c = cosine / tmp
354 RETURN
355*
356 END IF
357 END IF
358 RETURN
359*
360* End of SLAIC1
361*
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
Here is the caller graph for this function: