LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
June 2016
Further Details:
  CNST1 = 9/16

Definition at line 153 of file dlasq4.f.

153 *
154 * -- LAPACK computational routine (version 3.7.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  gam = dn1
244  IF( z( np-4 ) .GT. z( np-2 ) )
245  $ RETURN
246  a2 = z( np-4 ) / z( np-2 )
247  IF( z( nn-9 ) .GT. z( nn-11 ) )
248  $ RETURN
249  b2 = z( nn-9 ) / z( nn-11 )
250  np = nn - 13
251  END IF
252 *
253 * Approximate contribution to norm squared from I < NN-1.
254 *
255  a2 = a2 + b2
256  DO 10 i4 = np, 4*i0 - 1 + pp, -4
257  IF( b2.EQ.zero )
258  $ GO TO 20
259  b1 = b2
260  IF( z( i4 ) .GT. z( i4-2 ) )
261  $ RETURN
262  b2 = b2*( z( i4 ) / z( i4-2 ) )
263  a2 = a2 + b2
264  IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
265  $ GO TO 20
266  10 CONTINUE
267  20 CONTINUE
268  a2 = cnst3*a2
269 *
270 * Rayleigh quotient residual bound.
271 *
272  IF( a2.LT.cnst1 )
273  $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
274  END IF
275  ELSE IF( dmin.EQ.dn2 ) THEN
276 *
277 * Case 5.
278 *
279  ttype = -5
280  s = qurtr*dmin
281 *
282 * Compute contribution to norm squared from I > NN-2.
283 *
284  np = nn - 2*pp
285  b1 = z( np-2 )
286  b2 = z( np-6 )
287  gam = dn2
288  IF( z( np-8 ).GT.b2 .OR. z( np-4 ).GT.b1 )
289  $ RETURN
290  a2 = ( z( np-8 ) / b2 )*( one+z( np-4 ) / b1 )
291 *
292 * Approximate contribution to norm squared from I < NN-2.
293 *
294  IF( n0-i0.GT.2 ) THEN
295  b2 = z( nn-13 ) / z( nn-15 )
296  a2 = a2 + b2
297  DO 30 i4 = nn - 17, 4*i0 - 1 + pp, -4
298  IF( b2.EQ.zero )
299  $ GO TO 40
300  b1 = b2
301  IF( z( i4 ) .GT. z( i4-2 ) )
302  $ RETURN
303  b2 = b2*( z( i4 ) / z( i4-2 ) )
304  a2 = a2 + b2
305  IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
306  $ GO TO 40
307  30 CONTINUE
308  40 CONTINUE
309  a2 = cnst3*a2
310  END IF
311 *
312  IF( a2.LT.cnst1 )
313  $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
314  ELSE
315 *
316 * Case 6, no information to guide us.
317 *
318  IF( ttype.EQ.-6 ) THEN
319  g = g + third*( one-g )
320  ELSE IF( ttype.EQ.-18 ) THEN
321  g = qurtr*third
322  ELSE
323  g = qurtr
324  END IF
325  s = g*dmin
326  ttype = -6
327  END IF
328 *
329  ELSE IF( n0in.EQ.( n0+1 ) ) THEN
330 *
331 * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
332 *
333  IF( dmin1.EQ.dn1 .AND. dmin2.EQ.dn2 ) THEN
334 *
335 * Cases 7 and 8.
336 *
337  ttype = -7
338  s = third*dmin1
339  IF( z( nn-5 ).GT.z( nn-7 ) )
340  $ RETURN
341  b1 = z( nn-5 ) / z( nn-7 )
342  b2 = b1
343  IF( b2.EQ.zero )
344  $ GO TO 60
345  DO 50 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
346  a2 = b1
347  IF( z( i4 ).GT.z( i4-2 ) )
348  $ RETURN
349  b1 = b1*( z( i4 ) / z( i4-2 ) )
350  b2 = b2 + b1
351  IF( hundrd*max( b1, a2 ).LT.b2 )
352  $ GO TO 60
353  50 CONTINUE
354  60 CONTINUE
355  b2 = sqrt( cnst3*b2 )
356  a2 = dmin1 / ( one+b2**2 )
357  gap2 = half*dmin2 - a2
358  IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
359  s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
360  ELSE
361  s = max( s, a2*( one-cnst2*b2 ) )
362  ttype = -8
363  END IF
364  ELSE
365 *
366 * Case 9.
367 *
368  s = qurtr*dmin1
369  IF( dmin1.EQ.dn1 )
370  $ s = half*dmin1
371  ttype = -9
372  END IF
373 *
374  ELSE IF( n0in.EQ.( n0+2 ) ) THEN
375 *
376 * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
377 *
378 * Cases 10 and 11.
379 *
380  IF( dmin2.EQ.dn2 .AND. two*z( nn-5 ).LT.z( nn-7 ) ) THEN
381  ttype = -10
382  s = third*dmin2
383  IF( z( nn-5 ).GT.z( nn-7 ) )
384  $ RETURN
385  b1 = z( nn-5 ) / z( nn-7 )
386  b2 = b1
387  IF( b2.EQ.zero )
388  $ GO TO 80
389  DO 70 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
390  IF( z( i4 ).GT.z( i4-2 ) )
391  $ RETURN
392  b1 = b1*( z( i4 ) / z( i4-2 ) )
393  b2 = b2 + b1
394  IF( hundrd*b1.LT.b2 )
395  $ GO TO 80
396  70 CONTINUE
397  80 CONTINUE
398  b2 = sqrt( cnst3*b2 )
399  a2 = dmin2 / ( one+b2**2 )
400  gap2 = z( nn-7 ) + z( nn-9 ) -
401  $ sqrt( z( nn-11 ) )*sqrt( z( nn-9 ) ) - a2
402  IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
403  s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
404  ELSE
405  s = max( s, a2*( one-cnst2*b2 ) )
406  END IF
407  ELSE
408  s = qurtr*dmin2
409  ttype = -11
410  END IF
411  ELSE IF( n0in.GT.( n0+2 ) ) THEN
412 *
413 * Case 12, more than two eigenvalues deflated. No information.
414 *
415  s = zero
416  ttype = -12
417  END IF
418 *
419  tau = s
420  RETURN
421 *
422 * End of DLASQ4
423 *
Here is the caller graph for this function: