LAPACK 3.12.0
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 133 of file slaic1.f.

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 REAL C, GAMMA, S, SEST, SESTPR
142* ..
143* .. Array Arguments ..
144 REAL 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, 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 REAL SDOT, SLAMCH
164 EXTERNAL sdot, slamch
165* ..
166* .. Executable Statements ..
167*
168 eps = slamch( 'Epsilon' )
169 alpha = sdot( 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 SLAIC1
363*
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: