1*> \brief \b DLAED3 used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLAED3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
22*                          CTOT, W, S, INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            INFO, K, LDQ, N, N1
26*       DOUBLE PRECISION   RHO
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            CTOT( * ), INDX( * )
30*       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
31*      $                   S( * ), W( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> DLAED3 finds the roots of the secular equation, as defined by the
41*> values in D, W, and RHO, between 1 and K.  It makes the
42*> appropriate calls to DLAED4 and then updates the eigenvectors by
43*> multiplying the matrix of eigenvectors of the pair of eigensystems
44*> being combined by the matrix of eigenvectors of the K-by-K system
45*> which is solved here.
46*>
47*> This code makes very mild assumptions about floating point
48*> arithmetic. It will work on machines with a guard digit in
49*> add/subtract, or on those binary machines without guard digits
50*> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
51*> It could conceivably fail on hexadecimal or decimal machines
52*> without guard digits, but we know of none.
53*> \endverbatim
54*
55*  Arguments:
56*  ==========
57*
58*> \param[in] K
59*> \verbatim
60*>          K is INTEGER
61*>          The number of terms in the rational function to be solved by
62*>          DLAED4.  K >= 0.
63*> \endverbatim
64*>
65*> \param[in] N
66*> \verbatim
67*>          N is INTEGER
68*>          The number of rows and columns in the Q matrix.
69*>          N >= K (deflation may result in N>K).
70*> \endverbatim
71*>
72*> \param[in] N1
73*> \verbatim
74*>          N1 is INTEGER
75*>          The location of the last eigenvalue in the leading submatrix.
76*>          min(1,N) <= N1 <= N/2.
77*> \endverbatim
78*>
79*> \param[out] D
80*> \verbatim
81*>          D is DOUBLE PRECISION array, dimension (N)
82*>          D(I) contains the updated eigenvalues for
83*>          1 <= I <= K.
84*> \endverbatim
85*>
86*> \param[out] Q
87*> \verbatim
88*>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
89*>          Initially the first K columns are used as workspace.
90*>          On output the columns 1 to K contain
91*>          the updated eigenvectors.
92*> \endverbatim
93*>
94*> \param[in] LDQ
95*> \verbatim
96*>          LDQ is INTEGER
97*>          The leading dimension of the array Q.  LDQ >= max(1,N).
98*> \endverbatim
99*>
100*> \param[in] RHO
101*> \verbatim
102*>          RHO is DOUBLE PRECISION
103*>          The value of the parameter in the rank one update equation.
104*>          RHO >= 0 required.
105*> \endverbatim
106*>
107*> \param[in,out] DLAMDA
108*> \verbatim
109*>          DLAMDA is DOUBLE PRECISION array, dimension (K)
110*>          The first K elements of this array contain the old roots
111*>          of the deflated updating problem.  These are the poles
112*>          of the secular equation. May be changed on output by
113*>          having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
114*>          Cray-2, or Cray C-90, as described above.
115*> \endverbatim
116*>
117*> \param[in] Q2
118*> \verbatim
119*>          Q2 is DOUBLE PRECISION array, dimension (LDQ2*N)
120*>          The first K columns of this matrix contain the non-deflated
121*>          eigenvectors for the split problem.
122*> \endverbatim
123*>
124*> \param[in] INDX
125*> \verbatim
126*>          INDX is INTEGER array, dimension (N)
127*>          The permutation used to arrange the columns of the deflated
128*>          Q matrix into three groups (see DLAED2).
129*>          The rows of the eigenvectors found by DLAED4 must be likewise
130*>          permuted before the matrix multiply can take place.
131*> \endverbatim
132*>
133*> \param[in] CTOT
134*> \verbatim
135*>          CTOT is INTEGER array, dimension (4)
136*>          A count of the total number of the various types of columns
137*>          in Q, as described in INDX.  The fourth column type is any
138*>          column which has been deflated.
139*> \endverbatim
140*>
141*> \param[in,out] W
142*> \verbatim
143*>          W is DOUBLE PRECISION array, dimension (K)
144*>          The first K elements of this array contain the components
145*>          of the deflation-adjusted updating vector. Destroyed on
146*>          output.
147*> \endverbatim
148*>
149*> \param[out] S
150*> \verbatim
151*>          S is DOUBLE PRECISION array, dimension (N1 + 1)*K
152*>          Will contain the eigenvectors of the repaired matrix which
153*>          will be multiplied by the previously accumulated eigenvectors
154*>          to update the system.
155*> \endverbatim
156*>
157*> \param[out] INFO
158*> \verbatim
159*>          INFO is INTEGER
160*>          = 0:  successful exit.
161*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
162*>          > 0:  if INFO = 1, an eigenvalue did not converge
163*> \endverbatim
164*
165*  Authors:
166*  ========
167*
168*> \author Univ. of Tennessee
169*> \author Univ. of California Berkeley
170*> \author Univ. of Colorado Denver
171*> \author NAG Ltd.
172*
173*> \date June 2017
174*
175*> \ingroup auxOTHERcomputational
176*
177*> \par Contributors:
178*  ==================
179*>
180*> Jeff Rutter, Computer Science Division, University of California
181*> at Berkeley, USA \n
182*>  Modified by Francoise Tisseur, University of Tennessee
183*>
184*  =====================================================================
185      SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
186     $                   CTOT, W, S, INFO )
187*
188*  -- LAPACK computational routine (version 3.7.1) --
189*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
190*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191*     June 2017
192*
193*     .. Scalar Arguments ..
194      INTEGER            INFO, K, LDQ, N, N1
195      DOUBLE PRECISION   RHO
196*     ..
197*     .. Array Arguments ..
198      INTEGER            CTOT( * ), INDX( * )
199      DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
200     $                   S( * ), W( * )
201*     ..
202*
203*  =====================================================================
204*
205*     .. Parameters ..
206      DOUBLE PRECISION   ONE, ZERO
207      PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
208*     ..
209*     .. Local Scalars ..
210      INTEGER            I, II, IQ2, J, N12, N2, N23
211      DOUBLE PRECISION   TEMP
212*     ..
213*     .. External Functions ..
214      DOUBLE PRECISION   DLAMC3, DNRM2
215      EXTERNAL           DLAMC3, DNRM2
216*     ..
217*     .. External Subroutines ..
218      EXTERNAL           DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA
219*     ..
220*     .. Intrinsic Functions ..
221      INTRINSIC          MAX, SIGN, SQRT
222*     ..
223*     .. Executable Statements ..
224*
225*     Test the input parameters.
226*
227      INFO = 0
228*
229      IF( K.LT.0 ) THEN
230         INFO = -1
231      ELSE IF( N.LT.K ) THEN
232         INFO = -2
233      ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
234         INFO = -6
235      END IF
236      IF( INFO.NE.0 ) THEN
237         CALL XERBLA( 'DLAED3', -INFO )
238         RETURN
239      END IF
240*
241*     Quick return if possible
242*
243      IF( K.EQ.0 )
244     $   RETURN
245*
246*     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
247*     be computed with high relative accuracy (barring over/underflow).
248*     This is a problem on machines without a guard digit in
249*     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
250*     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
251*     which on any of these machines zeros out the bottommost
252*     bit of DLAMDA(I) if it is 1; this makes the subsequent
253*     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
254*     occurs. On binary machines with a guard digit (almost all
255*     machines) it does not change DLAMDA(I) at all. On hexadecimal
256*     and decimal machines with a guard digit, it slightly
257*     changes the bottommost bits of DLAMDA(I). It does not account
258*     for hexadecimal or decimal machines without guard digits
259*     (we know of none). We use a subroutine call to compute
260*     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
261*     this code.
262*
263      DO 10 I = 1, K
264         DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
265   10 CONTINUE
266*
267      DO 20 J = 1, K
268         CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
269*
270*        If the zero finder fails, the computation is terminated.
271*
272         IF( INFO.NE.0 )
273     $      GO TO 120
274   20 CONTINUE
275*
276      IF( K.EQ.1 )
277     $   GO TO 110
278      IF( K.EQ.2 ) THEN
279         DO 30 J = 1, K
280            W( 1 ) = Q( 1, J )
281            W( 2 ) = Q( 2, J )
282            II = INDX( 1 )
283            Q( 1, J ) = W( II )
284            II = INDX( 2 )
285            Q( 2, J ) = W( II )
286   30    CONTINUE
287         GO TO 110
288      END IF
289*
290*     Compute updated W.
291*
292      CALL DCOPY( K, W, 1, S, 1 )
293*
294*     Initialize W(I) = Q(I,I)
295*
296      CALL DCOPY( K, Q, LDQ+1, W, 1 )
297      DO 60 J = 1, K
298         DO 40 I = 1, J - 1
299            W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
300   40    CONTINUE
301         DO 50 I = J + 1, K
302            W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
303   50    CONTINUE
304   60 CONTINUE
305      DO 70 I = 1, K
306         W( I ) = SIGN( SQRT( -W( I ) ), S( I ) )
307   70 CONTINUE
308*
309*     Compute eigenvectors of the modified rank-1 modification.
310*
311      DO 100 J = 1, K
312         DO 80 I = 1, K
313            S( I ) = W( I ) / Q( I, J )
314   80    CONTINUE
315         TEMP = DNRM2( K, S, 1 )
316         DO 90 I = 1, K
317            II = INDX( I )
318            Q( I, J ) = S( II ) / TEMP
319   90    CONTINUE
320  100 CONTINUE
321*
322*     Compute the updated eigenvectors.
323*
324  110 CONTINUE
325*
326      N2 = N - N1
327      N12 = CTOT( 1 ) + CTOT( 2 )
328      N23 = CTOT( 2 ) + CTOT( 3 )
329*
330      CALL DLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 )
331      IQ2 = N1*N12 + 1
332      IF( N23.NE.0 ) THEN
333         CALL DGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23,
334     $               ZERO, Q( N1+1, 1 ), LDQ )
335      ELSE
336         CALL DLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ )
337      END IF
338*
339      CALL DLACPY( 'A', N12, K, Q, LDQ, S, N12 )
340      IF( N12.NE.0 ) THEN
341         CALL DGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q,
342     $               LDQ )
343      ELSE
344         CALL DLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ )
345      END IF
346*
347*
348  120 CONTINUE
349      RETURN
350*
351*     End of DLAED3
352*
353      END
354