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