LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slaed7.f
Go to the documentation of this file.
1 *> \brief \b SLAED7 used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. 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 SLAED7 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed7.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed7.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed7.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
22 * LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR,
23 * PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK,
24 * INFO )
25 *
26 * .. Scalar Arguments ..
27 * INTEGER CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
28 * $ QSIZ, TLVLS
29 * REAL RHO
30 * ..
31 * .. Array Arguments ..
32 * INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
33 * $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
34 * REAL D( * ), GIVNUM( 2, * ), Q( LDQ, * ),
35 * $ QSTORE( * ), WORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> SLAED7 computes the updated eigensystem of a diagonal
45 *> matrix after modification by a rank-one symmetric matrix. This
46 *> routine is used only for the eigenproblem which requires all
47 *> eigenvalues and optionally eigenvectors of a dense symmetric matrix
48 *> that has been reduced to tridiagonal form. SLAED1 handles
49 *> the case in which all eigenvalues and eigenvectors of a symmetric
50 *> tridiagonal matrix are desired.
51 *>
52 *> T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
53 *>
54 *> where Z = Q**Tu, u is a vector of length N with ones in the
55 *> CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
56 *>
57 *> The eigenvectors of the original matrix are stored in Q, and the
58 *> eigenvalues are in D. The algorithm consists of three stages:
59 *>
60 *> The first stage consists of deflating the size of the problem
61 *> when there are multiple eigenvalues or if there is a zero in
62 *> the Z vector. For each such occurence the dimension of the
63 *> secular equation problem is reduced by one. This stage is
64 *> performed by the routine SLAED8.
65 *>
66 *> The second stage consists of calculating the updated
67 *> eigenvalues. This is done by finding the roots of the secular
68 *> equation via the routine SLAED4 (as called by SLAED9).
69 *> This routine also calculates the eigenvectors of the current
70 *> problem.
71 *>
72 *> The final stage consists of computing the updated eigenvectors
73 *> directly using the updated eigenvalues. The eigenvectors for
74 *> the current problem are multiplied with the eigenvectors from
75 *> the overall problem.
76 *> \endverbatim
77 *
78 * Arguments:
79 * ==========
80 *
81 *> \param[in] ICOMPQ
82 *> \verbatim
83 *> ICOMPQ is INTEGER
84 *> = 0: Compute eigenvalues only.
85 *> = 1: Compute eigenvectors of original dense symmetric matrix
86 *> also. On entry, Q contains the orthogonal matrix used
87 *> to reduce the original matrix to tridiagonal form.
88 *> \endverbatim
89 *>
90 *> \param[in] N
91 *> \verbatim
92 *> N is INTEGER
93 *> The dimension of the symmetric tridiagonal matrix. N >= 0.
94 *> \endverbatim
95 *>
96 *> \param[in] QSIZ
97 *> \verbatim
98 *> QSIZ is INTEGER
99 *> The dimension of the orthogonal matrix used to reduce
100 *> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
101 *> \endverbatim
102 *>
103 *> \param[in] TLVLS
104 *> \verbatim
105 *> TLVLS is INTEGER
106 *> The total number of merging levels in the overall divide and
107 *> conquer tree.
108 *> \endverbatim
109 *>
110 *> \param[in] CURLVL
111 *> \verbatim
112 *> CURLVL is INTEGER
113 *> The current level in the overall merge routine,
114 *> 0 <= CURLVL <= TLVLS.
115 *> \endverbatim
116 *>
117 *> \param[in] CURPBM
118 *> \verbatim
119 *> CURPBM is INTEGER
120 *> The current problem in the current level in the overall
121 *> merge routine (counting from upper left to lower right).
122 *> \endverbatim
123 *>
124 *> \param[in,out] D
125 *> \verbatim
126 *> D is REAL array, dimension (N)
127 *> On entry, the eigenvalues of the rank-1-perturbed matrix.
128 *> On exit, the eigenvalues of the repaired matrix.
129 *> \endverbatim
130 *>
131 *> \param[in,out] Q
132 *> \verbatim
133 *> Q is REAL array, dimension (LDQ, N)
134 *> On entry, the eigenvectors of the rank-1-perturbed matrix.
135 *> On exit, the eigenvectors of the repaired tridiagonal matrix.
136 *> \endverbatim
137 *>
138 *> \param[in] LDQ
139 *> \verbatim
140 *> LDQ is INTEGER
141 *> The leading dimension of the array Q. LDQ >= max(1,N).
142 *> \endverbatim
143 *>
144 *> \param[out] INDXQ
145 *> \verbatim
146 *> INDXQ is INTEGER array, dimension (N)
147 *> The permutation which will reintegrate the subproblem just
148 *> solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
149 *> will be in ascending order.
150 *> \endverbatim
151 *>
152 *> \param[in] RHO
153 *> \verbatim
154 *> RHO is REAL
155 *> The subdiagonal element used to create the rank-1
156 *> modification.
157 *> \endverbatim
158 *>
159 *> \param[in] CUTPNT
160 *> \verbatim
161 *> CUTPNT is INTEGER
162 *> Contains the location of the last eigenvalue in the leading
163 *> sub-matrix. min(1,N) <= CUTPNT <= N.
164 *> \endverbatim
165 *>
166 *> \param[in,out] QSTORE
167 *> \verbatim
168 *> QSTORE is REAL array, dimension (N**2+1)
169 *> Stores eigenvectors of submatrices encountered during
170 *> divide and conquer, packed together. QPTR points to
171 *> beginning of the submatrices.
172 *> \endverbatim
173 *>
174 *> \param[in,out] QPTR
175 *> \verbatim
176 *> QPTR is INTEGER array, dimension (N+2)
177 *> List of indices pointing to beginning of submatrices stored
178 *> in QSTORE. The submatrices are numbered starting at the
179 *> bottom left of the divide and conquer tree, from left to
180 *> right and bottom to top.
181 *> \endverbatim
182 *>
183 *> \param[in] PRMPTR
184 *> \verbatim
185 *> PRMPTR is INTEGER array, dimension (N lg N)
186 *> Contains a list of pointers which indicate where in PERM a
187 *> level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
188 *> indicates the size of the permutation and also the size of
189 *> the full, non-deflated problem.
190 *> \endverbatim
191 *>
192 *> \param[in] PERM
193 *> \verbatim
194 *> PERM is INTEGER array, dimension (N lg N)
195 *> Contains the permutations (from deflation and sorting) to be
196 *> applied to each eigenblock.
197 *> \endverbatim
198 *>
199 *> \param[in] GIVPTR
200 *> \verbatim
201 *> GIVPTR is INTEGER array, dimension (N lg N)
202 *> Contains a list of pointers which indicate where in GIVCOL a
203 *> level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
204 *> indicates the number of Givens rotations.
205 *> \endverbatim
206 *>
207 *> \param[in] GIVCOL
208 *> \verbatim
209 *> GIVCOL is INTEGER array, dimension (2, N lg N)
210 *> Each pair of numbers indicates a pair of columns to take place
211 *> in a Givens rotation.
212 *> \endverbatim
213 *>
214 *> \param[in] GIVNUM
215 *> \verbatim
216 *> GIVNUM is REAL array, dimension (2, N lg N)
217 *> Each number indicates the S value to be used in the
218 *> corresponding Givens rotation.
219 *> \endverbatim
220 *>
221 *> \param[out] WORK
222 *> \verbatim
223 *> WORK is REAL array, dimension (3*N+2*QSIZ*N)
224 *> \endverbatim
225 *>
226 *> \param[out] IWORK
227 *> \verbatim
228 *> IWORK is INTEGER array, dimension (4*N)
229 *> \endverbatim
230 *>
231 *> \param[out] INFO
232 *> \verbatim
233 *> INFO is INTEGER
234 *> = 0: successful exit.
235 *> < 0: if INFO = -i, the i-th argument had an illegal value.
236 *> > 0: if INFO = 1, an eigenvalue did not converge
237 *> \endverbatim
238 *
239 * Authors:
240 * ========
241 *
242 *> \author Univ. of Tennessee
243 *> \author Univ. of California Berkeley
244 *> \author Univ. of Colorado Denver
245 *> \author NAG Ltd.
246 *
247 *> \date September 2012
248 *
249 *> \ingroup auxOTHERcomputational
250 *
251 *> \par Contributors:
252 * ==================
253 *>
254 *> Jeff Rutter, Computer Science Division, University of California
255 *> at Berkeley, USA
256 *
257 * =====================================================================
258  SUBROUTINE slaed7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
259  $ ldq, indxq, rho, cutpnt, qstore, qptr, prmptr,
260  $ perm, givptr, givcol, givnum, work, iwork,
261  $ info )
262 *
263 * -- LAPACK computational routine (version 3.4.2) --
264 * -- LAPACK is a software package provided by Univ. of Tennessee, --
265 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
266 * September 2012
267 *
268 * .. Scalar Arguments ..
269  INTEGER curlvl, curpbm, cutpnt, icompq, info, ldq, n,
270  $ qsiz, tlvls
271  REAL rho
272 * ..
273 * .. Array Arguments ..
274  INTEGER givcol( 2, * ), givptr( * ), indxq( * ),
275  $ iwork( * ), perm( * ), prmptr( * ), qptr( * )
276  REAL d( * ), givnum( 2, * ), q( ldq, * ),
277  $ qstore( * ), work( * )
278 * ..
279 *
280 * =====================================================================
281 *
282 * .. Parameters ..
283  REAL one, zero
284  parameter( one = 1.0e0, zero = 0.0e0 )
285 * ..
286 * .. Local Scalars ..
287  INTEGER coltyp, curr, i, idlmda, indx, indxc, indxp,
288  $ iq2, is, iw, iz, k, ldq2, n1, n2, ptr
289 * ..
290 * .. External Subroutines ..
291  EXTERNAL sgemm, slaed8, slaed9, slaeda, slamrg, xerbla
292 * ..
293 * .. Intrinsic Functions ..
294  INTRINSIC max, min
295 * ..
296 * .. Executable Statements ..
297 *
298 * Test the input parameters.
299 *
300  info = 0
301 *
302  IF( icompq.LT.0 .OR. icompq.GT.1 ) THEN
303  info = -1
304  ELSE IF( n.LT.0 ) THEN
305  info = -2
306  ELSE IF( icompq.EQ.1 .AND. qsiz.LT.n ) THEN
307  info = -4
308  ELSE IF( ldq.LT.max( 1, n ) ) THEN
309  info = -9
310  ELSE IF( min( 1, n ).GT.cutpnt .OR. n.LT.cutpnt ) THEN
311  info = -12
312  END IF
313  IF( info.NE.0 ) THEN
314  CALL xerbla( 'SLAED7', -info )
315  return
316  END IF
317 *
318 * Quick return if possible
319 *
320  IF( n.EQ.0 )
321  $ return
322 *
323 * The following values are for bookkeeping purposes only. They are
324 * integer pointers which indicate the portion of the workspace
325 * used by a particular array in SLAED8 and SLAED9.
326 *
327  IF( icompq.EQ.1 ) THEN
328  ldq2 = qsiz
329  ELSE
330  ldq2 = n
331  END IF
332 *
333  iz = 1
334  idlmda = iz + n
335  iw = idlmda + n
336  iq2 = iw + n
337  is = iq2 + n*ldq2
338 *
339  indx = 1
340  indxc = indx + n
341  coltyp = indxc + n
342  indxp = coltyp + n
343 *
344 * Form the z-vector which consists of the last row of Q_1 and the
345 * first row of Q_2.
346 *
347  ptr = 1 + 2**tlvls
348  DO 10 i = 1, curlvl - 1
349  ptr = ptr + 2**( tlvls-i )
350  10 continue
351  curr = ptr + curpbm
352  CALL slaeda( n, tlvls, curlvl, curpbm, prmptr, perm, givptr,
353  $ givcol, givnum, qstore, qptr, work( iz ),
354  $ work( iz+n ), info )
355 *
356 * When solving the final problem, we no longer need the stored data,
357 * so we will overwrite the data from this level onto the previously
358 * used storage space.
359 *
360  IF( curlvl.EQ.tlvls ) THEN
361  qptr( curr ) = 1
362  prmptr( curr ) = 1
363  givptr( curr ) = 1
364  END IF
365 *
366 * Sort and Deflate eigenvalues.
367 *
368  CALL slaed8( icompq, k, n, qsiz, d, q, ldq, indxq, rho, cutpnt,
369  $ work( iz ), work( idlmda ), work( iq2 ), ldq2,
370  $ work( iw ), perm( prmptr( curr ) ), givptr( curr+1 ),
371  $ givcol( 1, givptr( curr ) ),
372  $ givnum( 1, givptr( curr ) ), iwork( indxp ),
373  $ iwork( indx ), info )
374  prmptr( curr+1 ) = prmptr( curr ) + n
375  givptr( curr+1 ) = givptr( curr+1 ) + givptr( curr )
376 *
377 * Solve Secular Equation.
378 *
379  IF( k.NE.0 ) THEN
380  CALL slaed9( k, 1, k, n, d, work( is ), k, rho, work( idlmda ),
381  $ work( iw ), qstore( qptr( curr ) ), k, info )
382  IF( info.NE.0 )
383  $ go to 30
384  IF( icompq.EQ.1 ) THEN
385  CALL sgemm( 'N', 'N', qsiz, k, k, one, work( iq2 ), ldq2,
386  $ qstore( qptr( curr ) ), k, zero, q, ldq )
387  END IF
388  qptr( curr+1 ) = qptr( curr ) + k**2
389 *
390 * Prepare the INDXQ sorting permutation.
391 *
392  n1 = k
393  n2 = n - k
394  CALL slamrg( n1, n2, d, 1, -1, indxq )
395  ELSE
396  qptr( curr+1 ) = qptr( curr )
397  DO 20 i = 1, n
398  indxq( i ) = i
399  20 continue
400  END IF
401 *
402  30 continue
403  return
404 *
405 * End of SLAED7
406 *
407  END