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