LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slasq2()

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.
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 111 of file slasq2.f.

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