![[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.1 bertrand 1: SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER J, JOB
10: DOUBLE PRECISION SEST, SESTPR
11: COMPLEX*16 C, GAMMA, S
12: * ..
13: * .. Array Arguments ..
14: COMPLEX*16 W( J ), X( J )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * ZLAIC1 applies one step of incremental condition estimation in
21: * its simplest version:
22: *
23: * Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
24: * lower triangular matrix L, such that
25: * twonorm(L*x) = sest
26: * Then ZLAIC1 computes sestpr, s, c such that
27: * the vector
28: * [ s*x ]
29: * xhat = [ c ]
30: * is an approximate singular vector of
31: * [ L 0 ]
32: * Lhat = [ w' gamma ]
33: * in the sense that
34: * twonorm(Lhat*xhat) = sestpr.
35: *
36: * Depending on JOB, an estimate for the largest or smallest singular
37: * value is computed.
38: *
39: * Note that [s c]' and sestpr**2 is an eigenpair of the system
40: *
41: * diag(sest*sest, 0) + [alpha gamma] * [ conjg(alpha) ]
42: * [ conjg(gamma) ]
43: *
44: * where alpha = conjg(x)'*w.
45: *
46: * Arguments
47: * =========
48: *
49: * JOB (input) INTEGER
50: * = 1: an estimate for the largest singular value is computed.
51: * = 2: an estimate for the smallest singular value is computed.
52: *
53: * J (input) INTEGER
54: * Length of X and W
55: *
56: * X (input) COMPLEX*16 array, dimension (J)
57: * The j-vector x.
58: *
59: * SEST (input) DOUBLE PRECISION
60: * Estimated singular value of j by j matrix L
61: *
62: * W (input) COMPLEX*16 array, dimension (J)
63: * The j-vector w.
64: *
65: * GAMMA (input) COMPLEX*16
66: * The diagonal element gamma.
67: *
68: * SESTPR (output) DOUBLE PRECISION
69: * Estimated singular value of (j+1) by (j+1) matrix Lhat.
70: *
71: * S (output) COMPLEX*16
72: * Sine needed in forming xhat.
73: *
74: * C (output) COMPLEX*16
75: * Cosine needed in forming xhat.
76: *
77: * =====================================================================
78: *
79: * .. Parameters ..
80: DOUBLE PRECISION ZERO, ONE, TWO
81: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
82: DOUBLE PRECISION HALF, FOUR
83: PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 )
84: * ..
85: * .. Local Scalars ..
86: DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
87: $ SCL, T, TEST, TMP, ZETA1, ZETA2
88: COMPLEX*16 ALPHA, COSINE, SINE
89: * ..
90: * .. Intrinsic Functions ..
91: INTRINSIC ABS, DCONJG, MAX, SQRT
92: * ..
93: * .. External Functions ..
94: DOUBLE PRECISION DLAMCH
95: COMPLEX*16 ZDOTC
96: EXTERNAL DLAMCH, ZDOTC
97: * ..
98: * .. Executable Statements ..
99: *
100: EPS = DLAMCH( 'Epsilon' )
101: ALPHA = ZDOTC( J, X, 1, W, 1 )
102: *
103: ABSALP = ABS( ALPHA )
104: ABSGAM = ABS( GAMMA )
105: ABSEST = ABS( SEST )
106: *
107: IF( JOB.EQ.1 ) THEN
108: *
109: * Estimating largest singular value
110: *
111: * special cases
112: *
113: IF( SEST.EQ.ZERO ) THEN
114: S1 = MAX( ABSGAM, ABSALP )
115: IF( S1.EQ.ZERO ) THEN
116: S = ZERO
117: C = ONE
118: SESTPR = ZERO
119: ELSE
120: S = ALPHA / S1
121: C = GAMMA / S1
122: TMP = SQRT( S*DCONJG( S )+C*DCONJG( C ) )
123: S = S / TMP
124: C = C / TMP
125: SESTPR = S1*TMP
126: END IF
127: RETURN
128: ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
129: S = ONE
130: C = ZERO
131: TMP = MAX( ABSEST, ABSALP )
132: S1 = ABSEST / TMP
133: S2 = ABSALP / TMP
134: SESTPR = TMP*SQRT( S1*S1+S2*S2 )
135: RETURN
136: ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
137: S1 = ABSGAM
138: S2 = ABSEST
139: IF( S1.LE.S2 ) THEN
140: S = ONE
141: C = ZERO
142: SESTPR = S2
143: ELSE
144: S = ZERO
145: C = ONE
146: SESTPR = S1
147: END IF
148: RETURN
149: ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
150: S1 = ABSGAM
151: S2 = ABSALP
152: IF( S1.LE.S2 ) THEN
153: TMP = S1 / S2
154: SCL = SQRT( ONE+TMP*TMP )
155: SESTPR = S2*SCL
156: S = ( ALPHA / S2 ) / SCL
157: C = ( GAMMA / S2 ) / SCL
158: ELSE
159: TMP = S2 / S1
160: SCL = SQRT( ONE+TMP*TMP )
161: SESTPR = S1*SCL
162: S = ( ALPHA / S1 ) / SCL
163: C = ( GAMMA / S1 ) / SCL
164: END IF
165: RETURN
166: ELSE
167: *
168: * normal case
169: *
170: ZETA1 = ABSALP / ABSEST
171: ZETA2 = ABSGAM / ABSEST
172: *
173: B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
174: C = ZETA1*ZETA1
175: IF( B.GT.ZERO ) THEN
176: T = C / ( B+SQRT( B*B+C ) )
177: ELSE
178: T = SQRT( B*B+C ) - B
179: END IF
180: *
181: SINE = -( ALPHA / ABSEST ) / T
182: COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
183: TMP = SQRT( SINE*DCONJG( SINE )+COSINE*DCONJG( COSINE ) )
184: S = SINE / TMP
185: C = COSINE / TMP
186: SESTPR = SQRT( T+ONE )*ABSEST
187: RETURN
188: END IF
189: *
190: ELSE IF( JOB.EQ.2 ) THEN
191: *
192: * Estimating smallest singular value
193: *
194: * special cases
195: *
196: IF( SEST.EQ.ZERO ) THEN
197: SESTPR = ZERO
198: IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
199: SINE = ONE
200: COSINE = ZERO
201: ELSE
202: SINE = -DCONJG( GAMMA )
203: COSINE = DCONJG( ALPHA )
204: END IF
205: S1 = MAX( ABS( SINE ), ABS( COSINE ) )
206: S = SINE / S1
207: C = COSINE / S1
208: TMP = SQRT( S*DCONJG( S )+C*DCONJG( C ) )
209: S = S / TMP
210: C = C / TMP
211: RETURN
212: ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
213: S = ZERO
214: C = ONE
215: SESTPR = ABSGAM
216: RETURN
217: ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
218: S1 = ABSGAM
219: S2 = ABSEST
220: IF( S1.LE.S2 ) THEN
221: S = ZERO
222: C = ONE
223: SESTPR = S1
224: ELSE
225: S = ONE
226: C = ZERO
227: SESTPR = S2
228: END IF
229: RETURN
230: ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
231: S1 = ABSGAM
232: S2 = ABSALP
233: IF( S1.LE.S2 ) THEN
234: TMP = S1 / S2
235: SCL = SQRT( ONE+TMP*TMP )
236: SESTPR = ABSEST*( TMP / SCL )
237: S = -( DCONJG( GAMMA ) / S2 ) / SCL
238: C = ( DCONJG( ALPHA ) / S2 ) / SCL
239: ELSE
240: TMP = S2 / S1
241: SCL = SQRT( ONE+TMP*TMP )
242: SESTPR = ABSEST / SCL
243: S = -( DCONJG( GAMMA ) / S1 ) / SCL
244: C = ( DCONJG( ALPHA ) / S1 ) / SCL
245: END IF
246: RETURN
247: ELSE
248: *
249: * normal case
250: *
251: ZETA1 = ABSALP / ABSEST
252: ZETA2 = ABSGAM / ABSEST
253: *
254: NORMA = MAX( ONE+ZETA1*ZETA1+ZETA1*ZETA2,
255: $ ZETA1*ZETA2+ZETA2*ZETA2 )
256: *
257: * See if root is closer to zero or to ONE
258: *
259: TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
260: IF( TEST.GE.ZERO ) THEN
261: *
262: * root is close to zero, compute directly
263: *
264: B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
265: C = ZETA2*ZETA2
266: T = C / ( B+SQRT( ABS( B*B-C ) ) )
267: SINE = ( ALPHA / ABSEST ) / ( ONE-T )
268: COSINE = -( GAMMA / ABSEST ) / T
269: SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
270: ELSE
271: *
272: * root is closer to ONE, shift by that amount
273: *
274: B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
275: C = ZETA1*ZETA1
276: IF( B.GE.ZERO ) THEN
277: T = -C / ( B+SQRT( B*B+C ) )
278: ELSE
279: T = B - SQRT( B*B+C )
280: END IF
281: SINE = -( ALPHA / ABSEST ) / T
282: COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
283: SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
284: END IF
285: TMP = SQRT( SINE*DCONJG( SINE )+COSINE*DCONJG( COSINE ) )
286: S = SINE / TMP
287: C = COSINE / TMP
288: RETURN
289: *
290: END IF
291: END IF
292: RETURN
293: *
294: * End of ZLAIC1
295: *
296: END