LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlasq2 ( integer  N,
double precision, dimension( * )  Z,
integer  INFO 
)

DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.

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

Purpose: 
 DLASQ2 computes all the eigenvalues of the symmetric positive 
 definite tridiagonal matrix associated with the qd array Z to high
 relative accuracy are computed to high relative accuracy, in the
 absence of denormalization, underflow and overflow.

 To see the relation of Z to the tridiagonal matrix, let L be a
 unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
 let U be an upper bidiagonal matrix with 1's above and diagonal
 Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
 symmetric tridiagonal to which it is similar.

 Note : DLASQ2 defines a logical variable, IEEE, which is true
 on machines which follow ieee-754 floating-point standard in their
 handling of infinities and NaNs, and false otherwise. This variable
 is passed to DLASQ3.
Parameters
[in]N
          N is INTEGER
        The number of rows and columns in the matrix. N >= 0.
[in,out]Z
          Z is DOUBLE PRECISION array, dimension ( 4*N )
        On entry Z holds the qd array. On exit, entries 1 to N hold
        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
        shifts that failed.
[out]INFO
          INFO is INTEGER
        = 0: successful exit
        < 0: if the i-th argument is a scalar and had an illegal
             value, then INFO = -i, if the i-th argument is an
             array and the j-entry had an illegal value, then
             INFO = -(i*100+j)
        > 0: the algorithm failed
              = 1, a split was marked by a positive value in E
              = 2, current block of Z not diagonalized after 100*N
                   iterations (in inner while loop).  On exit Z holds
                   a qd array with the same eigenvalues as the given Z.
              = 3, termination criterion of outer while loop not met 
                   (program created more than N unreduced blocks)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details: 
  Local Variables: I0:N0 defines a current unreduced segment of Z.
  The shifts are accumulated in SIGMA. Iteration count is in ITER.
  Ping-pong is controlled by PP (alternates between 0 and 1).

Definition at line 114 of file dlasq2.f.

