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