 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ cstedc()

 subroutine cstedc ( character COMPZ, integer N, real, dimension( * ) D, real, dimension( * ) E, complex, dimension( ldz, * ) Z, integer LDZ, complex, dimension( * ) WORK, integer LWORK, real, dimension( * ) RWORK, integer LRWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO )

CSTEDC

Purpose:
``` CSTEDC computes all eigenvalues and, optionally, eigenvectors of a
symmetric tridiagonal matrix using the divide and conquer method.
The eigenvectors of a full or band complex Hermitian matrix can also
be found if CHETRD or CHPTRD or CHBTRD has been used to reduce this
matrix to tridiagonal form.

This code makes very mild assumptions about floating point
arithmetic. It will work on machines with a guard digit in
add/subtract, or on those binary machines without guard digits
which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
It could conceivably fail on hexadecimal or decimal machines
without guard digits, but we know of none.  See SLAED3 for details.```
Parameters
 [in] COMPZ ``` COMPZ is CHARACTER*1 = 'N': Compute eigenvalues only. = 'I': Compute eigenvectors of tridiagonal matrix also. = 'V': Compute eigenvectors of original Hermitian matrix also. On entry, Z contains the unitary matrix used to reduce the original matrix to tridiagonal form.``` [in] N ``` N is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.``` [in,out] D ``` D is REAL array, dimension (N) On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.``` [in,out] E ``` E is REAL array, dimension (N-1) On entry, the subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.``` [in,out] Z ``` Z is COMPLEX array, dimension (LDZ,N) On entry, if COMPZ = 'V', then Z contains the unitary matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original Hermitian matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.``` [in] LDZ ``` LDZ is INTEGER The leading dimension of the array Z. LDZ >= 1. If eigenvectors are desired, then LDZ >= max(1,N).``` [out] WORK ``` WORK is COMPLEX array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.``` [in] LWORK ``` LWORK is INTEGER The dimension of the array WORK. If COMPZ = 'N' or 'I', or N <= 1, LWORK must be at least 1. If COMPZ = 'V' and N > 1, LWORK must be at least N*N. Note that for COMPZ = 'V', then if N is less than or equal to the minimum divide size, usually 25, then LWORK need only be 1. If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal sizes of the WORK, RWORK and IWORK arrays, returns these values as the first entries of the WORK, RWORK and IWORK arrays, and no error message related to LWORK or LRWORK or LIWORK is issued by XERBLA.``` [out] RWORK ``` RWORK is REAL array, dimension (MAX(1,LRWORK)) On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.``` [in] LRWORK ``` LRWORK is INTEGER The dimension of the array RWORK. If COMPZ = 'N' or N <= 1, LRWORK must be at least 1. If COMPZ = 'V' and N > 1, LRWORK must be at least 1 + 3*N + 2*N*lg N + 4*N**2 , where lg( N ) = smallest integer k such that 2**k >= N. If COMPZ = 'I' and N > 1, LRWORK must be at least 1 + 4*N + 2*N**2 . Note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LRWORK need only be max(1,2*(N-1)). If LRWORK = -1, then a workspace query is assumed; the routine only calculates the optimal sizes of the WORK, RWORK and IWORK arrays, returns these values as the first entries of the WORK, RWORK and IWORK arrays, and no error message related to LWORK or LRWORK or LIWORK is issued by XERBLA.``` [out] IWORK ``` IWORK is INTEGER array, dimension (MAX(1,LIWORK)) On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.``` [in] LIWORK ``` LIWORK is INTEGER The dimension of the array IWORK. If COMPZ = 'N' or N <= 1, LIWORK must be at least 1. If COMPZ = 'V' or N > 1, LIWORK must be at least 6 + 6*N + 5*N*lg N. If COMPZ = 'I' or N > 1, LIWORK must be at least 3 + 5*N . Note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LIWORK need only be 1. If LIWORK = -1, then a workspace query is assumed; the routine only calculates the optimal sizes of the WORK, RWORK and IWORK arrays, returns these values as the first entries of the WORK, RWORK and IWORK arrays, and no error message related to LWORK or LRWORK or LIWORK is issued by XERBLA.``` [out] INFO ``` INFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1).```
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 210 of file cstedc.f.

212*
213* -- LAPACK computational routine --
214* -- LAPACK is a software package provided by Univ. of Tennessee, --
215* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216*
217* .. Scalar Arguments ..
218 CHARACTER COMPZ
219 INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
220* ..
221* .. Array Arguments ..
222 INTEGER IWORK( * )
223 REAL D( * ), E( * ), RWORK( * )
224 COMPLEX WORK( * ), Z( LDZ, * )
225* ..
226*
227* =====================================================================
228*
229* .. Parameters ..
230 REAL ZERO, ONE, TWO
231 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
232* ..
233* .. Local Scalars ..
234 LOGICAL LQUERY
235 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL,
236 \$ LRWMIN, LWMIN, M, SMLSIZ, START
237 REAL EPS, ORGNRM, P, TINY
238* ..
239* .. External Functions ..
240 LOGICAL LSAME
241 INTEGER ILAENV
242 REAL SLAMCH, SLANST
243 EXTERNAL ilaenv, lsame, slamch, slanst
244* ..
245* .. External Subroutines ..
246 EXTERNAL xerbla, clacpy, clacrm, claed0, csteqr, cswap,
248* ..
249* .. Intrinsic Functions ..
250 INTRINSIC abs, int, log, max, mod, real, sqrt
251* ..
252* .. Executable Statements ..
253*
254* Test the input parameters.
255*
256 info = 0
257 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
258*
259 IF( lsame( compz, 'N' ) ) THEN
260 icompz = 0
261 ELSE IF( lsame( compz, 'V' ) ) THEN
262 icompz = 1
263 ELSE IF( lsame( compz, 'I' ) ) THEN
264 icompz = 2
265 ELSE
266 icompz = -1
267 END IF
268 IF( icompz.LT.0 ) THEN
269 info = -1
270 ELSE IF( n.LT.0 ) THEN
271 info = -2
272 ELSE IF( ( ldz.LT.1 ) .OR.
273 \$ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
274 info = -6
275 END IF
276*
277 IF( info.EQ.0 ) THEN
278*
279* Compute the workspace requirements
280*
281 smlsiz = ilaenv( 9, 'CSTEDC', ' ', 0, 0, 0, 0 )
282 IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
283 lwmin = 1
284 liwmin = 1
285 lrwmin = 1
286 ELSE IF( n.LE.smlsiz ) THEN
287 lwmin = 1
288 liwmin = 1
289 lrwmin = 2*( n - 1 )
290 ELSE IF( icompz.EQ.1 ) THEN
291 lgn = int( log( real( n ) ) / log( two ) )
292 IF( 2**lgn.LT.n )
293 \$ lgn = lgn + 1
294 IF( 2**lgn.LT.n )
295 \$ lgn = lgn + 1
296 lwmin = n*n
297 lrwmin = 1 + 3*n + 2*n*lgn + 4*n**2
298 liwmin = 6 + 6*n + 5*n*lgn
299 ELSE IF( icompz.EQ.2 ) THEN
300 lwmin = 1
301 lrwmin = 1 + 4*n + 2*n**2
302 liwmin = 3 + 5*n
303 END IF
304 work( 1 ) = lwmin
305 rwork( 1 ) = lrwmin
306 iwork( 1 ) = liwmin
307*
308 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
309 info = -8
310 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
311 info = -10
312 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
313 info = -12
314 END IF
315 END IF
316*
317 IF( info.NE.0 ) THEN
318 CALL xerbla( 'CSTEDC', -info )
319 RETURN
320 ELSE IF( lquery ) THEN
321 RETURN
322 END IF
323*
324* Quick return if possible
325*
326 IF( n.EQ.0 )
327 \$ RETURN
328 IF( n.EQ.1 ) THEN
329 IF( icompz.NE.0 )
330 \$ z( 1, 1 ) = one
331 RETURN
332 END IF
333*
334* If the following conditional clause is removed, then the routine
335* will use the Divide and Conquer routine to compute only the
336* eigenvalues, which requires (3N + 3N**2) real workspace and
337* (2 + 5N + 2N lg(N)) integer workspace.
338* Since on many architectures SSTERF is much faster than any other
339* algorithm for finding eigenvalues only, it is used here
340* as the default. If the conditional clause is removed, then
341* information on the size of workspace needs to be changed.
342*
343* If COMPZ = 'N', use SSTERF to compute the eigenvalues.
344*
345 IF( icompz.EQ.0 ) THEN
346 CALL ssterf( n, d, e, info )
347 GO TO 70
348 END IF
349*
350* If N is smaller than the minimum divide size (SMLSIZ+1), then
351* solve the problem with another solver.
352*
353 IF( n.LE.smlsiz ) THEN
354*
355 CALL csteqr( compz, n, d, e, z, ldz, rwork, info )
356*
357 ELSE
358*
359* If COMPZ = 'I', we simply call SSTEDC instead.
360*
361 IF( icompz.EQ.2 ) THEN
362 CALL slaset( 'Full', n, n, zero, one, rwork, n )
363 ll = n*n + 1
364 CALL sstedc( 'I', n, d, e, rwork, n,
365 \$ rwork( ll ), lrwork-ll+1, iwork, liwork, info )
366 DO 20 j = 1, n
367 DO 10 i = 1, n
368 z( i, j ) = rwork( ( j-1 )*n+i )
369 10 CONTINUE
370 20 CONTINUE
371 GO TO 70
372 END IF
373*
374* From now on, only option left to be handled is COMPZ = 'V',
375* i.e. ICOMPZ = 1.
376*
377* Scale.
378*
379 orgnrm = slanst( 'M', n, d, e )
380 IF( orgnrm.EQ.zero )
381 \$ GO TO 70
382*
383 eps = slamch( 'Epsilon' )
384*
385 start = 1
386*
387* while ( START <= N )
388*
389 30 CONTINUE
390 IF( start.LE.n ) THEN
391*
392* Let FINISH be the position of the next subdiagonal entry
393* such that E( FINISH ) <= TINY or FINISH = N if no such
394* subdiagonal exists. The matrix identified by the elements
395* between START and FINISH constitutes an independent
396* sub-problem.
397*
398 finish = start
399 40 CONTINUE
400 IF( finish.LT.n ) THEN
401 tiny = eps*sqrt( abs( d( finish ) ) )*
402 \$ sqrt( abs( d( finish+1 ) ) )
403 IF( abs( e( finish ) ).GT.tiny ) THEN
404 finish = finish + 1
405 GO TO 40
406 END IF
407 END IF
408*
409* (Sub) Problem determined. Compute its size and solve it.
410*
411 m = finish - start + 1
412 IF( m.GT.smlsiz ) THEN
413*
414* Scale.
415*
416 orgnrm = slanst( 'M', m, d( start ), e( start ) )
417 CALL slascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
418 \$ info )
419 CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
420 \$ m-1, info )
421*
422 CALL claed0( n, m, d( start ), e( start ), z( 1, start ),
423 \$ ldz, work, n, rwork, iwork, info )
424 IF( info.GT.0 ) THEN
425 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
426 \$ mod( info, ( m+1 ) ) + start - 1
427 GO TO 70
428 END IF
429*
430* Scale back.
431*
432 CALL slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
433 \$ info )
434*
435 ELSE
436 CALL ssteqr( 'I', m, d( start ), e( start ), rwork, m,
437 \$ rwork( m*m+1 ), info )
438 CALL clacrm( n, m, z( 1, start ), ldz, rwork, m, work, n,
439 \$ rwork( m*m+1 ) )
440 CALL clacpy( 'A', n, m, work, n, z( 1, start ), ldz )
441 IF( info.GT.0 ) THEN
442 info = start*( n+1 ) + finish
443 GO TO 70
444 END IF
445 END IF
446*
447 start = finish + 1
448 GO TO 30
449 END IF
450*
451* endwhile
452*
453*
454* Use Selection Sort to minimize swaps of eigenvectors
455*
456 DO 60 ii = 2, n
457 i = ii - 1
458 k = i
459 p = d( i )
460 DO 50 j = ii, n
461 IF( d( j ).LT.p ) THEN
462 k = j
463 p = d( j )
464 END IF
465 50 CONTINUE
466 IF( k.NE.i ) THEN
467 d( k ) = d( i )
468 d( i ) = p
469 CALL cswap( n, z( 1, i ), 1, z( 1, k ), 1 )
470 END IF
471 60 CONTINUE
472 END IF
473*
474 70 CONTINUE
475 work( 1 ) = lwmin
476 rwork( 1 ) = lrwmin
477 iwork( 1 ) = liwmin
478*
479 RETURN
480*
481* End of CSTEDC
482*
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
real function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: slanst.f:100
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:131
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:188
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine clacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLACRM multiplies a complex matrix by a square real matrix.
Definition: clacrm.f:114
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine claed0(QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, IWORK, INFO)
CLAED0 used by CSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: claed0.f:145
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:132
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: