LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ zlaic1()

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.

Definition at line 134 of file zlaic1.f.

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