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

◆ dlasq4()

subroutine dlasq4 ( integer  i0,
integer  n0,
double precision, dimension( * )  z,
integer  pp,
integer  n0in,
double precision  dmin,
double precision  dmin1,
double precision  dmin2,
double precision  dn,
double precision  dn1,
double precision  dn2,
double precision  tau,
integer  ttype,
double precision  g 
)

DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr.

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

Purpose:
 DLASQ4 computes an approximation TAU to the smallest eigenvalue
 using values of d from the previous transform.
Parameters
[in]I0
          I0 is INTEGER
        First index.
[in]N0
          N0 is INTEGER
        Last index.
[in]Z
          Z is DOUBLE PRECISION array, dimension ( 4*N0 )
        Z holds the qd array.
[in]PP
          PP is INTEGER
        PP=0 for ping, PP=1 for pong.
[in]N0IN
          N0IN is INTEGER
        The value of N0 at start of EIGTEST.
[in]DMIN
          DMIN is DOUBLE PRECISION
        Minimum value of d.
[in]DMIN1
          DMIN1 is DOUBLE PRECISION
        Minimum value of d, excluding D( N0 ).
[in]DMIN2
          DMIN2 is DOUBLE PRECISION
        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
[in]DN
          DN is DOUBLE PRECISION
        d(N)
[in]DN1
          DN1 is DOUBLE PRECISION
        d(N-1)
[in]DN2
          DN2 is DOUBLE PRECISION
        d(N-2)
[out]TAU
          TAU is DOUBLE PRECISION
        This is the shift.
[out]TTYPE
          TTYPE is INTEGER
        Shift type.
[in,out]G
          G is DOUBLE PRECISION
        G is passed as an argument in order to save its value between
        calls to DLASQ4.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  CNST1 = 9/16

Definition at line 149 of file dlasq4.f.