114 *
115 * -- LAPACK computational routine (version 3.4.2) --
116 * -- LAPACK is a software package provided by Univ. of Tennessee, --
117 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118 * September 2012
119 *
120 * .. Scalar Arguments ..
121  INTEGER info, n
122 * ..
123 * .. Array Arguments ..
124  DOUBLE PRECISION z( * )
125 * ..
126 *
127 * =====================================================================
128 *
129 * .. Parameters ..
130  DOUBLE PRECISION cbias
131  parameter ( cbias = 1.50d0 )
132  DOUBLE PRECISION zero, half, one, two, four, hundrd
133  parameter ( zero = 0.0d0, half = 0.5d0, one = 1.0d0,
134  $ two = 2.0d0, four = 4.0d0, hundrd = 100.0d0 )
135 * ..
136 * .. Local Scalars ..
137  LOGICAL ieee
138  INTEGER i0, i1, i4, iinfo, ipn4, iter, iwhila, iwhilb,
139  $ k, kmin, n0, n1, nbig, ndiv, nfail, pp, splt,
140  $ ttype
141  DOUBLE PRECISION d, dee, deemin, desig, dmin, dmin1, dmin2, dn,
142  $ dn1, dn2, e, emax, emin, eps, g, oldemn, qmax,
143  $ qmin, s, safmin, sigma, t, tau, temp, tol,
144  $ tol2, trace, zmax, tempe, tempq
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL dlasq3, dlasrt, xerbla
148 * ..
149 * .. External Functions ..
150  INTEGER ilaenv
151  DOUBLE PRECISION dlamch
152  EXTERNAL dlamch, ilaenv
153 * ..
154 * .. Intrinsic Functions ..
155  INTRINSIC abs, dble, max, min, sqrt
156 * ..
157 * .. Executable Statements ..
158 *
159 * Test the input arguments.
160 * (in case DLASQ2 is not called by DLASQ1)
161 *
162  info = 0
163  eps = dlamch( 'Precision' )
164  safmin = dlamch( 'Safe minimum' )
165  tol = eps*hundrd
166  tol2 = tol**2
167 *
168  IF( n.LT.0 ) THEN
169  info = -1
170  CALL xerbla( 'DLASQ2', 1 )
171  RETURN
172  ELSE IF( n.EQ.0 ) THEN
173  RETURN
174  ELSE IF( n.EQ.1 ) THEN
175 *
176 * 1-by-1 case.
177 *
178  IF( z( 1 ).LT.zero ) THEN
179  info = -201
180  CALL xerbla( 'DLASQ2', 2 )
181  END IF
182  RETURN
183  ELSE IF( n.EQ.2 ) THEN
184 *
185 * 2-by-2 case.
186 *
187  IF( z( 2 ).LT.zero .OR. z( 3 ).LT.zero ) THEN
188  info = -2
189  CALL xerbla( 'DLASQ2', 2 )
190  RETURN
191  ELSE IF( z( 3 ).GT.z( 1 ) ) THEN
192  d = z( 3 )
193  z( 3 ) = z( 1 )
194  z( 1 ) = d
195  END IF
196  z( 5 ) = z( 1 ) + z( 2 ) + z( 3 )
197  IF( z( 2 ).GT.z( 3 )*tol2 ) THEN
198  t = half*( ( z( 1 )-z( 3 ) )+z( 2 ) )
199  s = z( 3 )*( z( 2 ) / t )
200  IF( s.LE.t ) THEN
201  s = z( 3 )*( z( 2 ) / ( t*( one+sqrt( one+s / t ) ) ) )
202  ELSE
203  s = z( 3 )*( z( 2 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
204  END IF
205  t = z( 1 ) + ( s+z( 2 ) )
206  z( 3 ) = z( 3 )*( z( 1 ) / t )
207  z( 1 ) = t
208  END IF
209  z( 2 ) = z( 3 )
210  z( 6 ) = z( 2 ) + z( 1 )
211  RETURN
212  END IF
213 *
214 * Check for negative data and compute sums of q's and e's.
215 *
216  z( 2*n ) = zero
217  emin = z( 2 )
218  qmax = zero
219  zmax = zero
220  d = zero
221  e = zero
222 *
223  DO 10 k = 1, 2*( n-1 ), 2
224  IF( z( k ).LT.zero ) THEN
225  info = -( 200+k )
226  CALL xerbla( 'DLASQ2', 2 )
227  RETURN
228  ELSE IF( z( k+1 ).LT.zero ) THEN
229  info = -( 200+k+1 )
230  CALL xerbla( 'DLASQ2', 2 )
231  RETURN
232  END IF
233  d = d + z( k )
234  e = e + z( k+1 )
235  qmax = max( qmax, z( k ) )
236  emin = min( emin, z( k+1 ) )
237  zmax = max( qmax, zmax, z( k+1 ) )
238  10 CONTINUE
239  IF( z( 2*n-1 ).LT.zero ) THEN
240  info = -( 200+2*n-1 )
241  CALL xerbla( 'DLASQ2', 2 )
242  RETURN
243  END IF
244  d = d + z( 2*n-1 )
245  qmax = max( qmax, z( 2*n-1 ) )
246  zmax = max( qmax, zmax )
247 *
248 * Check for diagonality.
249 *
250  IF( e.EQ.zero ) THEN
251  DO 20 k = 2, n
252  z( k ) = z( 2*k-1 )
253  20 CONTINUE
254  CALL dlasrt( 'D', n, z, iinfo )
255  z( 2*n-1 ) = d
256  RETURN
257  END IF
258 *
259  trace = d + e
260 *
261 * Check for zero data.
262 *
263  IF( trace.EQ.zero ) THEN
264  z( 2*n-1 ) = zero
265  RETURN
266  END IF
267 *
268 * Check whether the machine is IEEE conformable.
269 *
270  ieee = ilaenv( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
271  $ ilaenv( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1
272 *
273 * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
274 *
275  DO 30 k = 2*n, 2, -2
276  z( 2*k ) = zero
277  z( 2*k-1 ) = z( k )
278  z( 2*k-2 ) = zero
279  z( 2*k-3 ) = z( k-1 )
280  30 CONTINUE
281 *
282  i0 = 1
283  n0 = n
284 *
285 * Reverse the qd-array, if warranted.
286 *
287  IF( cbias*z( 4*i0-3 ).LT.z( 4*n0-3 ) ) THEN
288  ipn4 = 4*( i0+n0 )
289  DO 40 i4 = 4*i0, 2*( i0+n0-1 ), 4
290  temp = z( i4-3 )
291  z( i4-3 ) = z( ipn4-i4-3 )
292  z( ipn4-i4-3 ) = temp
293  temp = z( i4-1 )
294  z( i4-1 ) = z( ipn4-i4-5 )
295  z( ipn4-i4-5 ) = temp
296  40 CONTINUE
297  END IF
298 *
299 * Initial split checking via dqd and Li's test.
300 *
301  pp = 0
302 *
303  DO 80 k = 1, 2
304 *
305  d = z( 4*n0+pp-3 )
306  DO 50 i4 = 4*( n0-1 ) + pp, 4*i0 + pp, -4
307  IF( z( i4-1 ).LE.tol2*d ) THEN
308  z( i4-1 ) = -zero
309  d = z( i4-3 )
310  ELSE
311  d = z( i4-3 )*( d / ( d+z( i4-1 ) ) )
312  END IF
313  50 CONTINUE
314 *
315 * dqd maps Z to ZZ plus Li's test.
316 *
317  emin = z( 4*i0+pp+1 )
318  d = z( 4*i0+pp-3 )
319  DO 60 i4 = 4*i0 + pp, 4*( n0-1 ) + pp, 4
320  z( i4-2*pp-2 ) = d + z( i4-1 )
321  IF( z( i4-1 ).LE.tol2*d ) THEN
322  z( i4-1 ) = -zero
323  z( i4-2*pp-2 ) = d
324  z( i4-2*pp ) = zero
325  d = z( i4+1 )
326  ELSE IF( safmin*z( i4+1 ).LT.z( i4-2*pp-2 ) .AND.
327  $ safmin*z( i4-2*pp-2 ).LT.z( i4+1 ) ) THEN
328  temp = z( i4+1 ) / z( i4-2*pp-2 )
329  z( i4-2*pp ) = z( i4-1 )*temp
330  d = d*temp
331  ELSE
332  z( i4-2*pp ) = z( i4+1 )*( z( i4-1 ) / z( i4-2*pp-2 ) )
333  d = z( i4+1 )*( d / z( i4-2*pp-2 ) )
334  END IF
335  emin = min( emin, z( i4-2*pp ) )
336  60 CONTINUE
337  z( 4*n0-pp-2 ) = d
338 *
339 * Now find qmax.
340 *
341  qmax = z( 4*i0-pp-2 )
342  DO 70 i4 = 4*i0 - pp + 2, 4*n0 - pp - 2, 4
343  qmax = max( qmax, z( i4 ) )
344  70 CONTINUE
345 *
346 * Prepare for the next iteration on K.
347 *
348  pp = 1 - pp
349  80 CONTINUE
350 *
351 * Initialise variables to pass to DLASQ3.
352 *
353  ttype = 0
354  dmin1 = zero
355  dmin2 = zero
356  dn = zero
357  dn1 = zero
358  dn2 = zero
359  g = zero
360  tau = zero
361 *
362  iter = 2
363  nfail = 0
364  ndiv = 2*( n0-i0 )
365 *
366  DO 160 iwhila = 1, n + 1
367  IF( n0.LT.1 )
368  $ GO TO 170
369 *
370 * While array unfinished do
371 *
372 * E(N0) holds the value of SIGMA when submatrix in I0:N0
373 * splits from the rest of the array, but is negated.
374 *
375  desig = zero
376  IF( n0.EQ.n ) THEN
377  sigma = zero
378  ELSE
379  sigma = -z( 4*n0-1 )
380  END IF
381  IF( sigma.LT.zero ) THEN
382  info = 1
383  RETURN
384  END IF
385 *
386 * Find last unreduced submatrix's top index I0, find QMAX and
387 * EMIN. Find Gershgorin-type bound if Q's much greater than E's.
388 *
389  emax = zero
390  IF( n0.GT.i0 ) THEN
391  emin = abs( z( 4*n0-5 ) )
392  ELSE
393  emin = zero
394  END IF
395  qmin = z( 4*n0-3 )
396  qmax = qmin
397  DO 90 i4 = 4*n0, 8, -4
398  IF( z( i4-5 ).LE.zero )
399  $ GO TO 100
400  IF( qmin.GE.four*emax ) THEN
401  qmin = min( qmin, z( i4-3 ) )
402  emax = max( emax, z( i4-5 ) )
403  END IF
404  qmax = max( qmax, z( i4-7 )+z( i4-5 ) )
405  emin = min( emin, z( i4-5 ) )
406  90 CONTINUE
407  i4 = 4
408 *
409  100 CONTINUE
410  i0 = i4 / 4
411  pp = 0
412 *
413  IF( n0-i0.GT.1 ) THEN
414  dee = z( 4*i0-3 )
415  deemin = dee
416  kmin = i0
417  DO 110 i4 = 4*i0+1, 4*n0-3, 4
418  dee = z( i4 )*( dee /( dee+z( i4-2 ) ) )
419  IF( dee.LE.deemin ) THEN
420  deemin = dee
421  kmin = ( i4+3 )/4
422  END IF
423  110 CONTINUE
424  IF( (kmin-i0)*2.LT.n0-kmin .AND.
425  $ deemin.LE.half*z(4*n0-3) ) THEN
426  ipn4 = 4*( i0+n0 )
427  pp = 2
428  DO 120 i4 = 4*i0, 2*( i0+n0-1 ), 4
429  temp = z( i4-3 )
430  z( i4-3 ) = z( ipn4-i4-3 )
431  z( ipn4-i4-3 ) = temp
432  temp = z( i4-2 )
433  z( i4-2 ) = z( ipn4-i4-2 )
434  z( ipn4-i4-2 ) = temp
435  temp = z( i4-1 )
436  z( i4-1 ) = z( ipn4-i4-5 )
437  z( ipn4-i4-5 ) = temp
438  temp = z( i4 )
439  z( i4 ) = z( ipn4-i4-4 )
440  z( ipn4-i4-4 ) = temp
441  120 CONTINUE
442  END IF
443  END IF
444 *
445 * Put -(initial shift) into DMIN.
446 *
447  dmin = -max( zero, qmin-two*sqrt( qmin )*sqrt( emax ) )
448 *
449 * Now I0:N0 is unreduced.
450 * PP = 0 for ping, PP = 1 for pong.
451 * PP = 2 indicates that flipping was applied to the Z array and
452 * and that the tests for deflation upon entry in DLASQ3
453 * should not be performed.
454 *
455  nbig = 100*( n0-i0+1 )
456  DO 140 iwhilb = 1, nbig
457  IF( i0.GT.n0 )
458  $ GO TO 150
459 *
460 * While submatrix unfinished take a good dqds step.
461 *
462  CALL dlasq3( i0, n0, z, pp, dmin, sigma, desig, qmax, nfail,
463  $ iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1,
464  $ dn2, g, tau )
465 *
466  pp = 1 - pp
467 *
468 * When EMIN is very small check for splits.
469 *
470  IF( pp.EQ.0 .AND. n0-i0.GE.3 ) THEN
471  IF( z( 4*n0 ).LE.tol2*qmax .OR.
472  $ z( 4*n0-1 ).LE.tol2*sigma ) THEN
473  splt = i0 - 1
474  qmax = z( 4*i0-3 )
475  emin = z( 4*i0-1 )
476  oldemn = z( 4*i0 )
477  DO 130 i4 = 4*i0, 4*( n0-3 ), 4
478  IF( z( i4 ).LE.tol2*z( i4-3 ) .OR.
479  $ z( i4-1 ).LE.tol2*sigma ) THEN
480  z( i4-1 ) = -sigma
481  splt = i4 / 4
482  qmax = zero
483  emin = z( i4+3 )
484  oldemn = z( i4+4 )
485  ELSE
486  qmax = max( qmax, z( i4+1 ) )
487  emin = min( emin, z( i4-1 ) )
488  oldemn = min( oldemn, z( i4 ) )
489  END IF
490  130 CONTINUE
491  z( 4*n0-1 ) = emin
492  z( 4*n0 ) = oldemn
493  i0 = splt + 1
494  END IF
495  END IF
496 *
497  140 CONTINUE
498 *
499  info = 2
500 *
501 * Maximum number of iterations exceeded, restore the shift
502 * SIGMA and place the new d's and e's in a qd array.
503 * This might need to be done for several blocks
504 *
505  i1 = i0
506  n1 = n0
507  145 CONTINUE
508  tempq = z( 4*i0-3 )
509  z( 4*i0-3 ) = z( 4*i0-3 ) + sigma
510  DO k = i0+1, n0
511  tempe = z( 4*k-5 )
512  z( 4*k-5 ) = z( 4*k-5 ) * (tempq / z( 4*k-7 ))
513  tempq = z( 4*k-3 )
514  z( 4*k-3 ) = z( 4*k-3 ) + sigma + tempe - z( 4*k-5 )
515  END DO
516 *
517 * Prepare to do this on the previous block if there is one
518 *
519  IF( i1.GT.1 ) THEN
520  n1 = i1-1
521  DO WHILE( ( i1.GE.2 ) .AND. ( z(4*i1-5).GE.zero ) )
522  i1 = i1 - 1
523  END DO
524  sigma = -z(4*n1-1)
525  GO TO 145
526  END IF
527 
528  DO k = 1, n
529  z( 2*k-1 ) = z( 4*k-3 )
530 *
531 * Only the block 1..N0 is unfinished. The rest of the e's
532 * must be essentially zero, although sometimes other data
533 * has been stored in them.
534 *
535  IF( k.LT.n0 ) THEN
536  z( 2*k ) = z( 4*k-1 )
537  ELSE
538  z( 2*k ) = 0
539  END IF
540  END DO
541  RETURN
542 *
543 * end IWHILB
544 *
545  150 CONTINUE
546 *
547  160 CONTINUE
548 *
549  info = 3
550  RETURN
551 *
552 * end IWHILA
553 *
554  170 CONTINUE
555 *
556 * Move q's to the front.
557 *
558  DO 180 k = 2, n
559  z( k ) = z( 4*k-3 )
560  180 CONTINUE
561 *
562 * Sort and compute sum of eigenvalues.
563 *
564  CALL dlasrt( 'D', n, z, iinfo )
565 *
566  e = zero
567  DO 190 k = n, 1, -1
568  e = e + z( k )
569  190 CONTINUE
570 *
571 * Store trace, sum(eigenvalues) and information on performance.
572 *
573  z( 2*n+1 ) = trace
574  z( 2*n+2 ) = e
575  z( 2*n+3 ) = dble( iter )
576  z( 2*n+4 ) = dble( ndiv ) / dble( n**2 )
577  z( 2*n+5 ) = hundrd*nfail / dble( iter )
578  RETURN
579 *
580 * End of DLASQ2
581 *
dlasrt
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:90
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dlasq3
subroutine dlasq3(I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU)
DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
Definition: dlasq3.f:184
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function: