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: *> \ingroup doubleOTHERauxiliary
131: *
132: * =====================================================================
133: SUBROUTINE DLAIC1( 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: DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR
142: * ..
143: * .. Array Arguments ..
144: DOUBLE PRECISION W( J ), X( J )
145: * ..
146: *
147: * =====================================================================
148: *
149: * .. Parameters ..
150: DOUBLE PRECISION ZERO, ONE, TWO
151: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
152: DOUBLE PRECISION HALF, FOUR
153: PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 )
154: * ..
155: * .. Local Scalars ..
156: DOUBLE PRECISION 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: DOUBLE PRECISION DDOT, DLAMCH
164: EXTERNAL DDOT, DLAMCH
165: * ..
166: * .. Executable Statements ..
167: *
168: EPS = DLAMCH( 'Epsilon' )
169: ALPHA = DDOT( 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 DLAIC1
363: *
364: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaic1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack