3 /* Table of constant values */
5 static integer c__2 = 2;
6 static integer c__1 = 1;
7 static real c_b10 = 1.f;
8 static real c_b11 = 0.f;
9 static integer c_n1 = -1;
11 /* Subroutine */ int slaed7_(integer *icompq, integer *n, integer *qsiz,
12 integer *tlvls, integer *curlvl, integer *curpbm, real *d__, real *q,
13 integer *ldq, integer *indxq, real *rho, integer *cutpnt, real *
14 qstore, integer *qptr, integer *prmptr, integer *perm, integer *
15 givptr, integer *givcol, real *givnum, real *work, integer *iwork,
18 /* System generated locals */
19 integer q_dim1, q_offset, i__1, i__2;
21 /* Builtin functions */
22 integer pow_ii(integer *, integer *);
25 integer i__, k, n1, n2, is, iw, iz, iq2, ptr, ldq2, indx, curr, indxc;
26 extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
27 integer *, real *, real *, integer *, real *, integer *, real *,
30 extern /* Subroutine */ int slaed8_(integer *, integer *, integer *,
31 integer *, real *, real *, integer *, integer *, real *, integer *
32 , real *, real *, real *, integer *, real *, integer *, integer *,
33 integer *, real *, integer *, integer *, integer *), slaed9_(
34 integer *, integer *, integer *, integer *, real *, real *,
35 integer *, real *, real *, real *, real *, integer *, integer *),
36 slaeda_(integer *, integer *, integer *, integer *, integer *,
37 integer *, integer *, integer *, real *, real *, integer *, real *
40 extern /* Subroutine */ int xerbla_(char *, integer *), slamrg_(
41 integer *, integer *, real *, integer *, integer *, integer *);
45 /* -- LAPACK routine (version 3.1) -- */
46 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
49 /* .. Scalar Arguments .. */
51 /* .. Array Arguments .. */
57 /* SLAED7 computes the updated eigensystem of a diagonal */
58 /* matrix after modification by a rank-one symmetric matrix. This */
59 /* routine is used only for the eigenproblem which requires all */
60 /* eigenvalues and optionally eigenvectors of a dense symmetric matrix */
61 /* that has been reduced to tridiagonal form. SLAED1 handles */
62 /* the case in which all eigenvalues and eigenvectors of a symmetric */
63 /* tridiagonal matrix are desired. */
65 /* T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out) */
67 /* where Z = Q'u, u is a vector of length N with ones in the */
68 /* CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. */
70 /* The eigenvectors of the original matrix are stored in Q, and the */
71 /* eigenvalues are in D. The algorithm consists of three stages: */
73 /* The first stage consists of deflating the size of the problem */
74 /* when there are multiple eigenvalues or if there is a zero in */
75 /* the Z vector. For each such occurence the dimension of the */
76 /* secular equation problem is reduced by one. This stage is */
77 /* performed by the routine SLAED8. */
79 /* The second stage consists of calculating the updated */
80 /* eigenvalues. This is done by finding the roots of the secular */
81 /* equation via the routine SLAED4 (as called by SLAED9). */
82 /* This routine also calculates the eigenvectors of the current */
85 /* The final stage consists of computing the updated eigenvectors */
86 /* directly using the updated eigenvalues. The eigenvectors for */
87 /* the current problem are multiplied with the eigenvectors from */
88 /* the overall problem. */
93 /* ICOMPQ (input) INTEGER */
94 /* = 0: Compute eigenvalues only. */
95 /* = 1: Compute eigenvectors of original dense symmetric matrix */
96 /* also. On entry, Q contains the orthogonal matrix used */
97 /* to reduce the original matrix to tridiagonal form. */
99 /* N (input) INTEGER */
100 /* The dimension of the symmetric tridiagonal matrix. N >= 0. */
102 /* QSIZ (input) INTEGER */
103 /* The dimension of the orthogonal matrix used to reduce */
104 /* the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. */
106 /* TLVLS (input) INTEGER */
107 /* The total number of merging levels in the overall divide and */
110 /* CURLVL (input) INTEGER */
111 /* The current level in the overall merge routine, */
112 /* 0 <= CURLVL <= TLVLS. */
114 /* CURPBM (input) INTEGER */
115 /* The current problem in the current level in the overall */
116 /* merge routine (counting from upper left to lower right). */
118 /* D (input/output) REAL array, dimension (N) */
119 /* On entry, the eigenvalues of the rank-1-perturbed matrix. */
120 /* On exit, the eigenvalues of the repaired matrix. */
122 /* Q (input/output) REAL array, dimension (LDQ, N) */
123 /* On entry, the eigenvectors of the rank-1-perturbed matrix. */
124 /* On exit, the eigenvectors of the repaired tridiagonal matrix. */
126 /* LDQ (input) INTEGER */
127 /* The leading dimension of the array Q. LDQ >= max(1,N). */
129 /* INDXQ (output) INTEGER array, dimension (N) */
130 /* The permutation which will reintegrate the subproblem just */
131 /* solved back into sorted order, i.e., D( INDXQ( I = 1, N ) ) */
132 /* will be in ascending order. */
134 /* RHO (input) REAL */
135 /* The subdiagonal element used to create the rank-1 */
138 /* CUTPNT (input) INTEGER */
139 /* Contains the location of the last eigenvalue in the leading */
140 /* sub-matrix. min(1,N) <= CUTPNT <= N. */
142 /* QSTORE (input/output) REAL array, dimension (N**2+1) */
143 /* Stores eigenvectors of submatrices encountered during */
144 /* divide and conquer, packed together. QPTR points to */
145 /* beginning of the submatrices. */
147 /* QPTR (input/output) INTEGER array, dimension (N+2) */
148 /* List of indices pointing to beginning of submatrices stored */
149 /* in QSTORE. The submatrices are numbered starting at the */
150 /* bottom left of the divide and conquer tree, from left to */
151 /* right and bottom to top. */
153 /* PRMPTR (input) INTEGER array, dimension (N lg N) */
154 /* Contains a list of pointers which indicate where in PERM a */
155 /* level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) */
156 /* indicates the size of the permutation and also the size of */
157 /* the full, non-deflated problem. */
159 /* PERM (input) INTEGER array, dimension (N lg N) */
160 /* Contains the permutations (from deflation and sorting) to be */
161 /* applied to each eigenblock. */
163 /* GIVPTR (input) INTEGER array, dimension (N lg N) */
164 /* Contains a list of pointers which indicate where in GIVCOL a */
165 /* level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) */
166 /* indicates the number of Givens rotations. */
168 /* GIVCOL (input) INTEGER array, dimension (2, N lg N) */
169 /* Each pair of numbers indicates a pair of columns to take place */
170 /* in a Givens rotation. */
172 /* GIVNUM (input) REAL array, dimension (2, N lg N) */
173 /* Each number indicates the S value to be used in the */
174 /* corresponding Givens rotation. */
176 /* WORK (workspace) REAL array, dimension (3*N+QSIZ*N) */
178 /* IWORK (workspace) INTEGER array, dimension (4*N) */
180 /* INFO (output) INTEGER */
181 /* = 0: successful exit. */
182 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
183 /* > 0: if INFO = 1, an eigenvalue did not converge */
185 /* Further Details */
186 /* =============== */
188 /* Based on contributions by */
189 /* Jeff Rutter, Computer Science Division, University of California */
190 /* at Berkeley, USA */
192 /* ===================================================================== */
194 /* .. Parameters .. */
196 /* .. Local Scalars .. */
198 /* .. External Subroutines .. */
200 /* .. Intrinsic Functions .. */
202 /* .. Executable Statements .. */
204 /* Test the input parameters. */
206 /* Parameter adjustments */
209 q_offset = 1 + q_dim1;
225 if (*icompq < 0 || *icompq > 1) {
229 } else if (*icompq == 1 && *qsiz < *n) {
231 } else if (*ldq < max(1,*n)) {
233 } else if (min(1,*n) > *cutpnt || *n < *cutpnt) {
238 xerbla_("SLAED7", &i__1);
242 /* Quick return if possible */
248 /* The following values are for bookkeeping purposes only. They are */
249 /* integer pointers which indicate the portion of the workspace */
250 /* used by a particular array in SLAED8 and SLAED9. */
262 is = iq2 + *n * ldq2;
269 /* Form the z-vector which consists of the last row of Q_1 and the */
270 /* first row of Q_2. */
272 ptr = pow_ii(&c__2, tlvls) + 1;
274 for (i__ = 1; i__ <= i__1; ++i__) {
276 ptr += pow_ii(&c__2, &i__2);
279 curr = ptr + *curpbm;
280 slaeda_(n, tlvls, curlvl, curpbm, &prmptr[1], &perm[1], &givptr[1], &
281 givcol[3], &givnum[3], &qstore[1], &qptr[1], &work[iz], &work[iz
284 /* When solving the final problem, we no longer need the stored data, */
285 /* so we will overwrite the data from this level onto the previously */
286 /* used storage space. */
288 if (*curlvl == *tlvls) {
294 /* Sort and Deflate eigenvalues. */
296 slaed8_(icompq, &k, n, qsiz, &d__[1], &q[q_offset], ldq, &indxq[1], rho,
297 cutpnt, &work[iz], &work[idlmda], &work[iq2], &ldq2, &work[iw], &
298 perm[prmptr[curr]], &givptr[curr + 1], &givcol[(givptr[curr] << 1)
299 + 1], &givnum[(givptr[curr] << 1) + 1], &iwork[indxp], &iwork[
301 prmptr[curr + 1] = prmptr[curr] + *n;
302 givptr[curr + 1] += givptr[curr];
304 /* Solve Secular Equation. */
307 slaed9_(&k, &c__1, &k, n, &d__[1], &work[is], &k, rho, &work[idlmda],
308 &work[iw], &qstore[qptr[curr]], &k, info);
313 sgemm_("N", "N", qsiz, &k, &k, &c_b10, &work[iq2], &ldq2, &qstore[
314 qptr[curr]], &k, &c_b11, &q[q_offset], ldq);
316 /* Computing 2nd power */
318 qptr[curr + 1] = qptr[curr] + i__1 * i__1;
320 /* Prepare the INDXQ sorting permutation. */
324 slamrg_(&n1, &n2, &d__[1], &c__1, &c_n1, &indxq[1]);
326 qptr[curr + 1] = qptr[curr];
328 for (i__ = 1; i__ <= i__1; ++i__) {