 LAPACK  3.8.0 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.```
Date
December 2016

Definition at line 136 of file dlaic1.f.

136 *
137 * -- LAPACK auxiliary routine (version 3.7.0) --
138 * -- LAPACK is a software package provided by Univ. of Tennessee, --
139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140 * December 2016
141 *
142 * .. Scalar Arguments ..
143  INTEGER j, job
144  DOUBLE PRECISION c, gamma, s, sest, sestpr
145 * ..
146 * .. Array Arguments ..
147  DOUBLE PRECISION w( j ), x( j )
148 * ..
149 *
150 * =====================================================================
151 *
152 * .. Parameters ..
153  DOUBLE PRECISION zero, one, two
154  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
155  DOUBLE PRECISION half, four
156  parameter( half = 0.5d0, four = 4.0d0 )
157 * ..
158 * .. Local Scalars ..
159  DOUBLE PRECISION absalp, absest, absgam, alpha, b, cosine, eps,
160  \$ norma, s1, s2, sine, t, test, tmp, zeta1, zeta2
161 * ..
162 * .. Intrinsic Functions ..
163  INTRINSIC abs, max, sign, sqrt
164 * ..
165 * .. External Functions ..
166  DOUBLE PRECISION ddot, dlamch
167  EXTERNAL ddot, dlamch
168 * ..
169 * .. Executable Statements ..
170 *
171  eps = dlamch( 'Epsilon' )
172  alpha = ddot( j, x, 1, w, 1 )
173 *
174  absalp = abs( alpha )
175  absgam = abs( gamma )
176  absest = abs( sest )
177 *
178  IF( job.EQ.1 ) THEN
179 *
180 * Estimating largest singular value
181 *
182 * special cases
183 *
184  IF( sest.EQ.zero ) THEN
185  s1 = max( absgam, absalp )
186  IF( s1.EQ.zero ) THEN
187  s = zero
188  c = one
189  sestpr = zero
190  ELSE
191  s = alpha / s1
192  c = gamma / s1
193  tmp = sqrt( s*s+c*c )
194  s = s / tmp
195  c = c / tmp
196  sestpr = s1*tmp
197  END IF
198  RETURN
199  ELSE IF( absgam.LE.eps*absest ) THEN
200  s = one
201  c = zero
202  tmp = max( absest, absalp )
203  s1 = absest / tmp
204  s2 = absalp / tmp
205  sestpr = tmp*sqrt( s1*s1+s2*s2 )
206  RETURN
207  ELSE IF( absalp.LE.eps*absest ) THEN
208  s1 = absgam
209  s2 = absest
210  IF( s1.LE.s2 ) THEN
211  s = one
212  c = zero
213  sestpr = s2
214  ELSE
215  s = zero
216  c = one
217  sestpr = s1
218  END IF
219  RETURN
220  ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam ) THEN
221  s1 = absgam
222  s2 = absalp
223  IF( s1.LE.s2 ) THEN
224  tmp = s1 / s2
225  s = sqrt( one+tmp*tmp )
226  sestpr = s2*s
227  c = ( gamma / s2 ) / s
228  s = sign( one, alpha ) / s
229  ELSE
230  tmp = s2 / s1
231  c = sqrt( one+tmp*tmp )
232  sestpr = s1*c
233  s = ( alpha / s1 ) / c
234  c = sign( one, gamma ) / c
235  END IF
236  RETURN
237  ELSE
238 *
239 * normal case
240 *
241  zeta1 = alpha / absest
242  zeta2 = gamma / absest
243 *
244  b = ( one-zeta1*zeta1-zeta2*zeta2 )*half
245  c = zeta1*zeta1
246  IF( b.GT.zero ) THEN
247  t = c / ( b+sqrt( b*b+c ) )
248  ELSE
249  t = sqrt( b*b+c ) - b
250  END IF
251 *
252  sine = -zeta1 / t
253  cosine = -zeta2 / ( one+t )
254  tmp = sqrt( sine*sine+cosine*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 = -gamma
274  cosine = alpha
275  END IF
276  s1 = max( abs( sine ), abs( cosine ) )
277  s = sine / s1
278  c = cosine / s1
279  tmp = sqrt( s*s+c*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  c = sqrt( one+tmp*tmp )
307  sestpr = absest*( tmp / c )
308  s = -( gamma / s2 ) / c
309  c = sign( one, alpha ) / c
310  ELSE
311  tmp = s2 / s1
312  s = sqrt( one+tmp*tmp )
313  sestpr = absest / s
314  c = ( alpha / s1 ) / s
315  s = -sign( one, gamma ) / s
316  END IF
317  RETURN
318  ELSE
319 *
320 * normal case
321 *
322  zeta1 = alpha / absest
323  zeta2 = gamma / absest
324 *
325  norma = max( one+zeta1*zeta1+abs( zeta1*zeta2 ),
326  \$ abs( 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 = c / ( b+sqrt( abs( b*b-c ) ) )
338  sine = zeta1 / ( one-t )
339  cosine = -zeta2 / 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 = -c / ( b+sqrt( b*b+c ) )
349  ELSE
350  t = b - sqrt( b*b+c )
351  END IF
352  sine = -zeta1 / t
353  cosine = -zeta2 / ( one+t )
354  sestpr = sqrt( one+t+four*eps*eps*norma )*absest
355  END IF
356  tmp = sqrt( sine*sine+cosine*cosine )
357  s = sine / tmp
358  c = cosine / tmp
359  RETURN
360 *
361  END IF
362  END IF
363  RETURN
364 *
365 * End of DLAIC1
366 *
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:84
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
Here is the caller graph for this function: