LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slaed9.f
Go to the documentation of this file.
1 *> \brief \b SLAED9 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 SLAED9 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed9.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed9.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed9.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAED9( 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 * REAL RHO
27 * ..
28 * .. Array Arguments ..
29 * REAL D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
30 * $ W( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SLAED9 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 SLAED4 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 *> SLAED4. 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 REAL 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 REAL 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 REAL
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 REAL 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 REAL 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 REAL 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 September 2012
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 slaed9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
157  $ s, lds, info )
158 *
159 * -- LAPACK computational routine (version 3.4.2) --
160 * -- LAPACK is a software package provided by Univ. of Tennessee, --
161 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162 * September 2012
163 *
164 * .. Scalar Arguments ..
165  INTEGER info, k, kstart, kstop, ldq, lds, n
166  REAL rho
167 * ..
168 * .. Array Arguments ..
169  REAL d( * ), dlamda( * ), q( ldq, * ), s( lds, * ),
170  $ w( * )
171 * ..
172 *
173 * =====================================================================
174 *
175 * .. Local Scalars ..
176  INTEGER i, j
177  REAL temp
178 * ..
179 * .. External Functions ..
180  REAL slamc3, snrm2
181  EXTERNAL slamc3, snrm2
182 * ..
183 * .. External Subroutines ..
184  EXTERNAL scopy, slaed4, 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( 'SLAED9', -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 ) = slamc3( dlamda( i ), dlamda( i ) ) - dlamda( i )
238  10 continue
239 *
240  DO 20 j = kstart, kstop
241  CALL slaed4( 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 scopy( k, w, 1, s, 1 )
261 *
262 * Initialize W(I) = Q(I,I)
263 *
264  CALL scopy( 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 = snrm2( 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 SLAED9
293 *
294  END