001:       SUBROUTINE SSTEDC( 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:       REAL               D( * ), E( * ), WORK( * ), Z( LDZ, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SSTEDC 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 SSYTRD or SSPTRD or SSBTRD 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 SLAED3 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) REAL 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) REAL 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) REAL 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) REAL array, dimension (MAX(1,LWORK))
070: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
071: *
072: *  LWORK   (input) INTEGER
073: *          The dimension of the array WORK.
074: *          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
075: *          If COMPZ = 'V' and N > 1 then LWORK must be at least
076: *                         ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
077: *                         where lg( N ) = smallest integer k such
078: *                         that 2**k >= N.
079: *          If COMPZ = 'I' and N > 1 then LWORK must be at least
080: *                         ( 1 + 4*N + N**2 ).
081: *          Note that for COMPZ = 'I' or 'V', then if N is less than or
082: *          equal to the minimum divide size, usually 25, then LWORK need
083: *          only be max(1,2*(N-1)).
084: *
085: *          If LWORK = -1, then a workspace query is assumed; the routine
086: *          only calculates the optimal size of the WORK array, returns
087: *          this value as the first entry of the WORK array, and no error
088: *          message related to LWORK is issued by XERBLA.
089: *
090: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
091: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
092: *
093: *  LIWORK  (input) INTEGER
094: *          The dimension of the array IWORK.
095: *          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
096: *          If COMPZ = 'V' and N > 1 then LIWORK must be at least
097: *                         ( 6 + 6*N + 5*N*lg N ).
098: *          If COMPZ = 'I' and N > 1 then LIWORK must be at least
099: *                         ( 3 + 5*N ).
100: *          Note that for COMPZ = 'I' or 'V', then if N is less than or
101: *          equal to the minimum divide size, usually 25, then LIWORK
102: *          need only be 1.
103: *
104: *          If LIWORK = -1, then a workspace query is assumed; the
105: *          routine only calculates the optimal size of the IWORK array,
106: *          returns this value as the first entry of the IWORK array, and
107: *          no error message related to LIWORK is issued by XERBLA.
108: *
109: *  INFO    (output) INTEGER
110: *          = 0:  successful exit.
111: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
112: *          > 0:  The algorithm failed to compute an eigenvalue while
113: *                working on the submatrix lying in rows and columns
114: *                INFO/(N+1) through mod(INFO,N+1).
115: *
116: *  Further Details
117: *  ===============
118: *
119: *  Based on contributions by
120: *     Jeff Rutter, Computer Science Division, University of California
121: *     at Berkeley, USA
122: *  Modified by Francoise Tisseur, University of Tennessee.
123: *
124: *  =====================================================================
125: *
126: *     .. Parameters ..
127:       REAL               ZERO, ONE, TWO
128:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
129: *     ..
130: *     .. Local Scalars ..
131:       LOGICAL            LQUERY
132:       INTEGER            FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
133:      $                   LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
134:       REAL               EPS, ORGNRM, P, TINY
135: *     ..
136: *     .. External Functions ..
137:       LOGICAL            LSAME
138:       INTEGER            ILAENV
139:       REAL               SLAMCH, SLANST
140:       EXTERNAL           ILAENV, LSAME, SLAMCH, SLANST
141: *     ..
142: *     .. External Subroutines ..
143:       EXTERNAL           SGEMM, SLACPY, SLAED0, SLASCL, SLASET, SLASRT,
144:      $                   SSTEQR, SSTERF, SSWAP, XERBLA
145: *     ..
146: *     .. Intrinsic Functions ..
147:       INTRINSIC          ABS, INT, LOG, MAX, MOD, REAL, SQRT
148: *     ..
149: *     .. Executable Statements ..
150: *
151: *     Test the input parameters.
152: *
153:       INFO = 0
154:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
155: *
156:       IF( LSAME( COMPZ, 'N' ) ) THEN
157:          ICOMPZ = 0
158:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
159:          ICOMPZ = 1
160:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
161:          ICOMPZ = 2
162:       ELSE
163:          ICOMPZ = -1
164:       END IF
165:       IF( ICOMPZ.LT.0 ) THEN
166:          INFO = -1
167:       ELSE IF( N.LT.0 ) THEN
168:          INFO = -2
169:       ELSE IF( ( LDZ.LT.1 ) .OR.
170:      $         ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
171:          INFO = -6
172:       END IF
173: *
174:       IF( INFO.EQ.0 ) THEN
175: *
176: *        Compute the workspace requirements
177: *
178:          SMLSIZ = ILAENV( 9, 'SSTEDC', ' ', 0, 0, 0, 0 )
179:          IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
180:             LIWMIN = 1
181:             LWMIN = 1
182:          ELSE IF( N.LE.SMLSIZ ) THEN
183:             LIWMIN = 1
184:             LWMIN = 2*( N - 1 )
185:          ELSE
186:             LGN = INT( LOG( REAL( N ) )/LOG( TWO ) )
187:             IF( 2**LGN.LT.N )
188:      $         LGN = LGN + 1
189:             IF( 2**LGN.LT.N )
190:      $         LGN = LGN + 1
191:             IF( ICOMPZ.EQ.1 ) THEN
192:                LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2
193:                LIWMIN = 6 + 6*N + 5*N*LGN
194:             ELSE IF( ICOMPZ.EQ.2 ) THEN
195:                LWMIN = 1 + 4*N + N**2
196:                LIWMIN = 3 + 5*N
197:             END IF
198:          END IF
199:          WORK( 1 ) = LWMIN
200:          IWORK( 1 ) = LIWMIN
201: *
202:          IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
203:             INFO = -8
204:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
205:             INFO = -10
206:          END IF
207:       END IF
208: *
209:       IF( INFO.NE.0 ) THEN
210:          CALL XERBLA( 'SSTEDC', -INFO )
211:          RETURN
212:       ELSE IF (LQUERY) THEN
213:          RETURN
214:       END IF
215: *
216: *     Quick return if possible
217: *
218:       IF( N.EQ.0 )
219:      $   RETURN
220:       IF( N.EQ.1 ) THEN
221:          IF( ICOMPZ.NE.0 )
222:      $      Z( 1, 1 ) = ONE
223:          RETURN
224:       END IF
225: *
226: *     If the following conditional clause is removed, then the routine
227: *     will use the Divide and Conquer routine to compute only the
228: *     eigenvalues, which requires (3N + 3N**2) real workspace and
229: *     (2 + 5N + 2N lg(N)) integer workspace.
230: *     Since on many architectures SSTERF is much faster than any other
231: *     algorithm for finding eigenvalues only, it is used here
232: *     as the default. If the conditional clause is removed, then
233: *     information on the size of workspace needs to be changed.
234: *
235: *     If COMPZ = 'N', use SSTERF to compute the eigenvalues.
236: *
237:       IF( ICOMPZ.EQ.0 ) THEN
238:          CALL SSTERF( N, D, E, INFO )
239:          GO TO 50
240:       END IF
241: *
242: *     If N is smaller than the minimum divide size (SMLSIZ+1), then
243: *     solve the problem with another solver.
244: *
245:       IF( N.LE.SMLSIZ ) THEN
246: *
247:          CALL SSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
248: *
249:       ELSE
250: *
251: *        If COMPZ = 'V', the Z matrix must be stored elsewhere for later
252: *        use.
253: *
254:          IF( ICOMPZ.EQ.1 ) THEN
255:             STOREZ = 1 + N*N
256:          ELSE
257:             STOREZ = 1
258:          END IF
259: *
260:          IF( ICOMPZ.EQ.2 ) THEN
261:             CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
262:          END IF
263: *
264: *        Scale.
265: *
266:          ORGNRM = SLANST( 'M', N, D, E )
267:          IF( ORGNRM.EQ.ZERO )
268:      $      GO TO 50
269: *
270:          EPS = SLAMCH( 'Epsilon' )
271: *
272:          START = 1
273: *
274: *        while ( START <= N )
275: *
276:    10    CONTINUE
277:          IF( START.LE.N ) THEN
278: *
279: *           Let FINISH be the position of the next subdiagonal entry
280: *           such that E( FINISH ) <= TINY or FINISH = N if no such
281: *           subdiagonal exists.  The matrix identified by the elements
282: *           between START and FINISH constitutes an independent
283: *           sub-problem.
284: *
285:             FINISH = START
286:    20       CONTINUE
287:             IF( FINISH.LT.N ) THEN
288:                TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
289:      $                    SQRT( ABS( D( FINISH+1 ) ) )
290:                IF( ABS( E( FINISH ) ).GT.TINY ) THEN
291:                   FINISH = FINISH + 1
292:                   GO TO 20
293:                END IF
294:             END IF
295: *
296: *           (Sub) Problem determined.  Compute its size and solve it.
297: *
298:             M = FINISH - START + 1
299:             IF( M.EQ.1 ) THEN
300:                START = FINISH + 1
301:                GO TO 10
302:             END IF
303:             IF( M.GT.SMLSIZ ) THEN
304: *
305: *              Scale.
306: *
307:                ORGNRM = SLANST( 'M', M, D( START ), E( START ) )
308:                CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
309:      $                      INFO )
310:                CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
311:      $                      M-1, INFO )
312: *
313:                IF( ICOMPZ.EQ.1 ) THEN
314:                   STRTRW = 1
315:                ELSE
316:                   STRTRW = START
317:                END IF
318:                CALL SLAED0( ICOMPZ, N, M, D( START ), E( START ),
319:      $                      Z( STRTRW, START ), LDZ, WORK( 1 ), N,
320:      $                      WORK( STOREZ ), IWORK, INFO )
321:                IF( INFO.NE.0 ) THEN
322:                   INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
323:      $                   MOD( INFO, ( M+1 ) ) + START - 1
324:                   GO TO 50
325:                END IF
326: *
327: *              Scale back.
328: *
329:                CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
330:      $                      INFO )
331: *
332:             ELSE
333:                IF( ICOMPZ.EQ.1 ) THEN
334: *
335: *                 Since QR won't update a Z matrix which is larger than
336: *                 the length of D, we must solve the sub-problem in a
337: *                 workspace and then multiply back into Z.
338: *
339:                   CALL SSTEQR( 'I', M, D( START ), E( START ), WORK, M,
340:      $                         WORK( M*M+1 ), INFO )
341:                   CALL SLACPY( 'A', N, M, Z( 1, START ), LDZ,
342:      $                         WORK( STOREZ ), N )
343:                   CALL SGEMM( 'N', 'N', N, M, M, ONE,
344:      $                        WORK( STOREZ ), N, WORK, M, ZERO,
345:      $                        Z( 1, START ), LDZ )
346:                ELSE IF( ICOMPZ.EQ.2 ) THEN
347:                   CALL SSTEQR( 'I', M, D( START ), E( START ),
348:      $                         Z( START, START ), LDZ, WORK, INFO )
349:                ELSE
350:                   CALL SSTERF( M, D( START ), E( START ), INFO )
351:                END IF
352:                IF( INFO.NE.0 ) THEN
353:                   INFO = START*( N+1 ) + FINISH
354:                   GO TO 50
355:                END IF
356:             END IF
357: *
358:             START = FINISH + 1
359:             GO TO 10
360:          END IF
361: *
362: *        endwhile
363: *
364: *        If the problem split any number of times, then the eigenvalues
365: *        will not be properly ordered.  Here we permute the eigenvalues
366: *        (and the associated eigenvectors) into ascending order.
367: *
368:          IF( M.NE.N ) THEN
369:             IF( ICOMPZ.EQ.0 ) THEN
370: *
371: *              Use Quick Sort
372: *
373:                CALL SLASRT( 'I', N, D, INFO )
374: *
375:             ELSE
376: *
377: *              Use Selection Sort to minimize swaps of eigenvectors
378: *
379:                DO 40 II = 2, N
380:                   I = II - 1
381:                   K = I
382:                   P = D( I )
383:                   DO 30 J = II, N
384:                      IF( D( J ).LT.P ) THEN
385:                         K = J
386:                         P = D( J )
387:                      END IF
388:    30             CONTINUE
389:                   IF( K.NE.I ) THEN
390:                      D( K ) = D( I )
391:                      D( I ) = P
392:                      CALL SSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
393:                   END IF
394:    40          CONTINUE
395:             END IF
396:          END IF
397:       END IF
398: *
399:    50 CONTINUE
400:       WORK( 1 ) = LWMIN
401:       IWORK( 1 ) = LIWMIN
402: *
403:       RETURN
404: *
405: *     End of SSTEDC
406: *
407:       END
408: