LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 136 of file dlaic1.f.

136 *
137 * -- LAPACK auxiliary routine (version 3.4.2) --
138 * -- LAPACK is a software package provided by Univ. of Tennessee, --
139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140 * September 2012
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 *
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53

Here is the caller graph for this function: