1       SUBROUTINE DLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
  2      $                   LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR,
  3      $                   PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK,
  4      $                   INFO )
  5 *
  6 *  -- LAPACK routine (version 3.3.1) --
  7 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  8 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  9 *  -- April 2011                                                      --
 10 *
 11 *     .. Scalar Arguments ..
 12       INTEGER            CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
 13      $                   QSIZ, TLVLS
 14       DOUBLE PRECISION   RHO
 15 *     ..
 16 *     .. Array Arguments ..
 17       INTEGER            GIVCOL( 2* ), GIVPTR( * ), INDXQ( * ),
 18      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
 19       DOUBLE PRECISION   D( * ), GIVNUM( 2* ), Q( LDQ, * ),
 20      $                   QSTORE( * ), WORK( * )
 21 *     ..
 22 *
 23 *  Purpose
 24 *  =======
 25 *
 26 *  DLAED7 computes the updated eigensystem of a diagonal
 27 *  matrix after modification by a rank-one symmetric matrix. This
 28 *  routine is used only for the eigenproblem which requires all
 29 *  eigenvalues and optionally eigenvectors of a dense symmetric matrix
 30 *  that has been reduced to tridiagonal form.  DLAED1 handles
 31 *  the case in which all eigenvalues and eigenvectors of a symmetric
 32 *  tridiagonal matrix are desired.
 33 *
 34 *    T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
 35 *
 36 *     where Z = Q**Tu, u is a vector of length N with ones in the
 37 *     CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
 38 *
 39 *     The eigenvectors of the original matrix are stored in Q, and the
 40 *     eigenvalues are in D.  The algorithm consists of three stages:
 41 *
 42 *        The first stage consists of deflating the size of the problem
 43 *        when there are multiple eigenvalues or if there is a zero in
 44 *        the Z vector.  For each such occurence the dimension of the
 45 *        secular equation problem is reduced by one.  This stage is
 46 *        performed by the routine DLAED8.
 47 *
 48 *        The second stage consists of calculating the updated
 49 *        eigenvalues. This is done by finding the roots of the secular
 50 *        equation via the routine DLAED4 (as called by DLAED9).
 51 *        This routine also calculates the eigenvectors of the current
 52 *        problem.
 53 *
 54 *        The final stage consists of computing the updated eigenvectors
 55 *        directly using the updated eigenvalues.  The eigenvectors for
 56 *        the current problem are multiplied with the eigenvectors from
 57 *        the overall problem.
 58 *
 59 *  Arguments
 60 *  =========
 61 *
 62 *  ICOMPQ  (input) INTEGER
 63 *          = 0:  Compute eigenvalues only.
 64 *          = 1:  Compute eigenvectors of original dense symmetric matrix
 65 *                also.  On entry, Q contains the orthogonal matrix used
 66 *                to reduce the original matrix to tridiagonal form.
 67 *
 68 *  N      (input) INTEGER
 69 *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
 70 *
 71 *  QSIZ   (input) INTEGER
 72 *         The dimension of the orthogonal matrix used to reduce
 73 *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
 74 *
 75 *  TLVLS  (input) INTEGER
 76 *         The total number of merging levels in the overall divide and
 77 *         conquer tree.
 78 *
 79 *  CURLVL (input) INTEGER
 80 *         The current level in the overall merge routine,
 81 *         0 <= CURLVL <= TLVLS.
 82 *
 83 *  CURPBM (input) INTEGER
 84 *         The current problem in the current level in the overall
 85 *         merge routine (counting from upper left to lower right).
 86 *
 87 *  D      (input/output) DOUBLE PRECISION array, dimension (N)
 88 *         On entry, the eigenvalues of the rank-1-perturbed matrix.
 89 *         On exit, the eigenvalues of the repaired matrix.
 90 *
 91 *  Q      (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
 92 *         On entry, the eigenvectors of the rank-1-perturbed matrix.
 93 *         On exit, the eigenvectors of the repaired tridiagonal matrix.
 94 *
 95 *  LDQ    (input) INTEGER
 96 *         The leading dimension of the array Q.  LDQ >= max(1,N).
 97 *
 98 *  INDXQ  (output) INTEGER array, dimension (N)
 99 *         The permutation which will reintegrate the subproblem just
100 *         solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
101 *         will be in ascending order.
102 *
103 *  RHO    (input) DOUBLE PRECISION
104 *         The subdiagonal element used to create the rank-1
105 *         modification.
106 *
107 *  CUTPNT (input) INTEGER
108 *         Contains the location of the last eigenvalue in the leading
109 *         sub-matrix.  min(1,N) <= CUTPNT <= N.
110 *
111 *  QSTORE (input/output) DOUBLE PRECISION array, dimension (N**2+1)
112 *         Stores eigenvectors of submatrices encountered during
113 *         divide and conquer, packed together. QPTR points to
114 *         beginning of the submatrices.
115 *
116 *  QPTR   (input/output) INTEGER array, dimension (N+2)
117 *         List of indices pointing to beginning of submatrices stored
118 *         in QSTORE. The submatrices are numbered starting at the
119 *         bottom left of the divide and conquer tree, from left to
120 *         right and bottom to top.
121 *
122 *  PRMPTR (input) INTEGER array, dimension (N lg N)
123 *         Contains a list of pointers which indicate where in PERM a
124 *         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
125 *         indicates the size of the permutation and also the size of
126 *         the full, non-deflated problem.
127 *
128 *  PERM   (input) INTEGER array, dimension (N lg N)
129 *         Contains the permutations (from deflation and sorting) to be
130 *         applied to each eigenblock.
131 *
132 *  GIVPTR (input) INTEGER array, dimension (N lg N)
133 *         Contains a list of pointers which indicate where in GIVCOL a
134 *         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
135 *         indicates the number of Givens rotations.
136 *
137 *  GIVCOL (input) INTEGER array, dimension (2, N lg N)
138 *         Each pair of numbers indicates a pair of columns to take place
139 *         in a Givens rotation.
140 *
141 *  GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
142 *         Each number indicates the S value to be used in the
143 *         corresponding Givens rotation.
144 *
145 *  WORK   (workspace) DOUBLE PRECISION array, dimension (3*N+QSIZ*N)
146 *
147 *  IWORK  (workspace) INTEGER array, dimension (4*N)
148 *
149 *  INFO   (output) INTEGER
150 *          = 0:  successful exit.
151 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
152 *          > 0:  if INFO = 1, an eigenvalue did not converge
153 *
154 *  Further Details
155 *  ===============
156 *
157 *  Based on contributions by
158 *     Jeff Rutter, Computer Science Division, University of California
159 *     at Berkeley, USA
160 *
161 *  =====================================================================
162 *
163 *     .. Parameters ..
164       DOUBLE PRECISION   ONE, ZERO
165       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
166 *     ..
167 *     .. Local Scalars ..
168       INTEGER            COLTYP, CURR, I, IDLMDA, INDX, INDXC, INDXP,
169      $                   IQ2, IS, IW, IZ, K, LDQ2, N1, N2, PTR
170 *     ..
171 *     .. External Subroutines ..
172       EXTERNAL           DGEMM, DLAED8, DLAED9, DLAEDA, DLAMRG, XERBLA
173 *     ..
174 *     .. Intrinsic Functions ..
175       INTRINSIC          MAXMIN
176 *     ..
177 *     .. Executable Statements ..
178 *
179 *     Test the input parameters.
180 *
181       INFO = 0
182 *
183       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
184          INFO = -1
185       ELSE IF( N.LT.0 ) THEN
186          INFO = -2
187       ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
188          INFO = -4
189       ELSE IF( LDQ.LT.MAX1, N ) ) THEN
190          INFO = -9
191       ELSE IFMIN1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
192          INFO = -12
193       END IF
194       IF( INFO.NE.0 ) THEN
195          CALL XERBLA( 'DLAED7'-INFO )
196          RETURN
197       END IF
198 *
199 *     Quick return if possible
200 *
201       IF( N.EQ.0 )
202      $   RETURN
203 *
204 *     The following values are for bookkeeping purposes only.  They are
205 *     integer pointers which indicate the portion of the workspace
206 *     used by a particular array in DLAED8 and DLAED9.
207 *
208       IF( ICOMPQ.EQ.1 ) THEN
209          LDQ2 = QSIZ
210       ELSE
211          LDQ2 = N
212       END IF
213 *
214       IZ = 1
215       IDLMDA = IZ + N
216       IW = IDLMDA + N
217       IQ2 = IW + N
218       IS = IQ2 + N*LDQ2
219 *
220       INDX = 1
221       INDXC = INDX + N
222       COLTYP = INDXC + N
223       INDXP = COLTYP + N
224 *
225 *     Form the z-vector which consists of the last row of Q_1 and the
226 *     first row of Q_2.
227 *
228       PTR = 1 + 2**TLVLS
229       DO 10 I = 1, CURLVL - 1
230          PTR = PTR + 2**( TLVLS-I )
231    10 CONTINUE
232       CURR = PTR + CURPBM
233       CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
234      $             GIVCOL, GIVNUM, QSTORE, QPTR, WORK( IZ ),
235      $             WORK( IZ+N ), INFO )
236 *
237 *     When solving the final problem, we no longer need the stored data,
238 *     so we will overwrite the data from this level onto the previously
239 *     used storage space.
240 *
241       IF( CURLVL.EQ.TLVLS ) THEN
242          QPTR( CURR ) = 1
243          PRMPTR( CURR ) = 1
244          GIVPTR( CURR ) = 1
245       END IF
246 *
247 *     Sort and Deflate eigenvalues.
248 *
249       CALL DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT,
250      $             WORK( IZ ), WORK( IDLMDA ), WORK( IQ2 ), LDQ2,
251      $             WORK( IW ), PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
252      $             GIVCOL( 1, GIVPTR( CURR ) ),
253      $             GIVNUM( 1, GIVPTR( CURR ) ), IWORK( INDXP ),
254      $             IWORK( INDX ), INFO )
255       PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
256       GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
257 *
258 *     Solve Secular Equation.
259 *
260       IF( K.NE.0 ) THEN
261          CALL DLAED9( K, 1, K, N, D, WORK( IS ), K, RHO, WORK( IDLMDA ),
262      $                WORK( IW ), QSTORE( QPTR( CURR ) ), K, INFO )
263          IF( INFO.NE.0 )
264      $      GO TO 30
265          IF( ICOMPQ.EQ.1 ) THEN
266             CALL DGEMM( 'N''N', QSIZ, K, K, ONE, WORK( IQ2 ), LDQ2,
267      $                  QSTORE( QPTR( CURR ) ), K, ZERO, Q, LDQ )
268          END IF
269          QPTR( CURR+1 ) = QPTR( CURR ) + K**2
270 *
271 *     Prepare the INDXQ sorting permutation.
272 *
273          N1 = K
274          N2 = N - K
275          CALL DLAMRG( N1, N2, D, 1-1, INDXQ )
276       ELSE
277          QPTR( CURR+1 ) = QPTR( CURR )
278          DO 20 I = 1, N
279             INDXQ( I ) = I
280    20    CONTINUE
281       END IF
282 *
283    30 CONTINUE
284       RETURN
285 *
286 *     End of DLAED7
287 *
288       END