132 SUBROUTINE zlaic1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
140 DOUBLE PRECISION SEST, SESTPR
141 COMPLEX*16 C, GAMMA, S
144 COMPLEX*16 W( J ), X( J )
150 DOUBLE PRECISION ZERO, ONE, TWO
151 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
152 DOUBLE PRECISION HALF, FOUR
153 parameter( half = 0.5d0, four = 4.0d0 )
156 DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
157 $ SCL, T, TEST, TMP, ZETA1, ZETA2
158 COMPLEX*16 ALPHA, COSINE, SINE
161 INTRINSIC abs, dconjg, max, sqrt
164 DOUBLE PRECISION DLAMCH
166 EXTERNAL dlamch, zdotc
170 eps = dlamch(
'Epsilon' )
171 alpha = zdotc( j, x, 1, w, 1 )
173 absalp = abs( alpha )
174 absgam = abs( gamma )
183 IF( sest.EQ.zero )
THEN
184 s1 = max( absgam, absalp )
185 IF( s1.EQ.zero )
THEN
192 tmp = dble( sqrt( s*dconjg( s )+c*dconjg( c ) ) )
198 ELSE IF( absgam.LE.eps*absest )
THEN
201 tmp = max( absest, absalp )
204 sestpr = tmp*sqrt( s1*s1+s2*s2 )
206 ELSE IF( absalp.LE.eps*absest )
THEN
219 ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam )
THEN
224 scl = sqrt( one+tmp*tmp )
226 s = ( alpha / s2 ) / scl
227 c = ( gamma / s2 ) / scl
230 scl = sqrt( one+tmp*tmp )
232 s = ( alpha / s1 ) / scl
233 c = ( gamma / s1 ) / scl
240 zeta1 = absalp / absest
241 zeta2 = absgam / absest
243 b = ( one-zeta1*zeta1-zeta2*zeta2 )*half
246 t = dble( c / ( b+sqrt( b*b+c ) ) )
248 t = dble( sqrt( b*b+c ) - b )
251 sine = -( alpha / absest ) / t
252 cosine = -( gamma / absest ) / ( one+t )
253 tmp = dble( sqrt( sine * dconjg( sine )
254 $ + cosine * dconjg( cosine ) ) )
258 sestpr = sqrt( t+one )*absest
262 ELSE IF( job.EQ.2 )
THEN
268 IF( sest.EQ.zero )
THEN
270 IF( max( absgam, absalp ).EQ.zero )
THEN
274 sine = -dconjg( gamma )
275 cosine = dconjg( alpha )
277 s1 = max( abs( sine ), abs( cosine ) )
280 tmp = dble( sqrt( s*dconjg( s )+c*dconjg( c ) ) )
284 ELSE IF( absgam.LE.eps*absest )
THEN
289 ELSE IF( absalp.LE.eps*absest )
THEN
302 ELSE IF( absest.LE.eps*absalp .OR. absest.LE.eps*absgam )
THEN
307 scl = sqrt( one+tmp*tmp )
308 sestpr = absest*( tmp / scl )
309 s = -( dconjg( gamma ) / s2 ) / scl
310 c = ( dconjg( alpha ) / s2 ) / scl
313 scl = sqrt( one+tmp*tmp )
314 sestpr = absest / scl
315 s = -( dconjg( gamma ) / s1 ) / scl
316 c = ( dconjg( alpha ) / s1 ) / scl
323 zeta1 = absalp / absest
324 zeta2 = absgam / absest
326 norma = max( one+zeta1*zeta1+zeta1*zeta2,
327 $ zeta1*zeta2+zeta2*zeta2 )
331 test = one + two*( zeta1-zeta2 )*( zeta1+zeta2 )
332 IF( test.GE.zero )
THEN
336 b = ( zeta1*zeta1+zeta2*zeta2+one )*half
338 t = dble( c / ( b+sqrt( abs( b*b-c ) ) ) )
339 sine = ( alpha / absest ) / ( one-t )
340 cosine = -( gamma / absest ) / t
341 sestpr = sqrt( t+four*eps*eps*norma )*absest
346 b = ( zeta2*zeta2+zeta1*zeta1-one )*half
349 t = dble( -c / ( b+sqrt( b*b+c ) ) )
351 t = dble( b - sqrt( b*b+c ) )
353 sine = -( alpha / absest ) / t
354 cosine = -( gamma / absest ) / ( one+t )
355 sestpr = sqrt( one+t+four*eps*eps*norma )*absest
357 tmp = dble( sqrt( sine * dconjg( sine )
358 $ + cosine * dconjg( cosine ) ) )
subroutine zlaic1(job, j, x, sest, w, gamma, sestpr, s, c)
ZLAIC1 applies one step of incremental condition estimation.