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