 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dlaic1()

 subroutine dlaic1 ( integer JOB, integer J, double precision, dimension( j ) X, double precision SEST, double precision, dimension( j ) W, double precision GAMMA, double precision SESTPR, double precision S, double precision C )

DLAIC1 applies one step of incremental condition estimation.

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

Purpose:
``` DLAIC1 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 DLAIC1 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 DOUBLE PRECISION array, dimension (J) The j-vector x.``` [in] SEST ``` SEST is DOUBLE PRECISION Estimated singular value of j by j matrix L``` [in] W ``` W is DOUBLE PRECISION array, dimension (J) The j-vector w.``` [in] GAMMA ``` GAMMA is DOUBLE PRECISION The diagonal element gamma.``` [out] SESTPR ``` SESTPR is DOUBLE PRECISION Estimated singular value of (j+1) by (j+1) matrix Lhat.``` [out] S ``` S is DOUBLE PRECISION Sine needed in forming xhat.``` [out] C ``` C is DOUBLE PRECISION Cosine needed in forming xhat.```

Definition at line 133 of file dlaic1.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  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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
Here is the caller graph for this function: