LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zlaic1 ( integer  JOB,
integer  J,
complex*16, dimension( j )  X,
double precision  SEST,
complex*16, dimension( j )  W,
complex*16  GAMMA,
double precision  SESTPR,
complex*16  S,
complex*16  C 
)

ZLAIC1 applies one step of incremental condition estimation.

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

Purpose: 
 ZLAIC1 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 ZLAIC1 computes sestpr, s, c such that
 the vector
                 [ s*x ]
          xhat = [  c  ]
 is an approximate singular vector of
                 [ L       0  ]
          Lhat = [ w**H 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]**H and sestpr**2 is an eigenpair of the system

     diag(sest*sest, 0) + [alpha  gamma] * [ conjg(alpha) ]
                                           [ conjg(gamma) ]

 where  alpha =  x**H * 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 COMPLEX*16 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 COMPLEX*16 array, dimension (J)
          The j-vector w.
[in]GAMMA
          GAMMA is COMPLEX*16
          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 COMPLEX*16
          Sine needed in forming xhat.
[out]C
          C is COMPLEX*16
          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 137 of file zlaic1.f.

137 *
138 * -- LAPACK auxiliary routine (version 3.4.2) --
139 * -- LAPACK is a software package provided by Univ. of Tennessee, --
140 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
141 * September 2012
142 *
143 * .. Scalar Arguments ..
144  INTEGER j, job
145  DOUBLE PRECISION sest, sestpr
146  COMPLEX*16 c, gamma, s
147 * ..
148 * .. Array Arguments ..
149  COMPLEX*16 w( j ), x( j )
150 * ..
151 *
152 * =====================================================================
153 *
154 * .. Parameters ..
155  DOUBLE PRECISION zero, one, two
156  parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
157  DOUBLE PRECISION half, four
158  parameter ( half = 0.5d0, four = 4.0d0 )
159 * ..
160 * .. Local Scalars ..
161  DOUBLE PRECISION absalp, absest, absgam, b, eps, norma, s1, s2,
162  $ scl, t, test, tmp, zeta1, zeta2
163  COMPLEX*16 alpha, cosine, sine
164 * ..
165 * .. Intrinsic Functions ..
166  INTRINSIC abs, dconjg, max, sqrt
167 * ..
168 * .. External Functions ..
169  DOUBLE PRECISION dlamch
170  COMPLEX*16 zdotc
171  EXTERNAL dlamch, zdotc
172 * ..
173 * .. Executable Statements ..
174 *
175  eps = dlamch( 'Epsilon' )
176  alpha = zdotc( j, x, 1, w, 1 )
177 *
178  absalp = abs( alpha )
179  absgam = abs( gamma )
180  absest = abs( sest )
181 *
182  IF( job.EQ.1 ) THEN
183 *
184 * Estimating largest singular value
185 *
186 * special cases
187 *
188  IF( sest.EQ.zero ) THEN
189  s1 = max( absgam, absalp )
190  IF( s1.EQ.zero ) THEN
191  s = zero
192  c = one
193  sestpr = zero
194  ELSE
195  s = alpha / s1
196  c = gamma / s1
197  tmp = sqrt( s*dconjg( s )+c*dconjg( c ) )
198  s = s / tmp
199  c = c / tmp
200  sestpr = s1*tmp
201  END IF
202  RETURN
203  ELSE IF( absgam.LE.eps*absest ) THEN
204  s = one
205  c = zero
206  tmp = max( absest, absalp )
207  s1 = absest / tmp
208  s2 = absalp / tmp
209  sestpr = tmp*sqrt( s1*s1+s2*s2 )
210  RETURN
211  ELSE IF( absalp.LE.eps*absest ) THEN
212  s1 = absgam
213  s2 = absest
214  IF( s1.LE.s2 ) THEN
215  s = one
216  c = zero
217  sestpr = s2
218  ELSE
219  s = zero
220  c = one
221  sestpr = s1
222  END IF
223  RETURN
224  ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam ) THEN
225  s1 = absgam
226  s2 = absalp
227  IF( s1.LE.s2 ) THEN
228  tmp = s1 / s2
229  scl = sqrt( one+tmp*tmp )
230  sestpr = s2*scl
231  s = ( alpha / s2 ) / scl
232  c = ( gamma / s2 ) / scl
233  ELSE
234  tmp = s2 / s1
235  scl = sqrt( one+tmp*tmp )
236  sestpr = s1*scl
237  s = ( alpha / s1 ) / scl
238  c = ( gamma / s1 ) / scl
239  END IF
240  RETURN
241  ELSE
242 *
243 * normal case
244 *
245  zeta1 = absalp / absest
246  zeta2 = absgam / absest
247 *
248  b = ( one-zeta1*zeta1-zeta2*zeta2 )*half
249  c = zeta1*zeta1
250  IF( b.GT.zero ) THEN
251  t = c / ( b+sqrt( b*b+c ) )
252  ELSE
253  t = sqrt( b*b+c ) - b
254  END IF
255 *
256  sine = -( alpha / absest ) / t
257  cosine = -( gamma / absest ) / ( one+t )
258  tmp = sqrt( sine*dconjg( sine )+cosine*dconjg( cosine ) )
259  s = sine / tmp
260  c = cosine / tmp
261  sestpr = sqrt( t+one )*absest
262  RETURN
263  END IF
264 *
265  ELSE IF( job.EQ.2 ) THEN
266 *
267 * Estimating smallest singular value
268 *
269 * special cases
270 *
271  IF( sest.EQ.zero ) THEN
272  sestpr = zero
273  IF( max( absgam, absalp ).EQ.zero ) THEN
274  sine = one
275  cosine = zero
276  ELSE
277  sine = -dconjg( gamma )
278  cosine = dconjg( alpha )
279  END IF
280  s1 = max( abs( sine ), abs( cosine ) )
281  s = sine / s1
282  c = cosine / s1
283  tmp = sqrt( s*dconjg( s )+c*dconjg( c ) )
284  s = s / tmp
285  c = c / tmp
286  RETURN
287  ELSE IF( absgam.LE.eps*absest ) THEN
288  s = zero
289  c = one
290  sestpr = absgam
291  RETURN
292  ELSE IF( absalp.LE.eps*absest ) THEN
293  s1 = absgam
294  s2 = absest
295  IF( s1.LE.s2 ) THEN
296  s = zero
297  c = one
298  sestpr = s1
299  ELSE
300  s = one
301  c = zero
302  sestpr = s2
303  END IF
304  RETURN
305  ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam ) THEN
306  s1 = absgam
307  s2 = absalp
308  IF( s1.LE.s2 ) THEN
309  tmp = s1 / s2
310  scl = sqrt( one+tmp*tmp )
311  sestpr = absest*( tmp / scl )
312  s = -( dconjg( gamma ) / s2 ) / scl
313  c = ( dconjg( alpha ) / s2 ) / scl
314  ELSE
315  tmp = s2 / s1
316  scl = sqrt( one+tmp*tmp )
317  sestpr = absest / scl
318  s = -( dconjg( gamma ) / s1 ) / scl
319  c = ( dconjg( alpha ) / s1 ) / scl
320  END IF
321  RETURN
322  ELSE
323 *
324 * normal case
325 *
326  zeta1 = absalp / absest
327  zeta2 = absgam / absest
328 *
329  norma = max( one+zeta1*zeta1+zeta1*zeta2,
330  $ zeta1*zeta2+zeta2*zeta2 )
331 *
332 * See if root is closer to zero or to ONE
333 *
334  test = one + two*( zeta1-zeta2 )*( zeta1+zeta2 )
335  IF( test.GE.zero ) THEN
336 *
337 * root is close to zero, compute directly
338 *
339  b = ( zeta1*zeta1+zeta2*zeta2+one )*half
340  c = zeta2*zeta2
341  t = c / ( b+sqrt( abs( b*b-c ) ) )
342  sine = ( alpha / absest ) / ( one-t )
343  cosine = -( gamma / absest ) / t
344  sestpr = sqrt( t+four*eps*eps*norma )*absest
345  ELSE
346 *
347 * root is closer to ONE, shift by that amount
348 *
349  b = ( zeta2*zeta2+zeta1*zeta1-one )*half
350  c = zeta1*zeta1
351  IF( b.GE.zero ) THEN
352  t = -c / ( b+sqrt( b*b+c ) )
353  ELSE
354  t = b - sqrt( b*b+c )
355  END IF
356  sine = -( alpha / absest ) / t
357  cosine = -( gamma / absest ) / ( one+t )
358  sestpr = sqrt( one+t+four*eps*eps*norma )*absest
359  END IF
360  tmp = sqrt( sine*dconjg( sine )+cosine*dconjg( cosine ) )
361  s = sine / tmp
362  c = cosine / tmp
363  RETURN
364 *
365  END IF
366  END IF
367  RETURN
368 *
369 * End of ZLAIC1
370 *
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zdotc
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:54

Here is the caller graph for this function: