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