1*> \brief \b DLAED9 used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. 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 DLAED9 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed9.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed9.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed9.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
22*                          S, LDS, INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
26*       DOUBLE PRECISION   RHO
27*       ..
28*       .. Array Arguments ..
29*       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
30*      $                   W( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> DLAED9 finds the roots of the secular equation, as defined by the
40*> values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
41*> appropriate calls to DLAED4 and then stores the new matrix of
42*> eigenvectors for use in calculating the next level of Z vectors.
43*> \endverbatim
44*
45*  Arguments:
46*  ==========
47*
48*> \param[in] K
49*> \verbatim
50*>          K is INTEGER
51*>          The number of terms in the rational function to be solved by
52*>          DLAED4.  K >= 0.
53*> \endverbatim
54*>
55*> \param[in] KSTART
56*> \verbatim
57*>          KSTART is INTEGER
58*> \endverbatim
59*>
60*> \param[in] KSTOP
61*> \verbatim
62*>          KSTOP is INTEGER
63*>          The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
64*>          are to be computed.  1 <= KSTART <= KSTOP <= K.
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*>          N is INTEGER
70*>          The number of rows and columns in the Q matrix.
71*>          N >= K (delation may result in N > K).
72*> \endverbatim
73*>
74*> \param[out] D
75*> \verbatim
76*>          D is DOUBLE PRECISION array, dimension (N)
77*>          D(I) contains the updated eigenvalues
78*>          for KSTART <= I <= KSTOP.
79*> \endverbatim
80*>
81*> \param[out] Q
82*> \verbatim
83*>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
84*> \endverbatim
85*>
86*> \param[in] LDQ
87*> \verbatim
88*>          LDQ is INTEGER
89*>          The leading dimension of the array Q.  LDQ >= max( 1, N ).
90*> \endverbatim
91*>
92*> \param[in] RHO
93*> \verbatim
94*>          RHO is DOUBLE PRECISION
95*>          The value of the parameter in the rank one update equation.
96*>          RHO >= 0 required.
97*> \endverbatim
98*>
99*> \param[in] DLAMDA
100*> \verbatim
101*>          DLAMDA is DOUBLE PRECISION array, dimension (K)
102*>          The first K elements of this array contain the old roots
103*>          of the deflated updating problem.  These are the poles
104*>          of the secular equation.
105*> \endverbatim
106*>
107*> \param[in] W
108*> \verbatim
109*>          W is DOUBLE PRECISION array, dimension (K)
110*>          The first K elements of this array contain the components
111*>          of the deflation-adjusted updating vector.
112*> \endverbatim
113*>
114*> \param[out] S
115*> \verbatim
116*>          S is DOUBLE PRECISION array, dimension (LDS, K)
117*>          Will contain the eigenvectors of the repaired matrix which
118*>          will be stored for subsequent Z vector calculation and
119*>          multiplied by the previously accumulated eigenvectors
120*>          to update the system.
121*> \endverbatim
122*>
123*> \param[in] LDS
124*> \verbatim
125*>          LDS is INTEGER
126*>          The leading dimension of S.  LDS >= max( 1, K ).
127*> \endverbatim
128*>
129*> \param[out] INFO
130*> \verbatim
131*>          INFO is INTEGER
132*>          = 0:  successful exit.
133*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
134*>          > 0:  if INFO = 1, an eigenvalue did not converge
135*> \endverbatim
136*
137*  Authors:
138*  ========
139*
140*> \author Univ. of Tennessee
141*> \author Univ. of California Berkeley
142*> \author Univ. of Colorado Denver
143*> \author NAG Ltd.
144*
145*> \date December 2016
146*
147*> \ingroup auxOTHERcomputational
148*
149*> \par Contributors:
150*  ==================
151*>
152*> Jeff Rutter, Computer Science Division, University of California
153*> at Berkeley, USA
154*
155*  =====================================================================
156      SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
157     $                   S, LDS, INFO )
158*
159*  -- LAPACK computational routine (version 3.7.0) --
160*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
161*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162*     December 2016
163*
164*     .. Scalar Arguments ..
165      INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
166      DOUBLE PRECISION   RHO
167*     ..
168*     .. Array Arguments ..
169      DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
170     $                   W( * )
171*     ..
172*
173*  =====================================================================
174*
175*     .. Local Scalars ..
176      INTEGER            I, J
177      DOUBLE PRECISION   TEMP
178*     ..
179*     .. External Functions ..
180      DOUBLE PRECISION   DLAMC3, DNRM2
181      EXTERNAL           DLAMC3, DNRM2
182*     ..
183*     .. External Subroutines ..
184      EXTERNAL           DCOPY, DLAED4, XERBLA
185*     ..
186*     .. Intrinsic Functions ..
187      INTRINSIC          MAX, SIGN, SQRT
188*     ..
189*     .. Executable Statements ..
190*
191*     Test the input parameters.
192*
193      INFO = 0
194*
195      IF( K.LT.0 ) THEN
196         INFO = -1
197      ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
198         INFO = -2
199      ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
200     $          THEN
201         INFO = -3
202      ELSE IF( N.LT.K ) THEN
203         INFO = -4
204      ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
205         INFO = -7
206      ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
207         INFO = -12
208      END IF
209      IF( INFO.NE.0 ) THEN
210         CALL XERBLA( 'DLAED9', -INFO )
211         RETURN
212      END IF
213*
214*     Quick return if possible
215*
216      IF( K.EQ.0 )
217     $   RETURN
218*
219*     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
220*     be computed with high relative accuracy (barring over/underflow).
221*     This is a problem on machines without a guard digit in
222*     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
223*     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
224*     which on any of these machines zeros out the bottommost
225*     bit of DLAMDA(I) if it is 1; this makes the subsequent
226*     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
227*     occurs. On binary machines with a guard digit (almost all
228*     machines) it does not change DLAMDA(I) at all. On hexadecimal
229*     and decimal machines with a guard digit, it slightly
230*     changes the bottommost bits of DLAMDA(I). It does not account
231*     for hexadecimal or decimal machines without guard digits
232*     (we know of none). We use a subroutine call to compute
233*     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
234*     this code.
235*
236      DO 10 I = 1, N
237         DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
238   10 CONTINUE
239*
240      DO 20 J = KSTART, KSTOP
241         CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
242*
243*        If the zero finder fails, the computation is terminated.
244*
245         IF( INFO.NE.0 )
246     $      GO TO 120
247   20 CONTINUE
248*
249      IF( K.EQ.1 .OR. K.EQ.2 ) THEN
250         DO 40 I = 1, K
251            DO 30 J = 1, K
252               S( J, I ) = Q( J, I )
253   30       CONTINUE
254   40    CONTINUE
255         GO TO 120
256      END IF
257*
258*     Compute updated W.
259*
260      CALL DCOPY( K, W, 1, S, 1 )
261*
262*     Initialize W(I) = Q(I,I)
263*
264      CALL DCOPY( K, Q, LDQ+1, W, 1 )
265      DO 70 J = 1, K
266         DO 50 I = 1, J - 1
267            W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
268   50    CONTINUE
269         DO 60 I = J + 1, K
270            W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
271   60    CONTINUE
272   70 CONTINUE
273      DO 80 I = 1, K
274         W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
275   80 CONTINUE
276*
277*     Compute eigenvectors of the modified rank-1 modification.
278*
279      DO 110 J = 1, K
280         DO 90 I = 1, K
281            Q( I, J ) = W( I ) / Q( I, J )
282   90    CONTINUE
283         TEMP = DNRM2( K, Q( 1, J ), 1 )
284         DO 100 I = 1, K
285            S( I, J ) = Q( I, J ) / TEMP
286  100    CONTINUE
287  110 CONTINUE
288*
289  120 CONTINUE
290      RETURN
291*
292*     End of DLAED9
293*
294      END
295