001:       SUBROUTINE SLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
002:      $                   CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
003:      $                   GIVCOL, GIVNUM, INDXP, INDX, INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
011:      $                   QSIZ
012:       REAL               RHO
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            GIVCOL( 2, * ), INDX( * ), INDXP( * ),
016:      $                   INDXQ( * ), PERM( * )
017:       REAL               D( * ), DLAMDA( * ), GIVNUM( 2, * ),
018:      $                   Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  SLAED8 merges the two sets of eigenvalues together into a single
025: *  sorted set.  Then it tries to deflate the size of the problem.
026: *  There are two ways in which deflation can occur:  when two or more
027: *  eigenvalues are close together or if there is a tiny element in the
028: *  Z vector.  For each such occurrence the order of the related secular
029: *  equation problem is reduced by one.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  ICOMPQ  (input) INTEGER
035: *          = 0:  Compute eigenvalues only.
036: *          = 1:  Compute eigenvectors of original dense symmetric matrix
037: *                also.  On entry, Q contains the orthogonal matrix used
038: *                to reduce the original matrix to tridiagonal form.
039: *
040: *  K      (output) INTEGER
041: *         The number of non-deflated eigenvalues, and the order of the
042: *         related secular equation.
043: *
044: *  N      (input) INTEGER
045: *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
046: *
047: *  QSIZ   (input) INTEGER
048: *         The dimension of the orthogonal matrix used to reduce
049: *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
050: *
051: *  D      (input/output) REAL array, dimension (N)
052: *         On entry, the eigenvalues of the two submatrices to be
053: *         combined.  On exit, the trailing (N-K) updated eigenvalues
054: *         (those which were deflated) sorted into increasing order.
055: *
056: *  Q      (input/output) REAL array, dimension (LDQ,N)
057: *         If ICOMPQ = 0, Q is not referenced.  Otherwise,
058: *         on entry, Q contains the eigenvectors of the partially solved
059: *         system which has been previously updated in matrix
060: *         multiplies with other partially solved eigensystems.
061: *         On exit, Q contains the trailing (N-K) updated eigenvectors
062: *         (those which were deflated) in its last N-K columns.
063: *
064: *  LDQ    (input) INTEGER
065: *         The leading dimension of the array Q.  LDQ >= max(1,N).
066: *
067: *  INDXQ  (input) INTEGER array, dimension (N)
068: *         The permutation which separately sorts the two sub-problems
069: *         in D into ascending order.  Note that elements in the second
070: *         half of this permutation must first have CUTPNT added to
071: *         their values in order to be accurate.
072: *
073: *  RHO    (input/output) REAL
074: *         On entry, the off-diagonal element associated with the rank-1
075: *         cut which originally split the two submatrices which are now
076: *         being recombined.
077: *         On exit, RHO has been modified to the value required by
078: *         SLAED3.
079: *
080: *  CUTPNT (input) INTEGER
081: *         The location of the last eigenvalue in the leading
082: *         sub-matrix.  min(1,N) <= CUTPNT <= N.
083: *
084: *  Z      (input) REAL array, dimension (N)
085: *         On entry, Z contains the updating vector (the last row of
086: *         the first sub-eigenvector matrix and the first row of the
087: *         second sub-eigenvector matrix).
088: *         On exit, the contents of Z are destroyed by the updating
089: *         process.
090: *
091: *  DLAMDA (output) REAL array, dimension (N)
092: *         A copy of the first K eigenvalues which will be used by
093: *         SLAED3 to form the secular equation.
094: *
095: *  Q2     (output) REAL array, dimension (LDQ2,N)
096: *         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
097: *         a copy of the first K eigenvectors which will be used by
098: *         SLAED7 in a matrix multiply (SGEMM) to update the new
099: *         eigenvectors.
100: *
101: *  LDQ2   (input) INTEGER
102: *         The leading dimension of the array Q2.  LDQ2 >= max(1,N).
103: *
104: *  W      (output) REAL array, dimension (N)
105: *         The first k values of the final deflation-altered z-vector and
106: *         will be passed to SLAED3.
107: *
108: *  PERM   (output) INTEGER array, dimension (N)
109: *         The permutations (from deflation and sorting) to be applied
110: *         to each eigenblock.
111: *
112: *  GIVPTR (output) INTEGER
113: *         The number of Givens rotations which took place in this
114: *         subproblem.
115: *
116: *  GIVCOL (output) INTEGER array, dimension (2, N)
117: *         Each pair of numbers indicates a pair of columns to take place
118: *         in a Givens rotation.
119: *
120: *  GIVNUM (output) REAL array, dimension (2, N)
121: *         Each number indicates the S value to be used in the
122: *         corresponding Givens rotation.
123: *
124: *  INDXP  (workspace) INTEGER array, dimension (N)
125: *         The permutation used to place deflated values of D at the end
126: *         of the array.  INDXP(1:K) points to the nondeflated D-values
127: *         and INDXP(K+1:N) points to the deflated eigenvalues.
128: *
129: *  INDX   (workspace) INTEGER array, dimension (N)
130: *         The permutation used to sort the contents of D into ascending
131: *         order.
132: *
133: *  INFO   (output) INTEGER
134: *          = 0:  successful exit.
135: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
136: *
137: *  Further Details
138: *  ===============
139: *
140: *  Based on contributions by
141: *     Jeff Rutter, Computer Science Division, University of California
142: *     at Berkeley, USA
143: *
144: *  =====================================================================
145: *
146: *     .. Parameters ..
147:       REAL               MONE, ZERO, ONE, TWO, EIGHT
148:       PARAMETER          ( MONE = -1.0E0, ZERO = 0.0E0, ONE = 1.0E0,
149:      $                   TWO = 2.0E0, EIGHT = 8.0E0 )
150: *     ..
151: *     .. Local Scalars ..
152: *
153:       INTEGER            I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
154:       REAL               C, EPS, S, T, TAU, TOL
155: *     ..
156: *     .. External Functions ..
157:       INTEGER            ISAMAX
158:       REAL               SLAMCH, SLAPY2
159:       EXTERNAL           ISAMAX, SLAMCH, SLAPY2
160: *     ..
161: *     .. External Subroutines ..
162:       EXTERNAL           SCOPY, SLACPY, SLAMRG, SROT, SSCAL, XERBLA
163: *     ..
164: *     .. Intrinsic Functions ..
165:       INTRINSIC          ABS, MAX, MIN, SQRT
166: *     ..
167: *     .. Executable Statements ..
168: *
169: *     Test the input parameters.
170: *
171:       INFO = 0
172: *
173:       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
174:          INFO = -1
175:       ELSE IF( N.LT.0 ) THEN
176:          INFO = -3
177:       ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
178:          INFO = -4
179:       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
180:          INFO = -7
181:       ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
182:          INFO = -10
183:       ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
184:          INFO = -14
185:       END IF
186:       IF( INFO.NE.0 ) THEN
187:          CALL XERBLA( 'SLAED8', -INFO )
188:          RETURN
189:       END IF
190: *
191: *     Quick return if possible
192: *
193:       IF( N.EQ.0 )
194:      $   RETURN
195: *
196:       N1 = CUTPNT
197:       N2 = N - N1
198:       N1P1 = N1 + 1
199: *
200:       IF( RHO.LT.ZERO ) THEN
201:          CALL SSCAL( N2, MONE, Z( N1P1 ), 1 )
202:       END IF
203: *
204: *     Normalize z so that norm(z) = 1
205: *
206:       T = ONE / SQRT( TWO )
207:       DO 10 J = 1, N
208:          INDX( J ) = J
209:    10 CONTINUE
210:       CALL SSCAL( N, T, Z, 1 )
211:       RHO = ABS( TWO*RHO )
212: *
213: *     Sort the eigenvalues into increasing order
214: *
215:       DO 20 I = CUTPNT + 1, N
216:          INDXQ( I ) = INDXQ( I ) + CUTPNT
217:    20 CONTINUE
218:       DO 30 I = 1, N
219:          DLAMDA( I ) = D( INDXQ( I ) )
220:          W( I ) = Z( INDXQ( I ) )
221:    30 CONTINUE
222:       I = 1
223:       J = CUTPNT + 1
224:       CALL SLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
225:       DO 40 I = 1, N
226:          D( I ) = DLAMDA( INDX( I ) )
227:          Z( I ) = W( INDX( I ) )
228:    40 CONTINUE
229: *
230: *     Calculate the allowable deflation tolerence
231: *
232:       IMAX = ISAMAX( N, Z, 1 )
233:       JMAX = ISAMAX( N, D, 1 )
234:       EPS = SLAMCH( 'Epsilon' )
235:       TOL = EIGHT*EPS*ABS( D( JMAX ) )
236: *
237: *     If the rank-1 modifier is small enough, no more needs to be done
238: *     except to reorganize Q so that its columns correspond with the
239: *     elements in D.
240: *
241:       IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
242:          K = 0
243:          IF( ICOMPQ.EQ.0 ) THEN
244:             DO 50 J = 1, N
245:                PERM( J ) = INDXQ( INDX( J ) )
246:    50       CONTINUE
247:          ELSE
248:             DO 60 J = 1, N
249:                PERM( J ) = INDXQ( INDX( J ) )
250:                CALL SCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
251:    60       CONTINUE
252:             CALL SLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ),
253:      $                   LDQ )
254:          END IF
255:          RETURN
256:       END IF
257: *
258: *     If there are multiple eigenvalues then the problem deflates.  Here
259: *     the number of equal eigenvalues are found.  As each equal
260: *     eigenvalue is found, an elementary reflector is computed to rotate
261: *     the corresponding eigensubspace so that the corresponding
262: *     components of Z are zero in this new basis.
263: *
264:       K = 0
265:       GIVPTR = 0
266:       K2 = N + 1
267:       DO 70 J = 1, N
268:          IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
269: *
270: *           Deflate due to small z component.
271: *
272:             K2 = K2 - 1
273:             INDXP( K2 ) = J
274:             IF( J.EQ.N )
275:      $         GO TO 110
276:          ELSE
277:             JLAM = J
278:             GO TO 80
279:          END IF
280:    70 CONTINUE
281:    80 CONTINUE
282:       J = J + 1
283:       IF( J.GT.N )
284:      $   GO TO 100
285:       IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
286: *
287: *        Deflate due to small z component.
288: *
289:          K2 = K2 - 1
290:          INDXP( K2 ) = J
291:       ELSE
292: *
293: *        Check if eigenvalues are close enough to allow deflation.
294: *
295:          S = Z( JLAM )
296:          C = Z( J )
297: *
298: *        Find sqrt(a**2+b**2) without overflow or
299: *        destructive underflow.
300: *
301:          TAU = SLAPY2( C, S )
302:          T = D( J ) - D( JLAM )
303:          C = C / TAU
304:          S = -S / TAU
305:          IF( ABS( T*C*S ).LE.TOL ) THEN
306: *
307: *           Deflation is possible.
308: *
309:             Z( J ) = TAU
310:             Z( JLAM ) = ZERO
311: *
312: *           Record the appropriate Givens rotation
313: *
314:             GIVPTR = GIVPTR + 1
315:             GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
316:             GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
317:             GIVNUM( 1, GIVPTR ) = C
318:             GIVNUM( 2, GIVPTR ) = S
319:             IF( ICOMPQ.EQ.1 ) THEN
320:                CALL SROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
321:      $                    Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
322:             END IF
323:             T = D( JLAM )*C*C + D( J )*S*S
324:             D( J ) = D( JLAM )*S*S + D( J )*C*C
325:             D( JLAM ) = T
326:             K2 = K2 - 1
327:             I = 1
328:    90       CONTINUE
329:             IF( K2+I.LE.N ) THEN
330:                IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
331:                   INDXP( K2+I-1 ) = INDXP( K2+I )
332:                   INDXP( K2+I ) = JLAM
333:                   I = I + 1
334:                   GO TO 90
335:                ELSE
336:                   INDXP( K2+I-1 ) = JLAM
337:                END IF
338:             ELSE
339:                INDXP( K2+I-1 ) = JLAM
340:             END IF
341:             JLAM = J
342:          ELSE
343:             K = K + 1
344:             W( K ) = Z( JLAM )
345:             DLAMDA( K ) = D( JLAM )
346:             INDXP( K ) = JLAM
347:             JLAM = J
348:          END IF
349:       END IF
350:       GO TO 80
351:   100 CONTINUE
352: *
353: *     Record the last eigenvalue.
354: *
355:       K = K + 1
356:       W( K ) = Z( JLAM )
357:       DLAMDA( K ) = D( JLAM )
358:       INDXP( K ) = JLAM
359: *
360:   110 CONTINUE
361: *
362: *     Sort the eigenvalues and corresponding eigenvectors into DLAMDA
363: *     and Q2 respectively.  The eigenvalues/vectors which were not
364: *     deflated go into the first K slots of DLAMDA and Q2 respectively,
365: *     while those which were deflated go into the last N - K slots.
366: *
367:       IF( ICOMPQ.EQ.0 ) THEN
368:          DO 120 J = 1, N
369:             JP = INDXP( J )
370:             DLAMDA( J ) = D( JP )
371:             PERM( J ) = INDXQ( INDX( JP ) )
372:   120    CONTINUE
373:       ELSE
374:          DO 130 J = 1, N
375:             JP = INDXP( J )
376:             DLAMDA( J ) = D( JP )
377:             PERM( J ) = INDXQ( INDX( JP ) )
378:             CALL SCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
379:   130    CONTINUE
380:       END IF
381: *
382: *     The deflated eigenvalues and their corresponding vectors go back
383: *     into the last N - K slots of D and Q respectively.
384: *
385:       IF( K.LT.N ) THEN
386:          IF( ICOMPQ.EQ.0 ) THEN
387:             CALL SCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
388:          ELSE
389:             CALL SCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
390:             CALL SLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2,
391:      $                   Q( 1, K+1 ), LDQ )
392:          END IF
393:       END IF
394: *
395:       RETURN
396: *
397: *     End of SLAED8
398: *
399:       END
400: