LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine claic1 ( integer  JOB,
integer  J,
complex, dimension( j )  X,
real  SEST,
complex, dimension( j )  W,
complex  GAMMA,
real  SESTPR,
complex  S,
complex  C 
)

CLAIC1 applies one step of incremental condition estimation.

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

Purpose: 
 CLAIC1 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 CLAIC1 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 array, dimension (J)
          The j-vector x.
[in]SEST
          SEST is REAL
          Estimated singular value of j by j matrix L
[in]W
          W is COMPLEX array, dimension (J)
          The j-vector w.
[in]GAMMA
          GAMMA is COMPLEX
          The diagonal element gamma.
[out]SESTPR
          SESTPR is REAL
          Estimated singular value of (j+1) by (j+1) matrix Lhat.
[out]S
          S is COMPLEX
          Sine needed in forming xhat.
[out]C
          C is COMPLEX
          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 claic1.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  REAL sest, sestpr
146  COMPLEX c, gamma, s
147 * ..
148 * .. Array Arguments ..
149  COMPLEX w( j ), x( j )
150 * ..
151 *
152 * =====================================================================
153 *
154 * .. Parameters ..
155  REAL zero, one, two
156  parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
157  REAL half, four
158  parameter ( half = 0.5e0, four = 4.0e0 )
159 * ..
160 * .. Local Scalars ..
161  REAL absalp, absest, absgam, b, eps, norma, s1, s2,
162  $ scl, t, test, tmp, zeta1, zeta2
163  COMPLEX alpha, cosine, sine
164 * ..
165 * .. Intrinsic Functions ..
166  INTRINSIC abs, conjg, max, sqrt
167 * ..
168 * .. External Functions ..
169  REAL slamch
170  COMPLEX cdotc
171  EXTERNAL slamch, cdotc
172 * ..
173 * .. Executable Statements ..
174 *
175  eps = slamch( 'Epsilon' )
176  alpha = cdotc( 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*conjg( s )+c*conjg( 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*conjg( sine )+cosine*conjg( 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 = -conjg( gamma )
278  cosine = conjg( alpha )
279  END IF
280  s1 = max( abs( sine ), abs( cosine ) )
281  s = sine / s1
282  c = cosine / s1
283  tmp = sqrt( s*conjg( s )+c*conjg( 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 = -( conjg( gamma ) / s2 ) / scl
313  c = ( conjg( alpha ) / s2 ) / scl
314  ELSE
315  tmp = s2 / s1
316  scl = sqrt( one+tmp*tmp )
317  sestpr = absest / scl
318  s = -( conjg( gamma ) / s1 ) / scl
319  c = ( conjg( 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*conjg( sine )+cosine*conjg( cosine ) )
361  s = sine / tmp
362  c = cosine / tmp
363  RETURN
364 *
365  END IF
366  END IF
367  RETURN
368 *
369 * End of CLAIC1
370 *
cdotc
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:54
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69

Here is the caller graph for this function: