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