1*> \brief \b DLAED8 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLAED8 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed8.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed8.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed8.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
22*                          CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
23*                          GIVCOL, GIVNUM, INDXP, INDX, INFO )
24*
25*       .. Scalar Arguments ..
26*       INTEGER            CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
27*      $                   QSIZ
28*       DOUBLE PRECISION   RHO
29*       ..
30*       .. Array Arguments ..
31*       INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
32*      $                   INDXQ( * ), PERM( * )
33*       DOUBLE PRECISION   D( * ), DLAMDA( * ), GIVNUM( 2, * ),
34*      $                   Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
35*       ..
36*
37*
38*> \par Purpose:
39*  =============
40*>
41*> \verbatim
42*>
43*> DLAED8 merges the two sets of eigenvalues together into a single
44*> sorted set.  Then it tries to deflate the size of the problem.
45*> There are two ways in which deflation can occur:  when two or more
46*> eigenvalues are close together or if there is a tiny element in the
47*> Z vector.  For each such occurrence the order of the related secular
48*> equation problem is reduced by one.
49*> \endverbatim
50*
51*  Arguments:
52*  ==========
53*
54*> \param[in] ICOMPQ
55*> \verbatim
56*>          ICOMPQ is INTEGER
57*>          = 0:  Compute eigenvalues only.
58*>          = 1:  Compute eigenvectors of original dense symmetric matrix
59*>                also.  On entry, Q contains the orthogonal matrix used
60*>                to reduce the original matrix to tridiagonal form.
61*> \endverbatim
62*>
63*> \param[out] K
64*> \verbatim
65*>          K is INTEGER
66*>         The number of non-deflated eigenvalues, and the order of the
67*>         related secular equation.
68*> \endverbatim
69*>
70*> \param[in] N
71*> \verbatim
72*>          N is INTEGER
73*>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
74*> \endverbatim
75*>
76*> \param[in] QSIZ
77*> \verbatim
78*>          QSIZ is INTEGER
79*>         The dimension of the orthogonal matrix used to reduce
80*>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
81*> \endverbatim
82*>
83*> \param[in,out] D
84*> \verbatim
85*>          D is DOUBLE PRECISION array, dimension (N)
86*>         On entry, the eigenvalues of the two submatrices to be
87*>         combined.  On exit, the trailing (N-K) updated eigenvalues
88*>         (those which were deflated) sorted into increasing order.
89*> \endverbatim
90*>
91*> \param[in,out] Q
92*> \verbatim
93*>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
94*>         If ICOMPQ = 0, Q is not referenced.  Otherwise,
95*>         on entry, Q contains the eigenvectors of the partially solved
96*>         system which has been previously updated in matrix
97*>         multiplies with other partially solved eigensystems.
98*>         On exit, Q contains the trailing (N-K) updated eigenvectors
99*>         (those which were deflated) in its last N-K columns.
100*> \endverbatim
101*>
102*> \param[in] LDQ
103*> \verbatim
104*>          LDQ is INTEGER
105*>         The leading dimension of the array Q.  LDQ >= max(1,N).
106*> \endverbatim
107*>
108*> \param[in] INDXQ
109*> \verbatim
110*>          INDXQ is INTEGER array, dimension (N)
111*>         The permutation which separately sorts the two sub-problems
112*>         in D into ascending order.  Note that elements in the second
113*>         half of this permutation must first have CUTPNT added to
114*>         their values in order to be accurate.
115*> \endverbatim
116*>
117*> \param[in,out] RHO
118*> \verbatim
119*>          RHO is DOUBLE PRECISION
120*>         On entry, the off-diagonal element associated with the rank-1
121*>         cut which originally split the two submatrices which are now
122*>         being recombined.
123*>         On exit, RHO has been modified to the value required by
124*>         DLAED3.
125*> \endverbatim
126*>
127*> \param[in] CUTPNT
128*> \verbatim
129*>          CUTPNT is INTEGER
130*>         The location of the last eigenvalue in the leading
131*>         sub-matrix.  min(1,N) <= CUTPNT <= N.
132*> \endverbatim
133*>
134*> \param[in] Z
135*> \verbatim
136*>          Z is DOUBLE PRECISION array, dimension (N)
137*>         On entry, Z contains the updating vector (the last row of
138*>         the first sub-eigenvector matrix and the first row of the
139*>         second sub-eigenvector matrix).
140*>         On exit, the contents of Z are destroyed by the updating
141*>         process.
142*> \endverbatim
143*>
144*> \param[out] DLAMDA
145*> \verbatim
146*>          DLAMDA is DOUBLE PRECISION array, dimension (N)
147*>         A copy of the first K eigenvalues which will be used by
148*>         DLAED3 to form the secular equation.
149*> \endverbatim
150*>
151*> \param[out] Q2
152*> \verbatim
153*>          Q2 is DOUBLE PRECISION array, dimension (LDQ2,N)
154*>         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
155*>         a copy of the first K eigenvectors which will be used by
156*>         DLAED7 in a matrix multiply (DGEMM) to update the new
157*>         eigenvectors.
158*> \endverbatim
159*>
160*> \param[in] LDQ2
161*> \verbatim
162*>          LDQ2 is INTEGER
163*>         The leading dimension of the array Q2.  LDQ2 >= max(1,N).
164*> \endverbatim
165*>
166*> \param[out] W
167*> \verbatim
168*>          W is DOUBLE PRECISION array, dimension (N)
169*>         The first k values of the final deflation-altered z-vector and
170*>         will be passed to DLAED3.
171*> \endverbatim
172*>
173*> \param[out] PERM
174*> \verbatim
175*>          PERM is INTEGER array, dimension (N)
176*>         The permutations (from deflation and sorting) to be applied
177*>         to each eigenblock.
178*> \endverbatim
179*>
180*> \param[out] GIVPTR
181*> \verbatim
182*>          GIVPTR is INTEGER
183*>         The number of Givens rotations which took place in this
184*>         subproblem.
185*> \endverbatim
186*>
187*> \param[out] GIVCOL
188*> \verbatim
189*>          GIVCOL is INTEGER array, dimension (2, N)
190*>         Each pair of numbers indicates a pair of columns to take place
191*>         in a Givens rotation.
192*> \endverbatim
193*>
194*> \param[out] GIVNUM
195*> \verbatim
196*>          GIVNUM is DOUBLE PRECISION array, dimension (2, N)
197*>         Each number indicates the S value to be used in the
198*>         corresponding Givens rotation.
199*> \endverbatim
200*>
201*> \param[out] INDXP
202*> \verbatim
203*>          INDXP is INTEGER array, dimension (N)
204*>         The permutation used to place deflated values of D at the end
205*>         of the array.  INDXP(1:K) points to the nondeflated D-values
206*>         and INDXP(K+1:N) points to the deflated eigenvalues.
207*> \endverbatim
208*>
209*> \param[out] INDX
210*> \verbatim
211*>          INDX is INTEGER array, dimension (N)
212*>         The permutation used to sort the contents of D into ascending
213*>         order.
214*> \endverbatim
215*>
216*> \param[out] INFO
217*> \verbatim
218*>          INFO is INTEGER
219*>          = 0:  successful exit.
220*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
221*> \endverbatim
222*
223*  Authors:
224*  ========
225*
226*> \author Univ. of Tennessee
227*> \author Univ. of California Berkeley
228*> \author Univ. of Colorado Denver
229*> \author NAG Ltd.
230*
231*> \ingroup auxOTHERcomputational
232*
233*> \par Contributors:
234*  ==================
235*>
236*> Jeff Rutter, Computer Science Division, University of California
237*> at Berkeley, USA
238*
239*  =====================================================================
240      SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
241     $                   CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
242     $                   GIVCOL, GIVNUM, INDXP, INDX, INFO )
243*
244*  -- LAPACK computational routine --
245*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
246*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
247*
248*     .. Scalar Arguments ..
249      INTEGER            CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
250     $                   QSIZ
251      DOUBLE PRECISION   RHO
252*     ..
253*     .. Array Arguments ..
254      INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
255     $                   INDXQ( * ), PERM( * )
256      DOUBLE PRECISION   D( * ), DLAMDA( * ), GIVNUM( 2, * ),
257     $                   Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
258*     ..
259*
260*  =====================================================================
261*
262*     .. Parameters ..
263      DOUBLE PRECISION   MONE, ZERO, ONE, TWO, EIGHT
264      PARAMETER          ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
265     $                   TWO = 2.0D0, EIGHT = 8.0D0 )
266*     ..
267*     .. Local Scalars ..
268*
269      INTEGER            I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
270      DOUBLE PRECISION   C, EPS, S, T, TAU, TOL
271*     ..
272*     .. External Functions ..
273      INTEGER            IDAMAX
274      DOUBLE PRECISION   DLAMCH, DLAPY2
275      EXTERNAL           IDAMAX, DLAMCH, DLAPY2
276*     ..
277*     .. External Subroutines ..
278      EXTERNAL           DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA
279*     ..
280*     .. Intrinsic Functions ..
281      INTRINSIC          ABS, MAX, MIN, SQRT
282*     ..
283*     .. Executable Statements ..
284*
285*     Test the input parameters.
286*
287      INFO = 0
288*
289      IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
290         INFO = -1
291      ELSE IF( N.LT.0 ) THEN
292         INFO = -3
293      ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
294         INFO = -4
295      ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
296         INFO = -7
297      ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
298         INFO = -10
299      ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
300         INFO = -14
301      END IF
302      IF( INFO.NE.0 ) THEN
303         CALL XERBLA( 'DLAED8', -INFO )
304         RETURN
305      END IF
306*
307*     Need to initialize GIVPTR to O here in case of quick exit
308*     to prevent an unspecified code behavior (usually sigfault)
309*     when IWORK array on entry to *stedc is not zeroed
310*     (or at least some IWORK entries which used in *laed7 for GIVPTR).
311*
312      GIVPTR = 0
313*
314*     Quick return if possible
315*
316      IF( N.EQ.0 )
317     $   RETURN
318*
319      N1 = CUTPNT
320      N2 = N - N1
321      N1P1 = N1 + 1
322*
323      IF( RHO.LT.ZERO ) THEN
324         CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
325      END IF
326*
327*     Normalize z so that norm(z) = 1
328*
329      T = ONE / SQRT( TWO )
330      DO 10 J = 1, N
331         INDX( J ) = J
332   10 CONTINUE
333      CALL DSCAL( N, T, Z, 1 )
334      RHO = ABS( TWO*RHO )
335*
336*     Sort the eigenvalues into increasing order
337*
338      DO 20 I = CUTPNT + 1, N
339         INDXQ( I ) = INDXQ( I ) + CUTPNT
340   20 CONTINUE
341      DO 30 I = 1, N
342         DLAMDA( I ) = D( INDXQ( I ) )
343         W( I ) = Z( INDXQ( I ) )
344   30 CONTINUE
345      I = 1
346      J = CUTPNT + 1
347      CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
348      DO 40 I = 1, N
349         D( I ) = DLAMDA( INDX( I ) )
350         Z( I ) = W( INDX( I ) )
351   40 CONTINUE
352*
353*     Calculate the allowable deflation tolerance
354*
355      IMAX = IDAMAX( N, Z, 1 )
356      JMAX = IDAMAX( N, D, 1 )
357      EPS = DLAMCH( 'Epsilon' )
358      TOL = EIGHT*EPS*ABS( D( JMAX ) )
359*
360*     If the rank-1 modifier is small enough, no more needs to be done
361*     except to reorganize Q so that its columns correspond with the
362*     elements in D.
363*
364      IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
365         K = 0
366         IF( ICOMPQ.EQ.0 ) THEN
367            DO 50 J = 1, N
368               PERM( J ) = INDXQ( INDX( J ) )
369   50       CONTINUE
370         ELSE
371            DO 60 J = 1, N
372               PERM( J ) = INDXQ( INDX( J ) )
373               CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
374   60       CONTINUE
375            CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ),
376     $                   LDQ )
377         END IF
378         RETURN
379      END IF
380*
381*     If there are multiple eigenvalues then the problem deflates.  Here
382*     the number of equal eigenvalues are found.  As each equal
383*     eigenvalue is found, an elementary reflector is computed to rotate
384*     the corresponding eigensubspace so that the corresponding
385*     components of Z are zero in this new basis.
386*
387      K = 0
388      K2 = N + 1
389      DO 70 J = 1, N
390         IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
391*
392*           Deflate due to small z component.
393*
394            K2 = K2 - 1
395            INDXP( K2 ) = J
396            IF( J.EQ.N )
397     $         GO TO 110
398         ELSE
399            JLAM = J
400            GO TO 80
401         END IF
402   70 CONTINUE
403   80 CONTINUE
404      J = J + 1
405      IF( J.GT.N )
406     $   GO TO 100
407      IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
408*
409*        Deflate due to small z component.
410*
411         K2 = K2 - 1
412         INDXP( K2 ) = J
413      ELSE
414*
415*        Check if eigenvalues are close enough to allow deflation.
416*
417         S = Z( JLAM )
418         C = Z( J )
419*
420*        Find sqrt(a**2+b**2) without overflow or
421*        destructive underflow.
422*
423         TAU = DLAPY2( C, S )
424         T = D( J ) - D( JLAM )
425         C = C / TAU
426         S = -S / TAU
427         IF( ABS( T*C*S ).LE.TOL ) THEN
428*
429*           Deflation is possible.
430*
431            Z( J ) = TAU
432            Z( JLAM ) = ZERO
433*
434*           Record the appropriate Givens rotation
435*
436            GIVPTR = GIVPTR + 1
437            GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
438            GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
439            GIVNUM( 1, GIVPTR ) = C
440            GIVNUM( 2, GIVPTR ) = S
441            IF( ICOMPQ.EQ.1 ) THEN
442               CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
443     $                    Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
444            END IF
445            T = D( JLAM )*C*C + D( J )*S*S
446            D( J ) = D( JLAM )*S*S + D( J )*C*C
447            D( JLAM ) = T
448            K2 = K2 - 1
449            I = 1
450   90       CONTINUE
451            IF( K2+I.LE.N ) THEN
452               IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
453                  INDXP( K2+I-1 ) = INDXP( K2+I )
454                  INDXP( K2+I ) = JLAM
455                  I = I + 1
456                  GO TO 90
457               ELSE
458                  INDXP( K2+I-1 ) = JLAM
459               END IF
460            ELSE
461               INDXP( K2+I-1 ) = JLAM
462            END IF
463            JLAM = J
464         ELSE
465            K = K + 1
466            W( K ) = Z( JLAM )
467            DLAMDA( K ) = D( JLAM )
468            INDXP( K ) = JLAM
469            JLAM = J
470         END IF
471      END IF
472      GO TO 80
473  100 CONTINUE
474*
475*     Record the last eigenvalue.
476*
477      K = K + 1
478      W( K ) = Z( JLAM )
479      DLAMDA( K ) = D( JLAM )
480      INDXP( K ) = JLAM
481*
482  110 CONTINUE
483*
484*     Sort the eigenvalues and corresponding eigenvectors into DLAMDA
485*     and Q2 respectively.  The eigenvalues/vectors which were not
486*     deflated go into the first K slots of DLAMDA and Q2 respectively,
487*     while those which were deflated go into the last N - K slots.
488*
489      IF( ICOMPQ.EQ.0 ) THEN
490         DO 120 J = 1, N
491            JP = INDXP( J )
492            DLAMDA( J ) = D( JP )
493            PERM( J ) = INDXQ( INDX( JP ) )
494  120    CONTINUE
495      ELSE
496         DO 130 J = 1, N
497            JP = INDXP( J )
498            DLAMDA( J ) = D( JP )
499            PERM( J ) = INDXQ( INDX( JP ) )
500            CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
501  130    CONTINUE
502      END IF
503*
504*     The deflated eigenvalues and their corresponding vectors go back
505*     into the last N - K slots of D and Q respectively.
506*
507      IF( K.LT.N ) THEN
508         IF( ICOMPQ.EQ.0 ) THEN
509            CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
510         ELSE
511            CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
512            CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2,
513     $                   Q( 1, K+1 ), LDQ )
514         END IF
515      END IF
516*
517      RETURN
518*
519*     End of DLAED8
520*
521      END
522