LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlasq3 ( integer  I0,
integer  N0,
double precision, dimension( * )  Z,
integer  PP,
double precision  DMIN,
double precision  SIGMA,
double precision  DESIG,
double precision  QMAX,
integer  NFAIL,
integer  ITER,
integer  NDIV,
logical  IEEE,
integer  TTYPE,
double precision  DMIN1,
double precision  DMIN2,
double precision  DN,
double precision  DN1,
double precision  DN2,
double precision  G,
double precision  TAU 
)

DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.

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

Purpose: 
 DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
 In case of failure it changes shifts, and tries again until output
 is positive.
Parameters
[in]I0
          I0 is INTEGER
         First index.
[in,out]N0
          N0 is INTEGER
         Last index.
[in,out]Z
          Z is DOUBLE PRECISION array, dimension ( 4*N0 )
         Z holds the qd array.
[in,out]PP
          PP is INTEGER
         PP=0 for ping, PP=1 for pong.
         PP=2 indicates that flipping was applied to the Z array   
         and that the initial tests for deflation should not be 
         performed.
[out]DMIN
          DMIN is DOUBLE PRECISION
         Minimum value of d.
[out]SIGMA
          SIGMA is DOUBLE PRECISION
         Sum of shifts used in current segment.
[in,out]DESIG
          DESIG is DOUBLE PRECISION
         Lower order part of SIGMA
[in]QMAX
          QMAX is DOUBLE PRECISION
         Maximum value of q.
[in,out]NFAIL
          NFAIL is INTEGER
         Increment NFAIL by 1 each time the shift was too big.
[in,out]ITER
          ITER is INTEGER
         Increment ITER by 1 for each iteration.
[in,out]NDIV
          NDIV is INTEGER
         Increment NDIV by 1 for each division.
[in]IEEE
          IEEE is LOGICAL
         Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
[in,out]TTYPE
          TTYPE is INTEGER
         Shift type.
[in,out]DMIN1
          DMIN1 is DOUBLE PRECISION
[in,out]DMIN2
          DMIN2 is DOUBLE PRECISION
[in,out]DN
          DN is DOUBLE PRECISION
[in,out]DN1
          DN1 is DOUBLE PRECISION
[in,out]DN2
          DN2 is DOUBLE PRECISION
[in,out]G
          G is DOUBLE PRECISION
[in,out]TAU
          TAU is DOUBLE PRECISION

         These are passed as arguments in order to save their values
         between calls to DLASQ3.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 184 of file dlasq3.f.

184 *
185 * -- LAPACK computational routine (version 3.6.1) --
186 * -- LAPACK is a software package provided by Univ. of Tennessee, --
187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188 * June 2016
189 *
190 * .. Scalar Arguments ..
191  LOGICAL ieee
192  INTEGER i0, iter, n0, ndiv, nfail, pp
193  DOUBLE PRECISION desig, dmin, dmin1, dmin2, dn, dn1, dn2, g,
194  $ qmax, sigma, tau
195 * ..
196 * .. Array Arguments ..
197  DOUBLE PRECISION z( * )
198 * ..
199 *
200 * =====================================================================
201 *
202 * .. Parameters ..
203  DOUBLE PRECISION cbias
204  parameter ( cbias = 1.50d0 )
205  DOUBLE PRECISION zero, qurtr, half, one, two, hundrd
206  parameter ( zero = 0.0d0, qurtr = 0.250d0, half = 0.5d0,
207  $ one = 1.0d0, two = 2.0d0, hundrd = 100.0d0 )
208 * ..
209 * .. Local Scalars ..
210  INTEGER ipn4, j4, n0in, nn, ttype
211  DOUBLE PRECISION eps, s, t, temp, tol, tol2
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL dlasq4, dlasq5, dlasq6
215 * ..
216 * .. External Function ..
217  DOUBLE PRECISION dlamch
218  LOGICAL disnan
219  EXTERNAL disnan, dlamch
220 * ..
221 * .. Intrinsic Functions ..
222  INTRINSIC abs, max, min, sqrt
223 * ..
224 * .. Executable Statements ..
225 *
226  n0in = n0
227  eps = dlamch( 'Precision' )
228  tol = eps*hundrd
229  tol2 = tol**2
230 *
231 * Check for deflation.
232 *
233  10 CONTINUE
234 *
235  IF( n0.LT.i0 )
236  $ RETURN
237  IF( n0.EQ.i0 )
238  $ GO TO 20
239  nn = 4*n0 + pp
240  IF( n0.EQ.( i0+1 ) )
241  $ GO TO 40
242 *
243 * Check whether E(N0-1) is negligible, 1 eigenvalue.
244 *
245  IF( z( nn-5 ).GT.tol2*( sigma+z( nn-3 ) ) .AND.
246  $ z( nn-2*pp-4 ).GT.tol2*z( nn-7 ) )
247  $ GO TO 30
248 *
249  20 CONTINUE
250 *
251  z( 4*n0-3 ) = z( 4*n0+pp-3 ) + sigma
252  n0 = n0 - 1
253  GO TO 10
254 *
255 * Check whether E(N0-2) is negligible, 2 eigenvalues.
256 *
257  30 CONTINUE
258 *
259  IF( z( nn-9 ).GT.tol2*sigma .AND.
260  $ z( nn-2*pp-8 ).GT.tol2*z( nn-11 ) )
261  $ GO TO 50
262 *
263  40 CONTINUE
264 *
265  IF( z( nn-3 ).GT.z( nn-7 ) ) THEN
266  s = z( nn-3 )
267  z( nn-3 ) = z( nn-7 )
268  z( nn-7 ) = s
269  END IF
270  t = half*( ( z( nn-7 )-z( nn-3 ) )+z( nn-5 ) )
271  IF( z( nn-5 ).GT.z( nn-3 )*tol2.AND.t.NE.zero ) THEN
272  s = z( nn-3 )*( z( nn-5 ) / t )
273  IF( s.LE.t ) THEN
274  s = z( nn-3 )*( z( nn-5 ) /
275  $ ( t*( one+sqrt( one+s / t ) ) ) )
276  ELSE
277  s = z( nn-3 )*( z( nn-5 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
278  END IF
279  t = z( nn-7 ) + ( s+z( nn-5 ) )
280  z( nn-3 ) = z( nn-3 )*( z( nn-7 ) / t )
281  z( nn-7 ) = t
282  END IF
283  z( 4*n0-7 ) = z( nn-7 ) + sigma
284  z( 4*n0-3 ) = z( nn-3 ) + sigma
285  n0 = n0 - 2
286  GO TO 10
287 *
288  50 CONTINUE
289  IF( pp.EQ.2 )
290  $ pp = 0
291 *
292 * Reverse the qd-array, if warranted.
293 *
294  IF( dmin.LE.zero .OR. n0.LT.n0in ) THEN
295  IF( cbias*z( 4*i0+pp-3 ).LT.z( 4*n0+pp-3 ) ) THEN
296  ipn4 = 4*( i0+n0 )
297  DO 60 j4 = 4*i0, 2*( i0+n0-1 ), 4
298  temp = z( j4-3 )
299  z( j4-3 ) = z( ipn4-j4-3 )
300  z( ipn4-j4-3 ) = temp
301  temp = z( j4-2 )
302  z( j4-2 ) = z( ipn4-j4-2 )
303  z( ipn4-j4-2 ) = temp
304  temp = z( j4-1 )
305  z( j4-1 ) = z( ipn4-j4-5 )
306  z( ipn4-j4-5 ) = temp
307  temp = z( j4 )
308  z( j4 ) = z( ipn4-j4-4 )
309  z( ipn4-j4-4 ) = temp
310  60 CONTINUE
311  IF( n0-i0.LE.4 ) THEN
312  z( 4*n0+pp-1 ) = z( 4*i0+pp-1 )
313  z( 4*n0-pp ) = z( 4*i0-pp )
314  END IF
315  dmin2 = min( dmin2, z( 4*n0+pp-1 ) )
316  z( 4*n0+pp-1 ) = min( z( 4*n0+pp-1 ), z( 4*i0+pp-1 ),
317  $ z( 4*i0+pp+3 ) )
318  z( 4*n0-pp ) = min( z( 4*n0-pp ), z( 4*i0-pp ),
319  $ z( 4*i0-pp+4 ) )
320  qmax = max( qmax, z( 4*i0+pp-3 ), z( 4*i0+pp+1 ) )
321  dmin = -zero
322  END IF
323  END IF
324 *
325 * Choose a shift.
326 *
327  CALL dlasq4( i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1,
328  $ dn2, tau, ttype, g )
329 *
330 * Call dqds until DMIN > 0.
331 *
332  70 CONTINUE
333 *
334  CALL dlasq5( i0, n0, z, pp, tau, sigma, dmin, dmin1, dmin2, dn,
335  $ dn1, dn2, ieee, eps )
336 *
337  ndiv = ndiv + ( n0-i0+2 )
338  iter = iter + 1
339 *
340 * Check status.
341 *
342  IF( dmin.GE.zero .AND. dmin1.GE.zero ) THEN
343 *
344 * Success.
345 *
346  GO TO 90
347 *
348  ELSE IF( dmin.LT.zero .AND. dmin1.GT.zero .AND.
349  $ z( 4*( n0-1 )-pp ).LT.tol*( sigma+dn1 ) .AND.
350  $ abs( dn ).LT.tol*sigma ) THEN
351 *
352 * Convergence hidden by negative DN.
353 *
354  z( 4*( n0-1 )-pp+2 ) = zero
355  dmin = zero
356  GO TO 90
357  ELSE IF( dmin.LT.zero ) THEN
358 *
359 * TAU too big. Select new TAU and try again.
360 *
361  nfail = nfail + 1
362  IF( ttype.LT.-22 ) THEN
363 *
364 * Failed twice. Play it safe.
365 *
366  tau = zero
367  ELSE IF( dmin1.GT.zero ) THEN
368 *
369 * Late failure. Gives excellent shift.
370 *
371  tau = ( tau+dmin )*( one-two*eps )
372  ttype = ttype - 11
373  ELSE
374 *
375 * Early failure. Divide by 4.
376 *
377  tau = qurtr*tau
378  ttype = ttype - 12
379  END IF
380  GO TO 70
381  ELSE IF( disnan( dmin ) ) THEN
382 *
383 * NaN.
384 *
385  IF( tau.EQ.zero ) THEN
386  GO TO 80
387  ELSE
388  tau = zero
389  GO TO 70
390  END IF
391  ELSE
392 *
393 * Possible underflow. Play it safe.
394 *
395  GO TO 80
396  END IF
397 *
398 * Risk of underflow.
399 *
400  80 CONTINUE
401  CALL dlasq6( i0, n0, z, pp, dmin, dmin1, dmin2, dn, dn1, dn2 )
402  ndiv = ndiv + ( n0-i0+2 )
403  iter = iter + 1
404  tau = zero
405 *
406  90 CONTINUE
407  IF( tau.LT.sigma ) THEN
408  desig = desig + tau
409  t = sigma + desig
410  desig = desig - ( t-sigma )
411  ELSE
412  t = sigma + tau
413  desig = sigma - ( t-tau ) + desig
414  END IF
415  sigma = t
416 *
417  RETURN
418 *
419 * End of DLASQ3
420 *
disnan
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dlasq4
subroutine dlasq4(I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE, G)
DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous trans...
Definition: dlasq4.f:153
dlasq5
subroutine dlasq5(I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE, EPS)
DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
Definition: dlasq5.f:146
dlasq6
subroutine dlasq6(I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2)
DLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.
Definition: dlasq6.f:121

Here is the call graph for this function:

Here is the caller graph for this function: