1*> \brief \b 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.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLASQ2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLASQ2( N, Z, INFO )
22*
23*       .. Scalar Arguments ..
24*       INTEGER            INFO, N
25*       ..
26*       .. Array Arguments ..
27*       DOUBLE PRECISION   Z( * )
28*       ..
29*
30*
31*> \par Purpose:
32*  =============
33*>
34*> \verbatim
35*>
36*> DLASQ2 computes all the eigenvalues of the symmetric positive
37*> definite tridiagonal matrix associated with the qd array Z to high
38*> relative accuracy are computed to high relative accuracy, in the
39*> absence of denormalization, underflow and overflow.
40*>
41*> To see the relation of Z to the tridiagonal matrix, let L be a
42*> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
43*> let U be an upper bidiagonal matrix with 1's above and diagonal
44*> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
45*> symmetric tridiagonal to which it is similar.
46*>
47*> Note : DLASQ2 defines a logical variable, IEEE, which is true
48*> on machines which follow ieee-754 floating-point standard in their
49*> handling of infinities and NaNs, and false otherwise. This variable
50*> is passed to DLASQ3.
51*> \endverbatim
52*
53*  Arguments:
54*  ==========
55*
56*> \param[in] N
57*> \verbatim
58*>          N is INTEGER
59*>        The number of rows and columns in the matrix. N >= 0.
60*> \endverbatim
61*>
62*> \param[in,out] Z
63*> \verbatim
64*>          Z is DOUBLE PRECISION array, dimension ( 4*N )
65*>        On entry Z holds the qd array. On exit, entries 1 to N hold
66*>        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
67*>        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
68*>        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
69*>        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
70*>        shifts that failed.
71*> \endverbatim
72*>
73*> \param[out] INFO
74*> \verbatim
75*>          INFO is INTEGER
76*>        = 0: successful exit
77*>        < 0: if the i-th argument is a scalar and had an illegal
78*>             value, then INFO = -i, if the i-th argument is an
79*>             array and the j-entry had an illegal value, then
80*>             INFO = -(i*100+j)
81*>        > 0: the algorithm failed
82*>              = 1, a split was marked by a positive value in E
83*>              = 2, current block of Z not diagonalized after 100*N
84*>                   iterations (in inner while loop).  On exit Z holds
85*>                   a qd array with the same eigenvalues as the given Z.
86*>              = 3, termination criterion of outer while loop not met
87*>                   (program created more than N unreduced blocks)
88*> \endverbatim
89*
90*  Authors:
91*  ========
92*
93*> \author Univ. of Tennessee
94*> \author Univ. of California Berkeley
95*> \author Univ. of Colorado Denver
96*> \author NAG Ltd.
97*
98*> \ingroup auxOTHERcomputational
99*
100*> \par Further Details:
101*  =====================
102*>
103*> \verbatim
104*>
105*>  Local Variables: I0:N0 defines a current unreduced segment of Z.
106*>  The shifts are accumulated in SIGMA. Iteration count is in ITER.
107*>  Ping-pong is controlled by PP (alternates between 0 and 1).
108*> \endverbatim
109*>
110*  =====================================================================
111      SUBROUTINE DLASQ2( N, Z, INFO )
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      DOUBLE PRECISION   Z( * )
122*     ..
123*
124*  =====================================================================
125*
126*     .. Parameters ..
127      DOUBLE PRECISION   CBIAS
128      PARAMETER          ( CBIAS = 1.50D0 )
129      DOUBLE PRECISION   ZERO, HALF, ONE, TWO, FOUR, HUNDRD
130      PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
131     $                     TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 )
132*     ..
133*     .. Local Scalars ..
134      LOGICAL            IEEE
135      INTEGER            I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB,
136     $                   K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT,
137     $                   TTYPE
138      DOUBLE PRECISION   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           DLASQ3, DLASRT, XERBLA
145*     ..
146*     .. External Functions ..
147      INTEGER            ILAENV
148      DOUBLE PRECISION   DLAMCH
149      EXTERNAL           DLAMCH, ILAENV
150*     ..
151*     .. Intrinsic Functions ..
152      INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
153*     ..
154*     .. Executable Statements ..
155*
156*     Test the input arguments.
157*     (in case DLASQ2 is not called by DLASQ1)
158*
159      INFO = 0
160      EPS = DLAMCH( 'Precision' )
161      SAFMIN = DLAMCH( 'Safe minimum' )
162      TOL = EPS*HUNDRD
163      TOL2 = TOL**2
164*
165      IF( N.LT.0 ) THEN
166         INFO = -1
167         CALL XERBLA( 'DLASQ2', 1 )
168         RETURN
169      ELSE IF( N.EQ.0 ) THEN
170         RETURN
171      ELSE IF( N.EQ.1 ) THEN
172*
173*        1-by-1 case.
174*
175         IF( Z( 1 ).LT.ZERO ) THEN
176            INFO = -201
177            CALL XERBLA( 'DLASQ2', 2 )
178         END IF
179         RETURN
180      ELSE IF( N.EQ.2 ) THEN
181*
182*        2-by-2 case.
183*
184         IF( Z( 1 ).LT.ZERO ) THEN
185            INFO = -201
186            CALL XERBLA( 'DLASQ2', 2 )
187            RETURN
188         ELSE IF( Z( 2 ).LT.ZERO ) THEN
189            INFO = -202
190            CALL XERBLA( 'DLASQ2', 2 )
191            RETURN
192         ELSE IF( Z( 3 ).LT.ZERO ) THEN
193           INFO = -203
194           CALL XERBLA( 'DLASQ2', 2 )
195           RETURN
196         ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
197            D = Z( 3 )
198            Z( 3 ) = Z( 1 )
199            Z( 1 ) = D
200         END IF
201         Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
202         IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
203            T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) )
204            S = Z( 3 )*( Z( 2 ) / T )
205            IF( S.LE.T ) THEN
206               S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
207            ELSE
208               S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
209            END IF
210            T = Z( 1 ) + ( S+Z( 2 ) )
211            Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
212            Z( 1 ) = T
213         END IF
214         Z( 2 ) = Z( 3 )
215         Z( 6 ) = Z( 2 ) + Z( 1 )
216         RETURN
217      END IF
218*
219*     Check for negative data and compute sums of q's and e's.
220*
221      Z( 2*N ) = ZERO
222      EMIN = Z( 2 )
223      QMAX = ZERO
224      ZMAX = ZERO
225      D = ZERO
226      E = ZERO
227*
228      DO 10 K = 1, 2*( N-1 ), 2
229         IF( Z( K ).LT.ZERO ) THEN
230            INFO = -( 200+K )
231            CALL XERBLA( 'DLASQ2', 2 )
232            RETURN
233         ELSE IF( Z( K+1 ).LT.ZERO ) THEN
234            INFO = -( 200+K+1 )
235            CALL XERBLA( 'DLASQ2', 2 )
236            RETURN
237         END IF
238         D = D + Z( K )
239         E = E + Z( K+1 )
240         QMAX = MAX( QMAX, Z( K ) )
241         EMIN = MIN( EMIN, Z( K+1 ) )
242         ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
243   10 CONTINUE
244      IF( Z( 2*N-1 ).LT.ZERO ) THEN
245         INFO = -( 200+2*N-1 )
246         CALL XERBLA( 'DLASQ2', 2 )
247         RETURN
248      END IF
249      D = D + Z( 2*N-1 )
250      QMAX = MAX( QMAX, Z( 2*N-1 ) )
251      ZMAX = MAX( QMAX, ZMAX )
252*
253*     Check for diagonality.
254*
255      IF( E.EQ.ZERO ) THEN
256         DO 20 K = 2, N
257            Z( K ) = Z( 2*K-1 )
258   20    CONTINUE
259         CALL DLASRT( 'D', N, Z, IINFO )
260         Z( 2*N-1 ) = D
261         RETURN
262      END IF
263*
264      TRACE = D + E
265*
266*     Check for zero data.
267*
268      IF( TRACE.EQ.ZERO ) THEN
269         Z( 2*N-1 ) = ZERO
270         RETURN
271      END IF
272*
273*     Check whether the machine is IEEE conformable.
274*
275      IEEE = ( ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 )
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 DLASQ3.
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 DLASQ3
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 DLASQ3( 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            SIGMA = -Z(4*N1-1)
529            GO TO 145
530         END IF
531
532         DO K = 1, N
533            Z( 2*K-1 ) = Z( 4*K-3 )
534*
535*        Only the block 1..N0 is unfinished.  The rest of the e's
536*        must be essentially zero, although sometimes other data
537*        has been stored in them.
538*
539            IF( K.LT.N0 ) THEN
540               Z( 2*K ) = Z( 4*K-1 )
541            ELSE
542               Z( 2*K ) = 0
543            END IF
544         END DO
545         RETURN
546*
547*        end IWHILB
548*
549  150    CONTINUE
550*
551  160 CONTINUE
552*
553      INFO = 3
554      RETURN
555*
556*     end IWHILA
557*
558  170 CONTINUE
559*
560*     Move q's to the front.
561*
562      DO 180 K = 2, N
563         Z( K ) = Z( 4*K-3 )
564  180 CONTINUE
565*
566*     Sort and compute sum of eigenvalues.
567*
568      CALL DLASRT( 'D', N, Z, IINFO )
569*
570      E = ZERO
571      DO 190 K = N, 1, -1
572         E = E + Z( K )
573  190 CONTINUE
574*
575*     Store trace, sum(eigenvalues) and information on performance.
576*
577      Z( 2*N+1 ) = TRACE
578      Z( 2*N+2 ) = E
579      Z( 2*N+3 ) = DBLE( ITER )
580      Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 )
581      Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER )
582      RETURN
583*
584*     End of DLASQ2
585*
586      END
587