LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlasq3()

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.

Definition at line 179 of file dlasq3.f.

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