LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
cstedc.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine cstedc (COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
 CSTEDC More...
 

Function/Subroutine Documentation

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

Download CSTEDC + dependencies [TGZ] [ZIP] [TXT]

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).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 214 of file cstedc.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: