LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> \ingroup auxOTHERcomputational
146 *
147 *> \par Contributors:
148 * ==================
149 *>
150 *> Jeff Rutter, Computer Science Division, University of California
151 *> at Berkeley, USA
152 *
153 * =====================================================================
154  SUBROUTINE slaed9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
155  $ S, LDS, INFO )
156 *
157 * -- LAPACK computational routine --
158 * -- LAPACK is a software package provided by Univ. of Tennessee, --
159 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160 *
161 * .. Scalar Arguments ..
162  INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N
163  REAL RHO
164 * ..
165 * .. Array Arguments ..
166  REAL D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
167  $ w( * )
168 * ..
169 *
170 * =====================================================================
171 *
172 * .. Local Scalars ..
173  INTEGER I, J
174  REAL TEMP
175 * ..
176 * .. External Functions ..
177  REAL SLAMC3, SNRM2
178  EXTERNAL slamc3, snrm2
179 * ..
180 * .. External Subroutines ..
181  EXTERNAL scopy, slaed4, xerbla
182 * ..
183 * .. Intrinsic Functions ..
184  INTRINSIC max, sign, sqrt
185 * ..
186 * .. Executable Statements ..
187 *
188 * Test the input parameters.
189 *
190  info = 0
191 *
192  IF( k.LT.0 ) THEN
193  info = -1
194  ELSE IF( kstart.LT.1 .OR. kstart.GT.max( 1, k ) ) THEN
195  info = -2
196  ELSE IF( max( 1, kstop ).LT.kstart .OR. kstop.GT.max( 1, k ) )
197  $ THEN
198  info = -3
199  ELSE IF( n.LT.k ) THEN
200  info = -4
201  ELSE IF( ldq.LT.max( 1, k ) ) THEN
202  info = -7
203  ELSE IF( lds.LT.max( 1, k ) ) THEN
204  info = -12
205  END IF
206  IF( info.NE.0 ) THEN
207  CALL xerbla( 'SLAED9', -info )
208  RETURN
209  END IF
210 *
211 * Quick return if possible
212 *
213  IF( k.EQ.0 )
214  $ RETURN
215 *
216 * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
217 * be computed with high relative accuracy (barring over/underflow).
218 * This is a problem on machines without a guard digit in
219 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
220 * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
221 * which on any of these machines zeros out the bottommost
222 * bit of DLAMDA(I) if it is 1; this makes the subsequent
223 * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
224 * occurs. On binary machines with a guard digit (almost all
225 * machines) it does not change DLAMDA(I) at all. On hexadecimal
226 * and decimal machines with a guard digit, it slightly
227 * changes the bottommost bits of DLAMDA(I). It does not account
228 * for hexadecimal or decimal machines without guard digits
229 * (we know of none). We use a subroutine call to compute
230 * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
231 * this code.
232 *
233  DO 10 i = 1, n
234  dlamda( i ) = slamc3( dlamda( i ), dlamda( i ) ) - dlamda( i )
235  10 CONTINUE
236 *
237  DO 20 j = kstart, kstop
238  CALL slaed4( k, j, dlamda, w, q( 1, j ), rho, d( j ), info )
239 *
240 * If the zero finder fails, the computation is terminated.
241 *
242  IF( info.NE.0 )
243  $ GO TO 120
244  20 CONTINUE
245 *
246  IF( k.EQ.1 .OR. k.EQ.2 ) THEN
247  DO 40 i = 1, k
248  DO 30 j = 1, k
249  s( j, i ) = q( j, i )
250  30 CONTINUE
251  40 CONTINUE
252  GO TO 120
253  END IF
254 *
255 * Compute updated W.
256 *
257  CALL scopy( k, w, 1, s, 1 )
258 *
259 * Initialize W(I) = Q(I,I)
260 *
261  CALL scopy( k, q, ldq+1, w, 1 )
262  DO 70 j = 1, k
263  DO 50 i = 1, j - 1
264  w( i ) = w( i )*( q( i, j ) / ( dlamda( i )-dlamda( j ) ) )
265  50 CONTINUE
266  DO 60 i = j + 1, k
267  w( i ) = w( i )*( q( i, j ) / ( dlamda( i )-dlamda( j ) ) )
268  60 CONTINUE
269  70 CONTINUE
270  DO 80 i = 1, k
271  w( i ) = sign( sqrt( -w( i ) ), s( i, 1 ) )
272  80 CONTINUE
273 *
274 * Compute eigenvectors of the modified rank-1 modification.
275 *
276  DO 110 j = 1, k
277  DO 90 i = 1, k
278  q( i, j ) = w( i ) / q( i, j )
279  90 CONTINUE
280  temp = snrm2( k, q( 1, j ), 1 )
281  DO 100 i = 1, k
282  s( i, j ) = q( i, j ) / temp
283  100 CONTINUE
284  110 CONTINUE
285 *
286  120 CONTINUE
287  RETURN
288 *
289 * End of SLAED9
290 *
291  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slaed9(K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO)
SLAED9 used by SSTEDC. Finds the roots of the secular equation and updates the eigenvectors....
Definition: slaed9.f:156
subroutine slaed4(N, I, D, Z, DELTA, RHO, DLAM, INFO)
SLAED4 used by SSTEDC. Finds a single root of the secular equation.
Definition: slaed4.f:145
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82