001:       SUBROUTINE SLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W,
002:      $                   Q2, INDX, INDXC, INDXP, COLTYP, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            INFO, K, LDQ, N, N1
011:       REAL               RHO
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ),
015:      $                   INDXQ( * )
016:       REAL               D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
017:      $                   W( * ), Z( * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  SLAED2 merges the two sets of eigenvalues together into a single
024: *  sorted set.  Then it tries to deflate the size of the problem.
025: *  There are two ways in which deflation can occur:  when two or more
026: *  eigenvalues are close together or if there is a tiny entry in the
027: *  Z vector.  For each such occurrence the order of the related secular
028: *  equation problem is reduced by one.
029: *
030: *  Arguments
031: *  =========
032: *
033: *  K      (output) INTEGER
034: *         The number of non-deflated eigenvalues, and the order of the
035: *         related secular equation. 0 <= K <=N.
036: *
037: *  N      (input) INTEGER
038: *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
039: *
040: *  N1     (input) INTEGER
041: *         The location of the last eigenvalue in the leading sub-matrix.
042: *         min(1,N) <= N1 <= N/2.
043: *
044: *  D      (input/output) REAL array, dimension (N)
045: *         On entry, D contains the eigenvalues of the two submatrices to
046: *         be combined.
047: *         On exit, D contains the trailing (N-K) updated eigenvalues
048: *         (those which were deflated) sorted into increasing order.
049: *
050: *  Q      (input/output) REAL array, dimension (LDQ, N)
051: *         On entry, Q contains the eigenvectors of two submatrices in
052: *         the two square blocks with corners at (1,1), (N1,N1)
053: *         and (N1+1, N1+1), (N,N).
054: *         On exit, Q contains the trailing (N-K) updated eigenvectors
055: *         (those which were deflated) in its last N-K columns.
056: *
057: *  LDQ    (input) INTEGER
058: *         The leading dimension of the array Q.  LDQ >= max(1,N).
059: *
060: *  INDXQ  (input/output) INTEGER array, dimension (N)
061: *         The permutation which separately sorts the two sub-problems
062: *         in D into ascending order.  Note that elements in the second
063: *         half of this permutation must first have N1 added to their
064: *         values. Destroyed on exit.
065: *
066: *  RHO    (input/output) REAL
067: *         On entry, the off-diagonal element associated with the rank-1
068: *         cut which originally split the two submatrices which are now
069: *         being recombined.
070: *         On exit, RHO has been modified to the value required by
071: *         SLAED3.
072: *
073: *  Z      (input) REAL array, dimension (N)
074: *         On entry, Z 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).
077: *         On exit, the contents of Z have been destroyed by the updating
078: *         process.
079: *
080: *  DLAMDA (output) REAL array, dimension (N)
081: *         A copy of the first K eigenvalues which will be used by
082: *         SLAED3 to form the secular equation.
083: *
084: *  W      (output) REAL array, dimension (N)
085: *         The first k values of the final deflation-altered z-vector
086: *         which will be passed to SLAED3.
087: *
088: *  Q2     (output) REAL array, dimension (N1**2+(N-N1)**2)
089: *         A copy of the first K eigenvectors which will be used by
090: *         SLAED3 in a matrix multiply (SGEMM) to solve for the new
091: *         eigenvectors.
092: *
093: *  INDX   (workspace) INTEGER array, dimension (N)
094: *         The permutation used to sort the contents of DLAMDA into
095: *         ascending order.
096: *
097: *  INDXC  (output) INTEGER array, dimension (N)
098: *         The permutation used to arrange the columns of the deflated
099: *         Q matrix into three groups:  the first group contains non-zero
100: *         elements only at and above N1, the second contains
101: *         non-zero elements only below N1, and the third is dense.
102: *
103: *  INDXP  (workspace) INTEGER array, dimension (N)
104: *         The permutation used to place deflated values of D at the end
105: *         of the array.  INDXP(1:K) points to the nondeflated D-values
106: *         and INDXP(K+1:N) points to the deflated eigenvalues.
107: *
108: *  COLTYP (workspace/output) INTEGER array, dimension (N)
109: *         During execution, a label which will indicate which of the
110: *         following types a column in the Q2 matrix is:
111: *         1 : non-zero in the upper half only;
112: *         2 : dense;
113: *         3 : non-zero in the lower half only;
114: *         4 : deflated.
115: *         On exit, COLTYP(i) is the number of columns of type i,
116: *         for i=1 to 4 only.
117: *
118: *  INFO   (output) INTEGER
119: *          = 0:  successful exit.
120: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
121: *
122: *  Further Details
123: *  ===============
124: *
125: *  Based on contributions by
126: *     Jeff Rutter, Computer Science Division, University of California
127: *     at Berkeley, USA
128: *  Modified by Francoise Tisseur, University of Tennessee.
129: *
130: *  =====================================================================
131: *
132: *     .. Parameters ..
133:       REAL               MONE, ZERO, ONE, TWO, EIGHT
134:       PARAMETER          ( MONE = -1.0E0, ZERO = 0.0E0, ONE = 1.0E0,
135:      $                   TWO = 2.0E0, EIGHT = 8.0E0 )
136: *     ..
137: *     .. Local Arrays ..
138:       INTEGER            CTOT( 4 ), PSM( 4 )
139: *     ..
140: *     .. Local Scalars ..
141:       INTEGER            CT, I, IMAX, IQ1, IQ2, J, JMAX, JS, K2, N1P1,
142:      $                   N2, NJ, PJ
143:       REAL               C, EPS, S, T, TAU, TOL
144: *     ..
145: *     .. External Functions ..
146:       INTEGER            ISAMAX
147:       REAL               SLAMCH, SLAPY2
148:       EXTERNAL           ISAMAX, SLAMCH, SLAPY2
149: *     ..
150: *     .. External Subroutines ..
151:       EXTERNAL           SCOPY, SLACPY, SLAMRG, SROT, SSCAL, XERBLA
152: *     ..
153: *     .. Intrinsic Functions ..
154:       INTRINSIC          ABS, MAX, MIN, SQRT
155: *     ..
156: *     .. Executable Statements ..
157: *
158: *     Test the input parameters.
159: *
160:       INFO = 0
161: *
162:       IF( N.LT.0 ) THEN
163:          INFO = -2
164:       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
165:          INFO = -6
166:       ELSE IF( MIN( 1, ( N / 2 ) ).GT.N1 .OR. ( N / 2 ).LT.N1 ) THEN
167:          INFO = -3
168:       END IF
169:       IF( INFO.NE.0 ) THEN
170:          CALL XERBLA( 'SLAED2', -INFO )
171:          RETURN
172:       END IF
173: *
174: *     Quick return if possible
175: *
176:       IF( N.EQ.0 )
177:      $   RETURN
178: *
179:       N2 = N - N1
180:       N1P1 = N1 + 1
181: *
182:       IF( RHO.LT.ZERO ) THEN
183:          CALL SSCAL( N2, MONE, Z( N1P1 ), 1 )
184:       END IF
185: *
186: *     Normalize z so that norm(z) = 1.  Since z is the concatenation of
187: *     two normalized vectors, norm2(z) = sqrt(2).
188: *
189:       T = ONE / SQRT( TWO )
190:       CALL SSCAL( N, T, Z, 1 )
191: *
192: *     RHO = ABS( norm(z)**2 * RHO )
193: *
194:       RHO = ABS( TWO*RHO )
195: *
196: *     Sort the eigenvalues into increasing order
197: *
198:       DO 10 I = N1P1, N
199:          INDXQ( I ) = INDXQ( I ) + N1
200:    10 CONTINUE
201: *
202: *     re-integrate the deflated parts from the last pass
203: *
204:       DO 20 I = 1, N
205:          DLAMDA( I ) = D( INDXQ( I ) )
206:    20 CONTINUE
207:       CALL SLAMRG( N1, N2, DLAMDA, 1, 1, INDXC )
208:       DO 30 I = 1, N
209:          INDX( I ) = INDXQ( INDXC( I ) )
210:    30 CONTINUE
211: *
212: *     Calculate the allowable deflation tolerance
213: *
214:       IMAX = ISAMAX( N, Z, 1 )
215:       JMAX = ISAMAX( N, D, 1 )
216:       EPS = SLAMCH( 'Epsilon' )
217:       TOL = EIGHT*EPS*MAX( ABS( D( JMAX ) ), ABS( Z( IMAX ) ) )
218: *
219: *     If the rank-1 modifier is small enough, no more needs to be done
220: *     except to reorganize Q so that its columns correspond with the
221: *     elements in D.
222: *
223:       IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
224:          K = 0
225:          IQ2 = 1
226:          DO 40 J = 1, N
227:             I = INDX( J )
228:             CALL SCOPY( N, Q( 1, I ), 1, Q2( IQ2 ), 1 )
229:             DLAMDA( J ) = D( I )
230:             IQ2 = IQ2 + N
231:    40    CONTINUE
232:          CALL SLACPY( 'A', N, N, Q2, N, Q, LDQ )
233:          CALL SCOPY( N, DLAMDA, 1, D, 1 )
234:          GO TO 190
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:       DO 50 I = 1, N1
244:          COLTYP( I ) = 1
245:    50 CONTINUE
246:       DO 60 I = N1P1, N
247:          COLTYP( I ) = 3
248:    60 CONTINUE
249: *
250: *
251:       K = 0
252:       K2 = N + 1
253:       DO 70 J = 1, N
254:          NJ = INDX( J )
255:          IF( RHO*ABS( Z( NJ ) ).LE.TOL ) THEN
256: *
257: *           Deflate due to small z component.
258: *
259:             K2 = K2 - 1
260:             COLTYP( NJ ) = 4
261:             INDXP( K2 ) = NJ
262:             IF( J.EQ.N )
263:      $         GO TO 100
264:          ELSE
265:             PJ = NJ
266:             GO TO 80
267:          END IF
268:    70 CONTINUE
269:    80 CONTINUE
270:       J = J + 1
271:       NJ = INDX( J )
272:       IF( J.GT.N )
273:      $   GO TO 100
274:       IF( RHO*ABS( Z( NJ ) ).LE.TOL ) THEN
275: *
276: *        Deflate due to small z component.
277: *
278:          K2 = K2 - 1
279:          COLTYP( NJ ) = 4
280:          INDXP( K2 ) = NJ
281:       ELSE
282: *
283: *        Check if eigenvalues are close enough to allow deflation.
284: *
285:          S = Z( PJ )
286:          C = Z( NJ )
287: *
288: *        Find sqrt(a**2+b**2) without overflow or
289: *        destructive underflow.
290: *
291:          TAU = SLAPY2( C, S )
292:          T = D( NJ ) - D( PJ )
293:          C = C / TAU
294:          S = -S / TAU
295:          IF( ABS( T*C*S ).LE.TOL ) THEN
296: *
297: *           Deflation is possible.
298: *
299:             Z( NJ ) = TAU
300:             Z( PJ ) = ZERO
301:             IF( COLTYP( NJ ).NE.COLTYP( PJ ) )
302:      $         COLTYP( NJ ) = 2
303:             COLTYP( PJ ) = 4
304:             CALL SROT( N, Q( 1, PJ ), 1, Q( 1, NJ ), 1, C, S )
305:             T = D( PJ )*C**2 + D( NJ )*S**2
306:             D( NJ ) = D( PJ )*S**2 + D( NJ )*C**2
307:             D( PJ ) = T
308:             K2 = K2 - 1
309:             I = 1
310:    90       CONTINUE
311:             IF( K2+I.LE.N ) THEN
312:                IF( D( PJ ).LT.D( INDXP( K2+I ) ) ) THEN
313:                   INDXP( K2+I-1 ) = INDXP( K2+I )
314:                   INDXP( K2+I ) = PJ
315:                   I = I + 1
316:                   GO TO 90
317:                ELSE
318:                   INDXP( K2+I-1 ) = PJ
319:                END IF
320:             ELSE
321:                INDXP( K2+I-1 ) = PJ
322:             END IF
323:             PJ = NJ
324:          ELSE
325:             K = K + 1
326:             DLAMDA( K ) = D( PJ )
327:             W( K ) = Z( PJ )
328:             INDXP( K ) = PJ
329:             PJ = NJ
330:          END IF
331:       END IF
332:       GO TO 80
333:   100 CONTINUE
334: *
335: *     Record the last eigenvalue.
336: *
337:       K = K + 1
338:       DLAMDA( K ) = D( PJ )
339:       W( K ) = Z( PJ )
340:       INDXP( K ) = PJ
341: *
342: *     Count up the total number of the various types of columns, then
343: *     form a permutation which positions the four column types into
344: *     four uniform groups (although one or more of these groups may be
345: *     empty).
346: *
347:       DO 110 J = 1, 4
348:          CTOT( J ) = 0
349:   110 CONTINUE
350:       DO 120 J = 1, N
351:          CT = COLTYP( J )
352:          CTOT( CT ) = CTOT( CT ) + 1
353:   120 CONTINUE
354: *
355: *     PSM(*) = Position in SubMatrix (of types 1 through 4)
356: *
357:       PSM( 1 ) = 1
358:       PSM( 2 ) = 1 + CTOT( 1 )
359:       PSM( 3 ) = PSM( 2 ) + CTOT( 2 )
360:       PSM( 4 ) = PSM( 3 ) + CTOT( 3 )
361:       K = N - CTOT( 4 )
362: *
363: *     Fill out the INDXC array so that the permutation which it induces
364: *     will place all type-1 columns first, all type-2 columns next,
365: *     then all type-3's, and finally all type-4's.
366: *
367:       DO 130 J = 1, N
368:          JS = INDXP( J )
369:          CT = COLTYP( JS )
370:          INDX( PSM( CT ) ) = JS
371:          INDXC( PSM( CT ) ) = J
372:          PSM( CT ) = PSM( CT ) + 1
373:   130 CONTINUE
374: *
375: *     Sort the eigenvalues and corresponding eigenvectors into DLAMDA
376: *     and Q2 respectively.  The eigenvalues/vectors which were not
377: *     deflated go into the first K slots of DLAMDA and Q2 respectively,
378: *     while those which were deflated go into the last N - K slots.
379: *
380:       I = 1
381:       IQ1 = 1
382:       IQ2 = 1 + ( CTOT( 1 )+CTOT( 2 ) )*N1
383:       DO 140 J = 1, CTOT( 1 )
384:          JS = INDX( I )
385:          CALL SCOPY( N1, Q( 1, JS ), 1, Q2( IQ1 ), 1 )
386:          Z( I ) = D( JS )
387:          I = I + 1
388:          IQ1 = IQ1 + N1
389:   140 CONTINUE
390: *
391:       DO 150 J = 1, CTOT( 2 )
392:          JS = INDX( I )
393:          CALL SCOPY( N1, Q( 1, JS ), 1, Q2( IQ1 ), 1 )
394:          CALL SCOPY( N2, Q( N1+1, JS ), 1, Q2( IQ2 ), 1 )
395:          Z( I ) = D( JS )
396:          I = I + 1
397:          IQ1 = IQ1 + N1
398:          IQ2 = IQ2 + N2
399:   150 CONTINUE
400: *
401:       DO 160 J = 1, CTOT( 3 )
402:          JS = INDX( I )
403:          CALL SCOPY( N2, Q( N1+1, JS ), 1, Q2( IQ2 ), 1 )
404:          Z( I ) = D( JS )
405:          I = I + 1
406:          IQ2 = IQ2 + N2
407:   160 CONTINUE
408: *
409:       IQ1 = IQ2
410:       DO 170 J = 1, CTOT( 4 )
411:          JS = INDX( I )
412:          CALL SCOPY( N, Q( 1, JS ), 1, Q2( IQ2 ), 1 )
413:          IQ2 = IQ2 + N
414:          Z( I ) = D( JS )
415:          I = I + 1
416:   170 CONTINUE
417: *
418: *     The deflated eigenvalues and their corresponding vectors go back
419: *     into the last N - K slots of D and Q respectively.
420: *
421:       CALL SLACPY( 'A', N, CTOT( 4 ), Q2( IQ1 ), N, Q( 1, K+1 ), LDQ )
422:       CALL SCOPY( N-K, Z( K+1 ), 1, D( K+1 ), 1 )
423: *
424: *     Copy CTOT into COLTYP for referencing in SLAED3.
425: *
426:       DO 180 J = 1, 4
427:          COLTYP( J ) = CTOT( J )
428:   180 CONTINUE
429: *
430:   190 CONTINUE
431:       RETURN
432: *
433: *     End of SLAED2
434: *
435:       END
436: