LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlaic1()

subroutine zlaic1 ( integer  job,
integer  j,
complex*16, dimension( j )  x,
double precision  sest,
complex*16, dimension( j )  w,
complex*16  gamma,
double precision  sestpr,
complex*16  s,
complex*16  c 
)

ZLAIC1 applies one step of incremental condition estimation.

Download ZLAIC1 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZLAIC1 applies one step of incremental condition estimation in
 its simplest version:

 Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
 lower triangular matrix L, such that
          twonorm(L*x) = sest
 Then ZLAIC1 computes sestpr, s, c such that
 the vector
                 [ s*x ]
          xhat = [  c  ]
 is an approximate singular vector of
                 [ L       0  ]
          Lhat = [ w**H gamma ]
 in the sense that
          twonorm(Lhat*xhat) = sestpr.

 Depending on JOB, an estimate for the largest or smallest singular
 value is computed.

 Note that [s c]**H and sestpr**2 is an eigenpair of the system

     diag(sest*sest, 0) + [alpha  gamma] * [ conjg(alpha) ]
                                           [ conjg(gamma) ]

 where  alpha =  x**H * w.
Parameters
[in]JOB
          JOB is INTEGER
          = 1: an estimate for the largest singular value is computed.
          = 2: an estimate for the smallest singular value is computed.
[in]J
          J is INTEGER
          Length of X and W
[in]X
          X is COMPLEX*16 array, dimension (J)
          The j-vector x.
[in]SEST
          SEST is DOUBLE PRECISION
          Estimated singular value of j by j matrix L
[in]W
          W is COMPLEX*16 array, dimension (J)
          The j-vector w.
[in]GAMMA
          GAMMA is COMPLEX*16
          The diagonal element gamma.
[out]SESTPR
          SESTPR is DOUBLE PRECISION
          Estimated singular value of (j+1) by (j+1) matrix Lhat.
[out]S
          S is COMPLEX*16
          Sine needed in forming xhat.
[out]C
          C is COMPLEX*16
          Cosine needed in forming xhat.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 134 of file zlaic1.f.

135*
136* -- LAPACK auxiliary routine --
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
194 tmp = dble( sqrt( s*dconjg( s )+c*dconjg( c ) ) )
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
248 t = dble( c / ( b+sqrt( b*b+c ) ) )
249 ELSE
250 t = dble( sqrt( b*b+c ) - b )
251 END IF
252*
253 sine = -( alpha / absest ) / t
254 cosine = -( gamma / absest ) / ( one+t )
255 tmp = dble( sqrt( sine * dconjg( sine )
256 $ + cosine * dconjg( cosine ) ) )
257
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
282 tmp = dble( sqrt( s*dconjg( s )+c*dconjg( c ) ) )
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
340 t = dble( c / ( b+sqrt( abs( b*b-c ) ) ) )
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
351 t = dble( -c / ( b+sqrt( b*b+c ) ) )
352 ELSE
353 t = dble( b - sqrt( b*b+c ) )
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
359 tmp = dble( sqrt( sine * dconjg( sine )
360 $ + 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*
complex *16 function zdotc(n, zx, incx, zy, incy)
ZDOTC
Definition zdotc.f:83
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
Here is the caller graph for this function: