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