LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the caller graph for this function: