Annotation of rpl/lapack/lapack/zlaic1.f, revision 1.19
1.12 bertrand 1: *> \brief \b ZLAIC1 applies one step of incremental condition estimation.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 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">
1.9 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
1.16 bertrand 22: *
1.9 bertrand 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: * ..
1.16 bertrand 31: *
1.9 bertrand 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: *
1.16 bertrand 126: *> \author Univ. of Tennessee
127: *> \author Univ. of California Berkeley
128: *> \author Univ. of Colorado Denver
129: *> \author NAG Ltd.
1.9 bertrand 130: *
131: *> \ingroup complex16OTHERauxiliary
132: *
133: * =====================================================================
1.1 bertrand 134: SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
135: *
1.19 ! bertrand 136: * -- LAPACK auxiliary routine --
1.1 bertrand 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: DOUBLE PRECISION SEST, SESTPR
143: COMPLEX*16 C, GAMMA, S
144: * ..
145: * .. Array Arguments ..
146: COMPLEX*16 W( J ), X( J )
147: * ..
148: *
149: * =====================================================================
150: *
151: * .. Parameters ..
152: DOUBLE PRECISION ZERO, ONE, TWO
153: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
154: DOUBLE PRECISION HALF, FOUR
155: PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 )
156: * ..
157: * .. Local Scalars ..
158: DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
159: $ SCL, T, TEST, TMP, ZETA1, ZETA2
160: COMPLEX*16 ALPHA, COSINE, SINE
161: * ..
162: * .. Intrinsic Functions ..
163: INTRINSIC ABS, DCONJG, MAX, SQRT
164: * ..
165: * .. External Functions ..
166: DOUBLE PRECISION DLAMCH
167: COMPLEX*16 ZDOTC
168: EXTERNAL DLAMCH, ZDOTC
169: * ..
170: * .. Executable Statements ..
171: *
172: EPS = DLAMCH( 'Epsilon' )
173: ALPHA = ZDOTC( 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
1.19 ! bertrand 194: TMP = DBLE( SQRT( S*DCONJG( S )+C*DCONJG( C ) ) )
1.1 bertrand 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
1.19 ! bertrand 248: T = DBLE( C / ( B+SQRT( B*B+C ) ) )
1.1 bertrand 249: ELSE
1.19 ! bertrand 250: T = DBLE( SQRT( B*B+C ) - B )
1.1 bertrand 251: END IF
252: *
253: SINE = -( ALPHA / ABSEST ) / T
254: COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
1.19 ! bertrand 255: TMP = DBLE( SQRT( SINE * DCONJG( SINE )
! 256: $ + COSINE * DCONJG( COSINE ) ) )
! 257:
1.1 bertrand 258: S = SINE / TMP
259: C = COSINE / TMP
260: SESTPR = SQRT( T+ONE )*ABSEST
261: RETURN
262: END IF
263: *
264: ELSE IF( JOB.EQ.2 ) THEN
265: *
266: * Estimating smallest singular value
267: *
268: * special cases
269: *
270: IF( SEST.EQ.ZERO ) THEN
271: SESTPR = ZERO
272: IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
273: SINE = ONE
274: COSINE = ZERO
275: ELSE
276: SINE = -DCONJG( GAMMA )
277: COSINE = DCONJG( ALPHA )
278: END IF
279: S1 = MAX( ABS( SINE ), ABS( COSINE ) )
280: S = SINE / S1
281: C = COSINE / S1
1.19 ! bertrand 282: TMP = DBLE( SQRT( S*DCONJG( S )+C*DCONJG( C ) ) )
1.1 bertrand 283: S = S / TMP
284: C = C / TMP
285: RETURN
286: ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
287: S = ZERO
288: C = ONE
289: SESTPR = ABSGAM
290: RETURN
291: ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
292: S1 = ABSGAM
293: S2 = ABSEST
294: IF( S1.LE.S2 ) THEN
295: S = ZERO
296: C = ONE
297: SESTPR = S1
298: ELSE
299: S = ONE
300: C = ZERO
301: SESTPR = S2
302: END IF
303: RETURN
304: ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
305: S1 = ABSGAM
306: S2 = ABSALP
307: IF( S1.LE.S2 ) THEN
308: TMP = S1 / S2
309: SCL = SQRT( ONE+TMP*TMP )
310: SESTPR = ABSEST*( TMP / SCL )
311: S = -( DCONJG( GAMMA ) / S2 ) / SCL
312: C = ( DCONJG( ALPHA ) / S2 ) / SCL
313: ELSE
314: TMP = S2 / S1
315: SCL = SQRT( ONE+TMP*TMP )
316: SESTPR = ABSEST / SCL
317: S = -( DCONJG( GAMMA ) / S1 ) / SCL
318: C = ( DCONJG( ALPHA ) / S1 ) / SCL
319: END IF
320: RETURN
321: ELSE
322: *
323: * normal case
324: *
325: ZETA1 = ABSALP / ABSEST
326: ZETA2 = ABSGAM / ABSEST
327: *
328: NORMA = MAX( ONE+ZETA1*ZETA1+ZETA1*ZETA2,
329: $ ZETA1*ZETA2+ZETA2*ZETA2 )
330: *
331: * See if root is closer to zero or to ONE
332: *
333: TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
334: IF( TEST.GE.ZERO ) THEN
335: *
336: * root is close to zero, compute directly
337: *
338: B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
339: C = ZETA2*ZETA2
1.19 ! bertrand 340: T = DBLE( C / ( B+SQRT( ABS( B*B-C ) ) ) )
1.1 bertrand 341: SINE = ( ALPHA / ABSEST ) / ( ONE-T )
342: COSINE = -( GAMMA / ABSEST ) / T
343: SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
344: ELSE
345: *
346: * root is closer to ONE, shift by that amount
347: *
348: B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
349: C = ZETA1*ZETA1
350: IF( B.GE.ZERO ) THEN
1.19 ! bertrand 351: T = DBLE( -C / ( B+SQRT( B*B+C ) ) )
1.1 bertrand 352: ELSE
1.19 ! bertrand 353: T = DBLE( B - SQRT( B*B+C ) )
1.1 bertrand 354: END IF
355: SINE = -( ALPHA / ABSEST ) / T
356: COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
357: SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
358: END IF
1.19 ! bertrand 359: TMP = DBLE( SQRT( SINE * DCONJG( SINE )
! 360: $ + COSINE * DCONJG( COSINE ) ) )
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[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