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