151*
152* -- LAPACK computational routine --
153* -- LAPACK is a software package provided by Univ. of Tennessee, --
154* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155*
156* .. Scalar Arguments ..
157 INTEGER I0, N0, N0IN, PP, TTYPE
158 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
159* ..
160* .. Array Arguments ..
161 DOUBLE PRECISION Z( * )
162* ..
163*
164* =====================================================================
165*
166* .. Parameters ..
167 DOUBLE PRECISION CNST1, CNST2, CNST3
168 parameter( cnst1 = 0.5630d0, cnst2 = 1.010d0,
169 $ cnst3 = 1.050d0 )
170 DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
171 parameter( qurtr = 0.250d0, third = 0.3330d0,
172 $ half = 0.50d0, zero = 0.0d0, one = 1.0d0,
173 $ two = 2.0d0, hundrd = 100.0d0 )
174* ..
175* .. Local Scalars ..
176 INTEGER I4, NN, NP
177 DOUBLE PRECISION A2, B1, B2, GAM, GAP1, GAP2, S
178* ..
179* .. Intrinsic Functions ..
180 INTRINSIC max, min, sqrt
181* ..
182* .. Executable Statements ..
183*
184* A negative DMIN forces the shift to take that absolute value
185* TTYPE records the type of shift.
186*
187 IF( dmin.LE.zero ) THEN
188 tau = -dmin
189 ttype = -1
190 RETURN
191 END IF
192*
193 nn = 4*n0 + pp
194 IF( n0in.EQ.n0 ) THEN
195*
196* No eigenvalues deflated.
197*
198 IF( dmin.EQ.dn .OR. dmin.EQ.dn1 ) THEN
199*
200 b1 = sqrt( z( nn-3 ) )*sqrt( z( nn-5 ) )
201 b2 = sqrt( z( nn-7 ) )*sqrt( z( nn-9 ) )
202 a2 = z( nn-7 ) + z( nn-5 )
203*
204* Cases 2 and 3.
205*
206 IF( dmin.EQ.dn .AND. dmin1.EQ.dn1 ) THEN
207 gap2 = dmin2 - a2 - dmin2*qurtr
208 IF( gap2.GT.zero .AND. gap2.GT.b2 ) THEN
209 gap1 = a2 - dn - ( b2 / gap2 )*b2
210 ELSE
211 gap1 = a2 - dn - ( b1+b2 )
212 END IF
213 IF( gap1.GT.zero .AND. gap1.GT.b1 ) THEN
214 s = max( dn-( b1 / gap1 )*b1, half*dmin )
215 ttype = -2
216 ELSE
217 s = zero
218 IF( dn.GT.b1 )
219 $ s = dn - b1
220 IF( a2.GT.( b1+b2 ) )
221 $ s = min( s, a2-( b1+b2 ) )
222 s = max( s, third*dmin )
223 ttype = -3
224 END IF
225 ELSE
226*
227* Case 4.
228*
229 ttype = -4
230 s = qurtr*dmin
231 IF( dmin.EQ.dn ) THEN
232 gam = dn
233 a2 = zero
234 IF( z( nn-5 ) .GT. z( nn-7 ) )
235 $ RETURN
236 b2 = z( nn-5 ) / z( nn-7 )
237 np = nn - 9
238 ELSE
239 np = nn - 2*pp
240 gam = dn1
241 IF( z( np-4 ) .GT. z( np-2 ) )
242 $ RETURN
243 a2 = z( np-4 ) / z( np-2 )
244 IF( z( nn-9 ) .GT. z( nn-11 ) )
245 $ RETURN
246 b2 = z( nn-9 ) / z( nn-11 )
247 np = nn - 13
248 END IF
249*
250* Approximate contribution to norm squared from I < NN-1.
251*
252 a2 = a2 + b2
253 DO 10 i4 = np, 4*i0 - 1 + pp, -4
254 IF( b2.EQ.zero )
255 $ GO TO 20
256 b1 = b2
257 IF( z( i4 ) .GT. z( i4-2 ) )
258 $ RETURN
259 b2 = b2*( z( i4 ) / z( i4-2 ) )
260 a2 = a2 + b2
261 IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
262 $ GO TO 20
263 10 CONTINUE
264 20 CONTINUE
265 a2 = cnst3*a2
266*
267* Rayleigh quotient residual bound.
268*
269 IF( a2.LT.cnst1 )
270 $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
271 END IF
272 ELSE IF( dmin.EQ.dn2 ) THEN
273*
274* Case 5.
275*
276 ttype = -5
277 s = qurtr*dmin
278*
279* Compute contribution to norm squared from I > NN-2.
280*
281 np = nn - 2*pp
282 b1 = z( np-2 )
283 b2 = z( np-6 )
284 gam = dn2
285 IF( z( np-8 ).GT.b2 .OR. z( np-4 ).GT.b1 )
286 $ RETURN
287 a2 = ( z( np-8 ) / b2 )*( one+z( np-4 ) / b1 )
288*
289* Approximate contribution to norm squared from I < NN-2.
290*
291 IF( n0-i0.GT.2 ) THEN
292 b2 = z( nn-13 ) / z( nn-15 )
293 a2 = a2 + b2
294 DO 30 i4 = nn - 17, 4*i0 - 1 + pp, -4
295 IF( b2.EQ.zero )
296 $ GO TO 40
297 b1 = b2
298 IF( z( i4 ) .GT. z( i4-2 ) )
299 $ RETURN
300 b2 = b2*( z( i4 ) / z( i4-2 ) )
301 a2 = a2 + b2
302 IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
303 $ GO TO 40
304 30 CONTINUE
305 40 CONTINUE
306 a2 = cnst3*a2
307 END IF
308*
309 IF( a2.LT.cnst1 )
310 $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
311 ELSE
312*
313* Case 6, no information to guide us.
314*
315 IF( ttype.EQ.-6 ) THEN
316 g = g + third*( one-g )
317 ELSE IF( ttype.EQ.-18 ) THEN
318 g = qurtr*third
319 ELSE
320 g = qurtr
321 END IF
322 s = g*dmin
323 ttype = -6
324 END IF
325*
326 ELSE IF( n0in.EQ.( n0+1 ) ) THEN
327*
328* One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
329*
330 IF( dmin1.EQ.dn1 .AND. dmin2.EQ.dn2 ) THEN
331*
332* Cases 7 and 8.
333*
334 ttype = -7
335 s = third*dmin1
336 IF( z( nn-5 ).GT.z( nn-7 ) )
337 $ RETURN
338 b1 = z( nn-5 ) / z( nn-7 )
339 b2 = b1
340 IF( b2.EQ.zero )
341 $ GO TO 60
342 DO 50 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
343 a2 = b1
344 IF( z( i4 ).GT.z( i4-2 ) )
345 $ RETURN
346 b1 = b1*( z( i4 ) / z( i4-2 ) )
347 b2 = b2 + b1
348 IF( hundrd*max( b1, a2 ).LT.b2 )
349 $ GO TO 60
350 50 CONTINUE
351 60 CONTINUE
352 b2 = sqrt( cnst3*b2 )
353 a2 = dmin1 / ( one+b2**2 )
354 gap2 = half*dmin2 - a2
355 IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
356 s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
357 ELSE
358 s = max( s, a2*( one-cnst2*b2 ) )
359 ttype = -8
360 END IF
361 ELSE
362*
363* Case 9.
364*
365 s = qurtr*dmin1
366 IF( dmin1.EQ.dn1 )
367 $ s = half*dmin1
368 ttype = -9
369 END IF
370*
371 ELSE IF( n0in.EQ.( n0+2 ) ) THEN
372*
373* Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
374*
375* Cases 10 and 11.
376*
377 IF( dmin2.EQ.dn2 .AND. two*z( nn-5 ).LT.z( nn-7 ) ) THEN
378 ttype = -10
379 s = third*dmin2
380 IF( z( nn-5 ).GT.z( nn-7 ) )
381 $ RETURN
382 b1 = z( nn-5 ) / z( nn-7 )
383 b2 = b1
384 IF( b2.EQ.zero )
385 $ GO TO 80
386 DO 70 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
387 IF( z( i4 ).GT.z( i4-2 ) )
388 $ RETURN
389 b1 = b1*( z( i4 ) / z( i4-2 ) )
390 b2 = b2 + b1
391 IF( hundrd*b1.LT.b2 )
392 $ GO TO 80
393 70 CONTINUE
394 80 CONTINUE
395 b2 = sqrt( cnst3*b2 )
396 a2 = dmin2 / ( one+b2**2 )
397 gap2 = z( nn-7 ) + z( nn-9 ) -
398 $ sqrt( z( nn-11 ) )*sqrt( z( nn-9 ) ) - a2
399 IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
400 s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
401 ELSE
402 s = max( s, a2*( one-cnst2*b2 ) )
403 END IF
404 ELSE
405 s = qurtr*dmin2
406 ttype = -11
407 END IF
408 ELSE IF( n0in.GT.( n0+2 ) ) THEN
409*
410* Case 12, more than two eigenvalues deflated. No information.
411*
412 s = zero
413 ttype = -12
414 END IF
415*
416 tau = s
417 RETURN
418*
419* End of DLASQ4
420*
Here is the caller graph for this function: