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.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SLAED7 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed7.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed7.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed7.f">
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,
26 * .. Scalar Arguments ..
27 * INTEGER CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
31 * .. Array Arguments ..
32 * INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
33 * $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
34 * REAL D( * ), GIVNUM( 2, * ), Q( LDQ, * ),
35 * $ QSTORE( * ), WORK( * )
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.
52 *> T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
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.
57 *> The eigenvectors of the original matrix are stored in Q, and the
58 *> eigenvalues are in D. The algorithm consists of three stages:
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 occurrence the dimension of the
63 *> secular equation problem is reduced by one. This stage is
64 *> performed by the routine SLAED8.
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
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.
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.
93 *> The dimension of the symmetric tridiagonal matrix. N >= 0.
99 *> The dimension of the orthogonal matrix used to reduce
100 *> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
106 *> The total number of merging levels in the overall divide and
113 *> The current level in the overall merge routine,
114 *> 0 <= CURLVL <= TLVLS.
120 *> The current problem in the current level in the overall
121 *> merge routine (counting from upper left to lower right).
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.
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.
141 *> The leading dimension of the array Q. LDQ >= max(1,N).
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.
155 *> The subdiagonal element used to create the rank-1
162 *> Contains the location of the last eigenvalue in the leading
163 *> sub-matrix. min(1,N) <= CUTPNT <= N.
166 *> \param[in,out] QSTORE
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.
174 *> \param[in,out] QPTR
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.
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.
194 *> PERM is INTEGER array, dimension (N lg N)
195 *> Contains the permutations (from deflation and sorting) to be
196 *> applied to each eigenblock.
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.
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.
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.
223 *> WORK is REAL array, dimension (3*N+2*QSIZ*N)
228 *> IWORK is INTEGER array, dimension (4*N)
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
242 *> \author Univ. of Tennessee
243 *> \author Univ. of California Berkeley
244 *> \author Univ. of Colorado Denver
249 *> \ingroup auxOTHERcomputational
251 *> \par Contributors:
254 *> Jeff Rutter, Computer Science Division, University of California
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,
263 * -- LAPACK computational routine (version 3.6.1) --
264 * -- LAPACK is a software package provided by Univ. of Tennessee, --
265 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * .. Scalar Arguments ..
269 INTEGER CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
273 * .. Array Arguments ..
274 INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
275 $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
276 REAL D( * ), GIVNUM( 2, * ), Q( LDQ, * ),
277 $ QSTORE( * ), WORK( * )
280 * =====================================================================
284 PARAMETER ( ONE = 1.0E0, ZERO = 0.0E0 )
286 * .. Local Scalars ..
287 INTEGER COLTYP, CURR, I, IDLMDA, INDX, INDXC, INDXP,
288 $ IQ2, IS, IW, IZ, K, LDQ2, N1, N2, PTR
290 * .. External Subroutines ..
291 EXTERNAL SGEMM, SLAED8, SLAED9, SLAEDA, SLAMRG, XERBLA
293 * .. Intrinsic Functions ..
296 * .. Executable Statements ..
298 * Test the input parameters.
302 IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
304 ELSE IF( N.LT.0 ) THEN
306 ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
308 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
310 ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
314 CALL XERBLA( 'SLAED7', -INFO )
318 * Quick return if possible
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.
327 IF( ICOMPQ.EQ.1 ) THEN
344 * Form the z-vector which consists of the last row of Q_1 and the
348 DO 10 I = 1, CURLVL - 1
349 PTR = PTR + 2**( TLVLS-I )
352 CALL SLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
353 $ GIVCOL, GIVNUM, QSTORE, QPTR, WORK( IZ ),
354 $ WORK( IZ+N ), INFO )
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.
360 IF( CURLVL.EQ.TLVLS ) THEN
366 * Sort and Deflate eigenvalues.
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 )
377 * Solve Secular Equation.
380 CALL SLAED9( K, 1, K, N, D, WORK( IS ), K, RHO, WORK( IDLMDA ),
381 $ WORK( IW ), QSTORE( QPTR( CURR ) ), K, INFO )
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 )
388 QPTR( CURR+1 ) = QPTR( CURR ) + K**2
390 * Prepare the INDXQ sorting permutation.
394 CALL SLAMRG( N1, N2, D, 1, -1, INDXQ )
396 QPTR( CURR+1 ) = QPTR( CURR )