LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine slasq2 ( integer  N,
real, dimension( * )  Z,
integer  INFO 
)

SLASQ2 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 SLASQ2 + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 SLASQ2 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 : SLASQ2 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 SLASQ3.
Parameters
[in]N
          N is INTEGER
        The number of rows and columns in the matrix. N >= 0.
[in,out]Z
          Z is REAL 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 slasq2.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  REAL z( * )
125 * ..
126 *
127 * =====================================================================
128 *
129 * .. Parameters ..
130  REAL cbias
131  parameter ( cbias = 1.50e0 )
132  REAL zero, half, one, two, four, hundrd
133  parameter ( zero = 0.0e0, half = 0.5e0, one = 1.0e0,
134  $ two = 2.0e0, four = 4.0e0, hundrd = 100.0e0 )
135 * ..
136 * .. Local Scalars ..
137  LOGICAL ieee
138  INTEGER i0, i4, iinfo, ipn4, iter, iwhila, iwhilb, k,
139  $ kmin, n0, nbig, ndiv, nfail, pp, splt, ttype,
140  $ i1, n1
141  REAL 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 slasq3, slasrt, xerbla
148 * ..
149 * .. External Functions ..
150  INTEGER ilaenv
151  REAL slamch
152  EXTERNAL ilaenv, slamch
153 * ..
154 * .. Intrinsic Functions ..
155  INTRINSIC abs, max, min, REAL, sqrt
156 * ..
157 * .. Executable Statements ..
158 *
159 * Test the input arguments.
160 * (in case SLASQ2 is not called by SLASQ1)
161 *
162  info = 0
163  eps = slamch( 'Precision' )
164  safmin = slamch( 'Safe minimum' )
165  tol = eps*hundrd
166  tol2 = tol**2
167 *
168  IF( n.LT.0 ) THEN
169  info = -1
170  CALL xerbla( 'SLASQ2', 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( 'SLASQ2', 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( 'SLASQ2', 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( 'SLASQ2', 2 )
227  RETURN
228  ELSE IF( z( k+1 ).LT.zero ) THEN
229  info = -( 200+k+1 )
230  CALL xerbla( 'SLASQ2', 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( 'SLASQ2', 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 slasrt( '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, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
271 * $ ILAENV( 11, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1
272 *
273 * [11/15/2008] The case IEEE=.TRUE. has a problem in single precision with
274 * some the test matrices of type 16. The double precision code is fine.
275 *
276  ieee = .false.
277 *
278 * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
279 *
280  DO 30 k = 2*n, 2, -2
281  z( 2*k ) = zero
282  z( 2*k-1 ) = z( k )
283  z( 2*k-2 ) = zero
284  z( 2*k-3 ) = z( k-1 )
285  30 CONTINUE
286 *
287  i0 = 1
288  n0 = n
289 *
290 * Reverse the qd-array, if warranted.
291 *
292  IF( cbias*z( 4*i0-3 ).LT.z( 4*n0-3 ) ) THEN
293  ipn4 = 4*( i0+n0 )
294  DO 40 i4 = 4*i0, 2*( i0+n0-1 ), 4
295  temp = z( i4-3 )
296  z( i4-3 ) = z( ipn4-i4-3 )
297  z( ipn4-i4-3 ) = temp
298  temp = z( i4-1 )
299  z( i4-1 ) = z( ipn4-i4-5 )
300  z( ipn4-i4-5 ) = temp
301  40 CONTINUE
302  END IF
303 *
304 * Initial split checking via dqd and Li's test.
305 *
306  pp = 0
307 *
308  DO 80 k = 1, 2
309 *
310  d = z( 4*n0+pp-3 )
311  DO 50 i4 = 4*( n0-1 ) + pp, 4*i0 + pp, -4
312  IF( z( i4-1 ).LE.tol2*d ) THEN
313  z( i4-1 ) = -zero
314  d = z( i4-3 )
315  ELSE
316  d = z( i4-3 )*( d / ( d+z( i4-1 ) ) )
317  END IF
318  50 CONTINUE
319 *
320 * dqd maps Z to ZZ plus Li's test.
321 *
322  emin = z( 4*i0+pp+1 )
323  d = z( 4*i0+pp-3 )
324  DO 60 i4 = 4*i0 + pp, 4*( n0-1 ) + pp, 4
325  z( i4-2*pp-2 ) = d + z( i4-1 )
326  IF( z( i4-1 ).LE.tol2*d ) THEN
327  z( i4-1 ) = -zero
328  z( i4-2*pp-2 ) = d
329  z( i4-2*pp ) = zero
330  d = z( i4+1 )
331  ELSE IF( safmin*z( i4+1 ).LT.z( i4-2*pp-2 ) .AND.
332  $ safmin*z( i4-2*pp-2 ).LT.z( i4+1 ) ) THEN
333  temp = z( i4+1 ) / z( i4-2*pp-2 )
334  z( i4-2*pp ) = z( i4-1 )*temp
335  d = d*temp
336  ELSE
337  z( i4-2*pp ) = z( i4+1 )*( z( i4-1 ) / z( i4-2*pp-2 ) )
338  d = z( i4+1 )*( d / z( i4-2*pp-2 ) )
339  END IF
340  emin = min( emin, z( i4-2*pp ) )
341  60 CONTINUE
342  z( 4*n0-pp-2 ) = d
343 *
344 * Now find qmax.
345 *
346  qmax = z( 4*i0-pp-2 )
347  DO 70 i4 = 4*i0 - pp + 2, 4*n0 - pp - 2, 4
348  qmax = max( qmax, z( i4 ) )
349  70 CONTINUE
350 *
351 * Prepare for the next iteration on K.
352 *
353  pp = 1 - pp
354  80 CONTINUE
355 *
356 * Initialise variables to pass to SLASQ3.
357 *
358  ttype = 0
359  dmin1 = zero
360  dmin2 = zero
361  dn = zero
362  dn1 = zero
363  dn2 = zero
364  g = zero
365  tau = zero
366 *
367  iter = 2
368  nfail = 0
369  ndiv = 2*( n0-i0 )
370 *
371  DO 160 iwhila = 1, n + 1
372  IF( n0.LT.1 )
373  $ GO TO 170
374 *
375 * While array unfinished do
376 *
377 * E(N0) holds the value of SIGMA when submatrix in I0:N0
378 * splits from the rest of the array, but is negated.
379 *
380  desig = zero
381  IF( n0.EQ.n ) THEN
382  sigma = zero
383  ELSE
384  sigma = -z( 4*n0-1 )
385  END IF
386  IF( sigma.LT.zero ) THEN
387  info = 1
388  RETURN
389  END IF
390 *
391 * Find last unreduced submatrix's top index I0, find QMAX and
392 * EMIN. Find Gershgorin-type bound if Q's much greater than E's.
393 *
394  emax = zero
395  IF( n0.GT.i0 ) THEN
396  emin = abs( z( 4*n0-5 ) )
397  ELSE
398  emin = zero
399  END IF
400  qmin = z( 4*n0-3 )
401  qmax = qmin
402  DO 90 i4 = 4*n0, 8, -4
403  IF( z( i4-5 ).LE.zero )
404  $ GO TO 100
405  IF( qmin.GE.four*emax ) THEN
406  qmin = min( qmin, z( i4-3 ) )
407  emax = max( emax, z( i4-5 ) )
408  END IF
409  qmax = max( qmax, z( i4-7 )+z( i4-5 ) )
410  emin = min( emin, z( i4-5 ) )
411  90 CONTINUE
412  i4 = 4
413 *
414  100 CONTINUE
415  i0 = i4 / 4
416  pp = 0
417 *
418  IF( n0-i0.GT.1 ) THEN
419  dee = z( 4*i0-3 )
420  deemin = dee
421  kmin = i0
422  DO 110 i4 = 4*i0+1, 4*n0-3, 4
423  dee = z( i4 )*( dee /( dee+z( i4-2 ) ) )
424  IF( dee.LE.deemin ) THEN
425  deemin = dee
426  kmin = ( i4+3 )/4
427  END IF
428  110 CONTINUE
429  IF( (kmin-i0)*2.LT.n0-kmin .AND.
430  $ deemin.LE.half*z(4*n0-3) ) THEN
431  ipn4 = 4*( i0+n0 )
432  pp = 2
433  DO 120 i4 = 4*i0, 2*( i0+n0-1 ), 4
434  temp = z( i4-3 )
435  z( i4-3 ) = z( ipn4-i4-3 )
436  z( ipn4-i4-3 ) = temp
437  temp = z( i4-2 )
438  z( i4-2 ) = z( ipn4-i4-2 )
439  z( ipn4-i4-2 ) = temp
440  temp = z( i4-1 )
441  z( i4-1 ) = z( ipn4-i4-5 )
442  z( ipn4-i4-5 ) = temp
443  temp = z( i4 )
444  z( i4 ) = z( ipn4-i4-4 )
445  z( ipn4-i4-4 ) = temp
446  120 CONTINUE
447  END IF
448  END IF
449 *
450 * Put -(initial shift) into DMIN.
451 *
452  dmin = -max( zero, qmin-two*sqrt( qmin )*sqrt( emax ) )
453 *
454 * Now I0:N0 is unreduced.
455 * PP = 0 for ping, PP = 1 for pong.
456 * PP = 2 indicates that flipping was applied to the Z array and
457 * and that the tests for deflation upon entry in SLASQ3
458 * should not be performed.
459 *
460  nbig = 100*( n0-i0+1 )
461  DO 140 iwhilb = 1, nbig
462  IF( i0.GT.n0 )
463  $ GO TO 150
464 *
465 * While submatrix unfinished take a good dqds step.
466 *
467  CALL slasq3( i0, n0, z, pp, dmin, sigma, desig, qmax, nfail,
468  $ iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1,
469  $ dn2, g, tau )
470 *
471  pp = 1 - pp
472 *
473 * When EMIN is very small check for splits.
474 *
475  IF( pp.EQ.0 .AND. n0-i0.GE.3 ) THEN
476  IF( z( 4*n0 ).LE.tol2*qmax .OR.
477  $ z( 4*n0-1 ).LE.tol2*sigma ) THEN
478  splt = i0 - 1
479  qmax = z( 4*i0-3 )
480  emin = z( 4*i0-1 )
481  oldemn = z( 4*i0 )
482  DO 130 i4 = 4*i0, 4*( n0-3 ), 4
483  IF( z( i4 ).LE.tol2*z( i4-3 ) .OR.
484  $ z( i4-1 ).LE.tol2*sigma ) THEN
485  z( i4-1 ) = -sigma
486  splt = i4 / 4
487  qmax = zero
488  emin = z( i4+3 )
489  oldemn = z( i4+4 )
490  ELSE
491  qmax = max( qmax, z( i4+1 ) )
492  emin = min( emin, z( i4-1 ) )
493  oldemn = min( oldemn, z( i4 ) )
494  END IF
495  130 CONTINUE
496  z( 4*n0-1 ) = emin
497  z( 4*n0 ) = oldemn
498  i0 = splt + 1
499  END IF
500  END IF
501 *
502  140 CONTINUE
503 *
504  info = 2
505 *
506 * Maximum number of iterations exceeded, restore the shift
507 * SIGMA and place the new d's and e's in a qd array.
508 * This might need to be done for several blocks
509 *
510  i1 = i0
511  n1 = n0
512  145 CONTINUE
513  tempq = z( 4*i0-3 )
514  z( 4*i0-3 ) = z( 4*i0-3 ) + sigma
515  DO k = i0+1, n0
516  tempe = z( 4*k-5 )
517  z( 4*k-5 ) = z( 4*k-5 ) * (tempq / z( 4*k-7 ))
518  tempq = z( 4*k-3 )
519  z( 4*k-3 ) = z( 4*k-3 ) + sigma + tempe - z( 4*k-5 )
520  END DO
521 *
522 * Prepare to do this on the previous block if there is one
523 *
524  IF( i1.GT.1 ) THEN
525  n1 = i1-1
526  DO WHILE( ( i1.GE.2 ) .AND. ( z(4*i1-5).GE.zero ) )
527  i1 = i1 - 1
528  END DO
529  IF( i1.GE.1 ) THEN
530  sigma = -z(4*n1-1)
531  GO TO 145
532  END IF
533  END IF
534 
535  DO k = 1, n
536  z( 2*k-1 ) = z( 4*k-3 )
537 *
538 * Only the block 1..N0 is unfinished. The rest of the e's
539 * must be essentially zero, although sometimes other data
540 * has been stored in them.
541 *
542  IF( k.LT.n0 ) THEN
543  z( 2*k ) = z( 4*k-1 )
544  ELSE
545  z( 2*k ) = 0
546  END IF
547  END DO
548  RETURN
549 *
550 * end IWHILB
551 *
552  150 CONTINUE
553 *
554  160 CONTINUE
555 *
556  info = 3
557  RETURN
558 *
559 * end IWHILA
560 *
561  170 CONTINUE
562 *
563 * Move q's to the front.
564 *
565  DO 180 k = 2, n
566  z( k ) = z( 4*k-3 )
567  180 CONTINUE
568 *
569 * Sort and compute sum of eigenvalues.
570 *
571  CALL slasrt( 'D', n, z, iinfo )
572 *
573  e = zero
574  DO 190 k = n, 1, -1
575  e = e + z( k )
576  190 CONTINUE
577 *
578 * Store trace, sum(eigenvalues) and information on performance.
579 *
580  z( 2*n+1 ) = trace
581  z( 2*n+2 ) = e
582  z( 2*n+3 ) = REAL( iter )
583  z( 2*n+4 ) = REAL( NDIV ) / REAL( n**2 )
584  z( 2*n+5 ) = hundrd*nfail / REAL( iter )
585  RETURN
586 *
587 * End of SLASQ2
588 *
slasq3
subroutine slasq3(I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU)
SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
Definition: slasq3.f:184
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
slasrt
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:90
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69

Here is the call graph for this function:

Here is the caller graph for this function: