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