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