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