![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaic1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zlaic1.f CVS log | Up to [local] / rpl / lapack / lapack |
1.12 bertrand 1: *> \brief \b ZLAIC1 applies one step of incremental condition estimation.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLAIC1 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaic1.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaic1.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaic1.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER J, JOB
25: * DOUBLE PRECISION SEST, SESTPR
26: * COMPLEX*16 C, GAMMA, S
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 W( J ), X( J )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZLAIC1 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 ZLAIC1 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*16 array, dimension (J)
84: *> The j-vector x.
85: *> \endverbatim
86: *>
87: *> \param[in] SEST
88: *> \verbatim
89: *> SEST is DOUBLE PRECISION
90: *> Estimated singular value of j by j matrix L
91: *> \endverbatim
92: *>
93: *> \param[in] W
94: *> \verbatim
95: *> W is COMPLEX*16 array, dimension (J)
96: *> The j-vector w.
97: *> \endverbatim
98: *>
99: *> \param[in] GAMMA
100: *> \verbatim
101: *> GAMMA is COMPLEX*16
102: *> The diagonal element gamma.
103: *> \endverbatim
104: *>
105: *> \param[out] SESTPR
106: *> \verbatim
107: *> SESTPR is DOUBLE PRECISION
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*16
114: *> Sine needed in forming xhat.
115: *> \endverbatim
116: *>
117: *> \param[out] C
118: *> \verbatim
119: *> C is COMPLEX*16
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: *
1.12 bertrand 131: *> \date September 2012
1.9 bertrand 132: *
133: *> \ingroup complex16OTHERauxiliary
134: *
135: * =====================================================================
1.1 bertrand 136: SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
137: *
1.12 bertrand 138: * -- LAPACK auxiliary routine (version 3.4.2) --
1.1 bertrand 139: * -- LAPACK is a software package provided by Univ. of Tennessee, --
140: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.12 bertrand 141: * September 2012
1.1 bertrand 142: *
143: * .. Scalar Arguments ..
144: INTEGER J, JOB
145: DOUBLE PRECISION SEST, SESTPR
146: COMPLEX*16 C, GAMMA, S
147: * ..
148: * .. Array Arguments ..
149: COMPLEX*16 W( J ), X( J )
150: * ..
151: *
152: * =====================================================================
153: *
154: * .. Parameters ..
155: DOUBLE PRECISION ZERO, ONE, TWO
156: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
157: DOUBLE PRECISION HALF, FOUR
158: PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 )
159: * ..
160: * .. Local Scalars ..
161: DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
162: $ SCL, T, TEST, TMP, ZETA1, ZETA2
163: COMPLEX*16 ALPHA, COSINE, SINE
164: * ..
165: * .. Intrinsic Functions ..
166: INTRINSIC ABS, DCONJG, MAX, SQRT
167: * ..
168: * .. External Functions ..
169: DOUBLE PRECISION DLAMCH
170: COMPLEX*16 ZDOTC
171: EXTERNAL DLAMCH, ZDOTC
172: * ..
173: * .. Executable Statements ..
174: *
175: EPS = DLAMCH( 'Epsilon' )
176: ALPHA = ZDOTC( 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*DCONJG( S )+C*DCONJG( 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*DCONJG( SINE )+COSINE*DCONJG( 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 = -DCONJG( GAMMA )
278: COSINE = DCONJG( ALPHA )
279: END IF
280: S1 = MAX( ABS( SINE ), ABS( COSINE ) )
281: S = SINE / S1
282: C = COSINE / S1
283: TMP = SQRT( S*DCONJG( S )+C*DCONJG( 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 = -( DCONJG( GAMMA ) / S2 ) / SCL
313: C = ( DCONJG( ALPHA ) / S2 ) / SCL
314: ELSE
315: TMP = S2 / S1
316: SCL = SQRT( ONE+TMP*TMP )
317: SESTPR = ABSEST / SCL
318: S = -( DCONJG( GAMMA ) / S1 ) / SCL
319: C = ( DCONJG( 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*DCONJG( SINE )+COSINE*DCONJG( COSINE ) )
361: S = SINE / TMP
362: C = COSINE / TMP
363: RETURN
364: *
365: END IF
366: END IF
367: RETURN
368: *
369: * End of ZLAIC1
370: *
371: END