LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
claic1.f
Go to the documentation of this file.
1 *> \brief \b CLAIC1 applies one step of incremental condition estimation.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CLAIC1 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claic1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claic1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claic1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER J, JOB
25 * REAL SEST, SESTPR
26 * COMPLEX C, GAMMA, S
27 * ..
28 * .. Array Arguments ..
29 * COMPLEX W( J ), X( J )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> CLAIC1 applies one step of incremental condition estimation in
39 *> its simplest version:
40 *>
41 *> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
42 *> lower triangular matrix L, such that
43 *> twonorm(L*x) = sest
44 *> Then CLAIC1 computes sestpr, s, c such that
45 *> the vector
46 *> [ s*x ]
47 *> xhat = [ c ]
48 *> is an approximate singular vector of
49 *> [ L 0 ]
50 *> Lhat = [ w**H gamma ]
51 *> in the sense that
52 *> twonorm(Lhat*xhat) = sestpr.
53 *>
54 *> Depending on JOB, an estimate for the largest or smallest singular
55 *> value is computed.
56 *>
57 *> Note that [s c]**H and sestpr**2 is an eigenpair of the system
58 *>
59 *> diag(sest*sest, 0) + [alpha gamma] * [ conjg(alpha) ]
60 *> [ conjg(gamma) ]
61 *>
62 *> where alpha = x**H*w.
63 *> \endverbatim
64 *
65 * Arguments:
66 * ==========
67 *
68 *> \param[in] JOB
69 *> \verbatim
70 *> JOB is INTEGER
71 *> = 1: an estimate for the largest singular value is computed.
72 *> = 2: an estimate for the smallest singular value is computed.
73 *> \endverbatim
74 *>
75 *> \param[in] J
76 *> \verbatim
77 *> J is INTEGER
78 *> Length of X and W
79 *> \endverbatim
80 *>
81 *> \param[in] X
82 *> \verbatim
83 *> X is COMPLEX array, dimension (J)
84 *> The j-vector x.
85 *> \endverbatim
86 *>
87 *> \param[in] SEST
88 *> \verbatim
89 *> SEST is REAL
90 *> Estimated singular value of j by j matrix L
91 *> \endverbatim
92 *>
93 *> \param[in] W
94 *> \verbatim
95 *> W is COMPLEX array, dimension (J)
96 *> The j-vector w.
97 *> \endverbatim
98 *>
99 *> \param[in] GAMMA
100 *> \verbatim
101 *> GAMMA is COMPLEX
102 *> The diagonal element gamma.
103 *> \endverbatim
104 *>
105 *> \param[out] SESTPR
106 *> \verbatim
107 *> SESTPR is REAL
108 *> Estimated singular value of (j+1) by (j+1) matrix Lhat.
109 *> \endverbatim
110 *>
111 *> \param[out] S
112 *> \verbatim
113 *> S is COMPLEX
114 *> Sine needed in forming xhat.
115 *> \endverbatim
116 *>
117 *> \param[out] C
118 *> \verbatim
119 *> C is COMPLEX
120 *> Cosine needed in forming xhat.
121 *> \endverbatim
122 *
123 * Authors:
124 * ========
125 *
126 *> \author Univ. of Tennessee
127 *> \author Univ. of California Berkeley
128 *> \author Univ. of Colorado Denver
129 *> \author NAG Ltd.
130 *
131 *> \ingroup complexOTHERauxiliary
132 *
133 * =====================================================================
134  SUBROUTINE claic1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
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  REAL SEST, SESTPR
143  COMPLEX C, GAMMA, S
144 * ..
145 * .. Array Arguments ..
146  COMPLEX W( J ), X( J )
147 * ..
148 *
149 * =====================================================================
150 *
151 * .. Parameters ..
152  REAL ZERO, ONE, TWO
153  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
154  REAL HALF, FOUR
155  parameter( half = 0.5e0, four = 4.0e0 )
156 * ..
157 * .. Local Scalars ..
158  REAL ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
159  $ SCL, T, TEST, TMP, ZETA1, ZETA2
160  COMPLEX ALPHA, COSINE, SINE
161 * ..
162 * .. Intrinsic Functions ..
163  INTRINSIC abs, conjg, max, sqrt
164 * ..
165 * .. External Functions ..
166  REAL SLAMCH
167  COMPLEX CDOTC
168  EXTERNAL slamch, cdotc
169 * ..
170 * .. Executable Statements ..
171 *
172  eps = slamch( 'Epsilon' )
173  alpha = cdotc( 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 = real( sqrt( s*conjg( s )+c*conjg( 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 = real( c / ( b+sqrt( b*b+c ) ) )
249  ELSE
250  t = real( sqrt( b*b+c ) - b )
251  END IF
252 *
253  sine = -( alpha / absest ) / t
254  cosine = -( gamma / absest ) / ( one+t )
255  tmp = real( sqrt( sine * conjg( sine )
256  $ + cosine * conjg( cosine ) ) )
257  s = sine / tmp
258  c = cosine / tmp
259  sestpr = sqrt( t+one )*absest
260  RETURN
261  END IF
262 *
263  ELSE IF( job.EQ.2 ) THEN
264 *
265 * Estimating smallest singular value
266 *
267 * special cases
268 *
269  IF( sest.EQ.zero ) THEN
270  sestpr = zero
271  IF( max( absgam, absalp ).EQ.zero ) THEN
272  sine = one
273  cosine = zero
274  ELSE
275  sine = -conjg( gamma )
276  cosine = conjg( alpha )
277  END IF
278  s1 = max( abs( sine ), abs( cosine ) )
279  s = sine / s1
280  c = cosine / s1
281  tmp = real( sqrt( s*conjg( s )+c*conjg( c ) ) )
282  s = s / tmp
283  c = c / tmp
284  RETURN
285  ELSE IF( absgam.LE.eps*absest ) THEN
286  s = zero
287  c = one
288  sestpr = absgam
289  RETURN
290  ELSE IF( absalp.LE.eps*absest ) THEN
291  s1 = absgam
292  s2 = absest
293  IF( s1.LE.s2 ) THEN
294  s = zero
295  c = one
296  sestpr = s1
297  ELSE
298  s = one
299  c = zero
300  sestpr = s2
301  END IF
302  RETURN
303  ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam ) THEN
304  s1 = absgam
305  s2 = absalp
306  IF( s1.LE.s2 ) THEN
307  tmp = s1 / s2
308  scl = sqrt( one+tmp*tmp )
309  sestpr = absest*( tmp / scl )
310  s = -( conjg( gamma ) / s2 ) / scl
311  c = ( conjg( alpha ) / s2 ) / scl
312  ELSE
313  tmp = s2 / s1
314  scl = sqrt( one+tmp*tmp )
315  sestpr = absest / scl
316  s = -( conjg( gamma ) / s1 ) / scl
317  c = ( conjg( alpha ) / s1 ) / scl
318  END IF
319  RETURN
320  ELSE
321 *
322 * normal case
323 *
324  zeta1 = absalp / absest
325  zeta2 = absgam / absest
326 *
327  norma = max( one+zeta1*zeta1+zeta1*zeta2,
328  $ zeta1*zeta2+zeta2*zeta2 )
329 *
330 * See if root is closer to zero or to ONE
331 *
332  test = one + two*( zeta1-zeta2 )*( zeta1+zeta2 )
333  IF( test.GE.zero ) THEN
334 *
335 * root is close to zero, compute directly
336 *
337  b = ( zeta1*zeta1+zeta2*zeta2+one )*half
338  c = zeta2*zeta2
339  t = real( c / ( b+sqrt( abs( b*b-c ) ) ) )
340  sine = ( alpha / absest ) / ( one-t )
341  cosine = -( gamma / absest ) / t
342  sestpr = sqrt( t+four*eps*eps*norma )*absest
343  ELSE
344 *
345 * root is closer to ONE, shift by that amount
346 *
347  b = ( zeta2*zeta2+zeta1*zeta1-one )*half
348  c = zeta1*zeta1
349  IF( b.GE.zero ) THEN
350  t = real( -c / ( b+sqrt( b*b+c ) ) )
351  ELSE
352  t = real( b - sqrt( b*b+c ) )
353  END IF
354  sine = -( alpha / absest ) / t
355  cosine = -( gamma / absest ) / ( one+t )
356  sestpr = sqrt( one+t+four*eps*eps*norma )*absest
357  END IF
358  tmp = real( sqrt( sine * conjg( sine )
359  $ + cosine * conjg( cosine ) ) )
360  s = sine / tmp
361  c = cosine / tmp
362  RETURN
363 *
364  END IF
365  END IF
366  RETURN
367 *
368 * End of CLAIC1
369 *
370  END
subroutine claic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
CLAIC1 applies one step of incremental condition estimation.
Definition: claic1.f:135