LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
slaed8.f
Go to the documentation of this file.
1 *> \brief \b SLAED8 used by SSTEDC. 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 SLAED8 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed8.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed8.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed8.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAED8( 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 * REAL RHO
29 * ..
30 * .. Array Arguments ..
31 * INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
32 * $ INDXQ( * ), PERM( * )
33 * REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ),
34 * $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> SLAED8 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 REAL 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 REAL 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 REAL
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 *> SLAED3.
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 REAL 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 REAL array, dimension (N)
147 *> A copy of the first K eigenvalues which will be used by
148 *> SLAED3 to form the secular equation.
149 *> \endverbatim
150 *>
151 *> \param[out] Q2
152 *> \verbatim
153 *> Q2 is REAL 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 *> SLAED7 in a matrix multiply (SGEMM) 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 REAL array, dimension (N)
169 *> The first k values of the final deflation-altered z-vector and
170 *> will be passed to SLAED3.
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 REAL 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 slaed8( 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  REAL RHO
252 * ..
253 * .. Array Arguments ..
254  INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
255  $ INDXQ( * ), PERM( * )
256  REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ),
257  $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
258 * ..
259 *
260 * =====================================================================
261 *
262 * .. Parameters ..
263  REAL MONE, ZERO, ONE, TWO, EIGHT
264  PARAMETER ( MONE = -1.0e0, zero = 0.0e0, one = 1.0e0,
265  $ two = 2.0e0, eight = 8.0e0 )
266 * ..
267 * .. Local Scalars ..
268 *
269  INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
270  REAL C, EPS, S, T, TAU, TOL
271 * ..
272 * .. External Functions ..
273  INTEGER ISAMAX
274  REAL SLAMCH, SLAPY2
275  EXTERNAL isamax, slamch, slapy2
276 * ..
277 * .. External Subroutines ..
278  EXTERNAL scopy, slacpy, slamrg, srot, sscal, 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( 'SLAED8', -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 sscal( 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 sscal( 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 slamrg( 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 = isamax( n, z, 1 )
356  jmax = isamax( n, d, 1 )
357  eps = slamch( '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 scopy( qsiz, q( 1, perm( j ) ), 1, q2( 1, j ), 1 )
374  60 CONTINUE
375  CALL slacpy( '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 = slapy2( 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 srot( 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 scopy( 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 scopy( n-k, dlamda( k+1 ), 1, d( k+1 ), 1 )
510  ELSE
511  CALL scopy( n-k, dlamda( k+1 ), 1, d( k+1 ), 1 )
512  CALL slacpy( '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 SLAED8
520 *
521  END
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slamrg(N1, N2, A, STRD1, STRD2, INDEX)
SLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single...
Definition: slamrg.f:99
subroutine slaed8(ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, INFO)
SLAED8 used by SSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matri...
Definition: slaed8.f:243
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79