001:       SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
002:      $                   LIWORK, INFO )
003: *
004: *  -- LAPACK driver 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:       CHARACTER          COMPZ
011:       INTEGER            INFO, LDZ, LIWORK, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
022: *  symmetric tridiagonal matrix using the divide and conquer method.
023: *  The eigenvectors of a full or band real symmetric matrix can also be
024: *  found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
025: *  matrix to tridiagonal form.
026: *
027: *  This code makes very mild assumptions about floating point
028: *  arithmetic. It will work on machines with a guard digit in
029: *  add/subtract, or on those binary machines without guard digits
030: *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
031: *  It could conceivably fail on hexadecimal or decimal machines
032: *  without guard digits, but we know of none.  See DLAED3 for details.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  COMPZ   (input) CHARACTER*1
038: *          = 'N':  Compute eigenvalues only.
039: *          = 'I':  Compute eigenvectors of tridiagonal matrix also.
040: *          = 'V':  Compute eigenvectors of original dense symmetric
041: *                  matrix also.  On entry, Z contains the orthogonal
042: *                  matrix used to reduce the original matrix to
043: *                  tridiagonal form.
044: *
045: *  N       (input) INTEGER
046: *          The dimension of the symmetric tridiagonal matrix.  N >= 0.
047: *
048: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
049: *          On entry, the diagonal elements of the tridiagonal matrix.
050: *          On exit, if INFO = 0, the eigenvalues in ascending order.
051: *
052: *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
053: *          On entry, the subdiagonal elements of the tridiagonal matrix.
054: *          On exit, E has been destroyed.
055: *
056: *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
057: *          On entry, if COMPZ = 'V', then Z contains the orthogonal
058: *          matrix used in the reduction to tridiagonal form.
059: *          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
060: *          orthonormal eigenvectors of the original symmetric matrix,
061: *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
062: *          of the symmetric tridiagonal matrix.
063: *          If  COMPZ = 'N', then Z is not referenced.
064: *
065: *  LDZ     (input) INTEGER
066: *          The leading dimension of the array Z.  LDZ >= 1.
067: *          If eigenvectors are desired, then LDZ >= max(1,N).
068: *
069: *  WORK    (workspace/output) DOUBLE PRECISION array,
070: *                                         dimension (LWORK)
071: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
072: *
073: *  LWORK   (input) INTEGER
074: *          The dimension of the array WORK.
075: *          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
076: *          If COMPZ = 'V' and N > 1 then LWORK must be at least
077: *                         ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
078: *                         where lg( N ) = smallest integer k such
079: *                         that 2**k >= N.
080: *          If COMPZ = 'I' and N > 1 then LWORK must be at least
081: *                         ( 1 + 4*N + N**2 ).
082: *          Note that for COMPZ = 'I' or 'V', then if N is less than or
083: *          equal to the minimum divide size, usually 25, then LWORK need
084: *          only be max(1,2*(N-1)).
085: *
086: *          If LWORK = -1, then a workspace query is assumed; the routine
087: *          only calculates the optimal size of the WORK array, returns
088: *          this value as the first entry of the WORK array, and no error
089: *          message related to LWORK is issued by XERBLA.
090: *
091: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
092: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
093: *
094: *  LIWORK  (input) INTEGER
095: *          The dimension of the array IWORK.
096: *          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
097: *          If COMPZ = 'V' and N > 1 then LIWORK must be at least
098: *                         ( 6 + 6*N + 5*N*lg N ).
099: *          If COMPZ = 'I' and N > 1 then LIWORK must be at least
100: *                         ( 3 + 5*N ).
101: *          Note that for COMPZ = 'I' or 'V', then if N is less than or
102: *          equal to the minimum divide size, usually 25, then LIWORK
103: *          need only be 1.
104: *
105: *          If LIWORK = -1, then a workspace query is assumed; the
106: *          routine only calculates the optimal size of the IWORK array,
107: *          returns this value as the first entry of the IWORK array, and
108: *          no error message related to LIWORK is issued by XERBLA.
109: *
110: *  INFO    (output) INTEGER
111: *          = 0:  successful exit.
112: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
113: *          > 0:  The algorithm failed to compute an eigenvalue while
114: *                working on the submatrix lying in rows and columns
115: *                INFO/(N+1) through mod(INFO,N+1).
116: *
117: *  Further Details
118: *  ===============
119: *
120: *  Based on contributions by
121: *     Jeff Rutter, Computer Science Division, University of California
122: *     at Berkeley, USA
123: *  Modified by Francoise Tisseur, University of Tennessee.
124: *
125: *  =====================================================================
126: *
127: *     .. Parameters ..
128:       DOUBLE PRECISION   ZERO, ONE, TWO
129:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
130: *     ..
131: *     .. Local Scalars ..
132:       LOGICAL            LQUERY
133:       INTEGER            FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
134:      $                   LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
135:       DOUBLE PRECISION   EPS, ORGNRM, P, TINY
136: *     ..
137: *     .. External Functions ..
138:       LOGICAL            LSAME
139:       INTEGER            ILAENV
140:       DOUBLE PRECISION   DLAMCH, DLANST
141:       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANST
142: *     ..
143: *     .. External Subroutines ..
144:       EXTERNAL           DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
145:      $                   DSTEQR, DSTERF, DSWAP, XERBLA
146: *     ..
147: *     .. Intrinsic Functions ..
148:       INTRINSIC          ABS, DBLE, INT, LOG, MAX, MOD, SQRT
149: *     ..
150: *     .. Executable Statements ..
151: *
152: *     Test the input parameters.
153: *
154:       INFO = 0
155:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
156: *
157:       IF( LSAME( COMPZ, 'N' ) ) THEN
158:          ICOMPZ = 0
159:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
160:          ICOMPZ = 1
161:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
162:          ICOMPZ = 2
163:       ELSE
164:          ICOMPZ = -1
165:       END IF
166:       IF( ICOMPZ.LT.0 ) THEN
167:          INFO = -1
168:       ELSE IF( N.LT.0 ) THEN
169:          INFO = -2
170:       ELSE IF( ( LDZ.LT.1 ) .OR.
171:      $         ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
172:          INFO = -6
173:       END IF
174: *
175:       IF( INFO.EQ.0 ) THEN
176: *
177: *        Compute the workspace requirements
178: *
179:          SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
180:          IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
181:             LIWMIN = 1
182:             LWMIN = 1
183:          ELSE IF( N.LE.SMLSIZ ) THEN
184:             LIWMIN = 1
185:             LWMIN = 2*( N - 1 )
186:          ELSE
187:             LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
188:             IF( 2**LGN.LT.N )
189:      $         LGN = LGN + 1
190:             IF( 2**LGN.LT.N )
191:      $         LGN = LGN + 1
192:             IF( ICOMPZ.EQ.1 ) THEN
193:                LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2
194:                LIWMIN = 6 + 6*N + 5*N*LGN
195:             ELSE IF( ICOMPZ.EQ.2 ) THEN
196:                LWMIN = 1 + 4*N + N**2
197:                LIWMIN = 3 + 5*N
198:             END IF
199:          END IF
200:          WORK( 1 ) = LWMIN
201:          IWORK( 1 ) = LIWMIN
202: *
203:          IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
204:             INFO = -8
205:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
206:             INFO = -10
207:          END IF
208:       END IF
209: *
210:       IF( INFO.NE.0 ) THEN
211:          CALL XERBLA( 'DSTEDC', -INFO )
212:          RETURN
213:       ELSE IF (LQUERY) THEN
214:          RETURN
215:       END IF
216: *
217: *     Quick return if possible
218: *
219:       IF( N.EQ.0 )
220:      $   RETURN
221:       IF( N.EQ.1 ) THEN
222:          IF( ICOMPZ.NE.0 )
223:      $      Z( 1, 1 ) = ONE
224:          RETURN
225:       END IF
226: *
227: *     If the following conditional clause is removed, then the routine
228: *     will use the Divide and Conquer routine to compute only the
229: *     eigenvalues, which requires (3N + 3N**2) real workspace and
230: *     (2 + 5N + 2N lg(N)) integer workspace.
231: *     Since on many architectures DSTERF is much faster than any other
232: *     algorithm for finding eigenvalues only, it is used here
233: *     as the default. If the conditional clause is removed, then
234: *     information on the size of workspace needs to be changed.
235: *
236: *     If COMPZ = 'N', use DSTERF to compute the eigenvalues.
237: *
238:       IF( ICOMPZ.EQ.0 ) THEN
239:          CALL DSTERF( N, D, E, INFO )
240:          GO TO 50
241:       END IF
242: *
243: *     If N is smaller than the minimum divide size (SMLSIZ+1), then
244: *     solve the problem with another solver.
245: *
246:       IF( N.LE.SMLSIZ ) THEN
247: *
248:          CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
249: *
250:       ELSE
251: *
252: *        If COMPZ = 'V', the Z matrix must be stored elsewhere for later
253: *        use.
254: *
255:          IF( ICOMPZ.EQ.1 ) THEN
256:             STOREZ = 1 + N*N
257:          ELSE
258:             STOREZ = 1
259:          END IF
260: *
261:          IF( ICOMPZ.EQ.2 ) THEN
262:             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
263:          END IF
264: *
265: *        Scale.
266: *
267:          ORGNRM = DLANST( 'M', N, D, E )
268:          IF( ORGNRM.EQ.ZERO )
269:      $      GO TO 50
270: *
271:          EPS = DLAMCH( 'Epsilon' )
272: *
273:          START = 1
274: *
275: *        while ( START <= N )
276: *
277:    10    CONTINUE
278:          IF( START.LE.N ) THEN
279: *
280: *           Let FINISH be the position of the next subdiagonal entry
281: *           such that E( FINISH ) <= TINY or FINISH = N if no such
282: *           subdiagonal exists.  The matrix identified by the elements
283: *           between START and FINISH constitutes an independent
284: *           sub-problem.
285: *
286:             FINISH = START
287:    20       CONTINUE
288:             IF( FINISH.LT.N ) THEN
289:                TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
290:      $                    SQRT( ABS( D( FINISH+1 ) ) )
291:                IF( ABS( E( FINISH ) ).GT.TINY ) THEN
292:                   FINISH = FINISH + 1
293:                   GO TO 20
294:                END IF
295:             END IF
296: *
297: *           (Sub) Problem determined.  Compute its size and solve it.
298: *
299:             M = FINISH - START + 1
300:             IF( M.EQ.1 ) THEN
301:                START = FINISH + 1
302:                GO TO 10
303:             END IF
304:             IF( M.GT.SMLSIZ ) THEN
305: *
306: *              Scale.
307: *
308:                ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
309:                CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
310:      $                      INFO )
311:                CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
312:      $                      M-1, INFO )
313: *
314:                IF( ICOMPZ.EQ.1 ) THEN
315:                   STRTRW = 1
316:                ELSE
317:                   STRTRW = START
318:                END IF
319:                CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
320:      $                      Z( STRTRW, START ), LDZ, WORK( 1 ), N,
321:      $                      WORK( STOREZ ), IWORK, INFO )
322:                IF( INFO.NE.0 ) THEN
323:                   INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
324:      $                   MOD( INFO, ( M+1 ) ) + START - 1
325:                   GO TO 50
326:                END IF
327: *
328: *              Scale back.
329: *
330:                CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
331:      $                      INFO )
332: *
333:             ELSE
334:                IF( ICOMPZ.EQ.1 ) THEN
335: *
336: *                 Since QR won't update a Z matrix which is larger than
337: *                 the length of D, we must solve the sub-problem in a
338: *                 workspace and then multiply back into Z.
339: *
340:                   CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
341:      $                         WORK( M*M+1 ), INFO )
342:                   CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
343:      $                         WORK( STOREZ ), N )
344:                   CALL DGEMM( 'N', 'N', N, M, M, ONE,
345:      $                        WORK( STOREZ ), N, WORK, M, ZERO,
346:      $                        Z( 1, START ), LDZ )
347:                ELSE IF( ICOMPZ.EQ.2 ) THEN
348:                   CALL DSTEQR( 'I', M, D( START ), E( START ),
349:      $                         Z( START, START ), LDZ, WORK, INFO )
350:                ELSE
351:                   CALL DSTERF( M, D( START ), E( START ), INFO )
352:                END IF
353:                IF( INFO.NE.0 ) THEN
354:                   INFO = START*( N+1 ) + FINISH
355:                   GO TO 50
356:                END IF
357:             END IF
358: *
359:             START = FINISH + 1
360:             GO TO 10
361:          END IF
362: *
363: *        endwhile
364: *
365: *        If the problem split any number of times, then the eigenvalues
366: *        will not be properly ordered.  Here we permute the eigenvalues
367: *        (and the associated eigenvectors) into ascending order.
368: *
369:          IF( M.NE.N ) THEN
370:             IF( ICOMPZ.EQ.0 ) THEN
371: *
372: *              Use Quick Sort
373: *
374:                CALL DLASRT( 'I', N, D, INFO )
375: *
376:             ELSE
377: *
378: *              Use Selection Sort to minimize swaps of eigenvectors
379: *
380:                DO 40 II = 2, N
381:                   I = II - 1
382:                   K = I
383:                   P = D( I )
384:                   DO 30 J = II, N
385:                      IF( D( J ).LT.P ) THEN
386:                         K = J
387:                         P = D( J )
388:                      END IF
389:    30             CONTINUE
390:                   IF( K.NE.I ) THEN
391:                      D( K ) = D( I )
392:                      D( I ) = P
393:                      CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
394:                   END IF
395:    40          CONTINUE
396:             END IF
397:          END IF
398:       END IF
399: *
400:    50 CONTINUE
401:       WORK( 1 ) = LWMIN
402:       IWORK( 1 ) = LIWMIN
403: *
404:       RETURN
405: *
406: *     End of DSTEDC
407: *
408:       END
409: