001:       SUBROUTINE SLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
002:      $                   LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR,
003:      $                   PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK,
004:      $                   INFO )
005: *
006: *  -- LAPACK routine (version 3.2) --
007: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
008: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
009: *     November 2006
010: *
011: *     .. Scalar Arguments ..
012:       INTEGER            CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
013:      $                   QSIZ, TLVLS
014:       REAL               RHO
015: *     ..
016: *     .. Array Arguments ..
017:       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
018:      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
019:       REAL               D( * ), GIVNUM( 2, * ), Q( LDQ, * ),
020:      $                   QSTORE( * ), WORK( * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  SLAED7 computes the updated eigensystem of a diagonal
027: *  matrix after modification by a rank-one symmetric matrix. This
028: *  routine is used only for the eigenproblem which requires all
029: *  eigenvalues and optionally eigenvectors of a dense symmetric matrix
030: *  that has been reduced to tridiagonal form.  SLAED1 handles
031: *  the case in which all eigenvalues and eigenvectors of a symmetric
032: *  tridiagonal matrix are desired.
033: *
034: *    T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
035: *
036: *     where Z = Q'u, u is a vector of length N with ones in the
037: *     CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
038: *
039: *     The eigenvectors of the original matrix are stored in Q, and the
040: *     eigenvalues are in D.  The algorithm consists of three stages:
041: *
042: *        The first stage consists of deflating the size of the problem
043: *        when there are multiple eigenvalues or if there is a zero in
044: *        the Z vector.  For each such occurence the dimension of the
045: *        secular equation problem is reduced by one.  This stage is
046: *        performed by the routine SLAED8.
047: *
048: *        The second stage consists of calculating the updated
049: *        eigenvalues. This is done by finding the roots of the secular
050: *        equation via the routine SLAED4 (as called by SLAED9).
051: *        This routine also calculates the eigenvectors of the current
052: *        problem.
053: *
054: *        The final stage consists of computing the updated eigenvectors
055: *        directly using the updated eigenvalues.  The eigenvectors for
056: *        the current problem are multiplied with the eigenvectors from
057: *        the overall problem.
058: *
059: *  Arguments
060: *  =========
061: *
062: *  ICOMPQ  (input) INTEGER
063: *          = 0:  Compute eigenvalues only.
064: *          = 1:  Compute eigenvectors of original dense symmetric matrix
065: *                also.  On entry, Q contains the orthogonal matrix used
066: *                to reduce the original matrix to tridiagonal form.
067: *
068: *  N      (input) INTEGER
069: *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
070: *
071: *  QSIZ   (input) INTEGER
072: *         The dimension of the orthogonal matrix used to reduce
073: *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
074: *
075: *  TLVLS  (input) INTEGER
076: *         The total number of merging levels in the overall divide and
077: *         conquer tree.
078: *
079: *  CURLVL (input) INTEGER
080: *         The current level in the overall merge routine,
081: *         0 <= CURLVL <= TLVLS.
082: *
083: *  CURPBM (input) INTEGER
084: *         The current problem in the current level in the overall
085: *         merge routine (counting from upper left to lower right).
086: *
087: *  D      (input/output) REAL array, dimension (N)
088: *         On entry, the eigenvalues of the rank-1-perturbed matrix.
089: *         On exit, the eigenvalues of the repaired matrix.
090: *
091: *  Q      (input/output) REAL array, dimension (LDQ, N)
092: *         On entry, the eigenvectors of the rank-1-perturbed matrix.
093: *         On exit, the eigenvectors of the repaired tridiagonal matrix.
094: *
095: *  LDQ    (input) INTEGER
096: *         The leading dimension of the array Q.  LDQ >= max(1,N).
097: *
098: *  INDXQ  (output) INTEGER array, dimension (N)
099: *         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) REAL
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) REAL 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) REAL 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) REAL 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:       REAL               ONE, ZERO
165:       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
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           SGEMM, SLAED8, SLAED9, SLAEDA, SLAMRG, XERBLA
173: *     ..
174: *     .. Intrinsic Functions ..
175:       INTRINSIC          MAX, MIN
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.MAX( 1, N ) ) THEN
190:          INFO = -9
191:       ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
192:          INFO = -12
193:       END IF
194:       IF( INFO.NE.0 ) THEN
195:          CALL XERBLA( 'SLAED7', -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 SLAED8 and SLAED9.
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 SLAEDA( 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 SLAED8( 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 SLAED9( 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 SGEMM( '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 SLAMRG( 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 SLAED7
287: *
288:       END
289: