LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zstedc ( character  COMPZ,
integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
complex*16, dimension( ldz, * )  Z,
integer  LDZ,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  LRWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
integer  INFO 
)

ZSTEDC

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

Purpose:
 ZSTEDC 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 ZHETRD or ZHPTRD or ZHBTRD 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 DLAED3 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 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 DOUBLE PRECISION array,
                                         dimension (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 2015
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 215 of file zstedc.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: