1 SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**H 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]**H 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 = x**H * 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
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**H 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]**H 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 = x**H * 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