LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
June 2016
Further Details: 
  CNST1 = 9/16

Definition at line 153 of file dlasq4.f.

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

Here is the caller graph for this function: