LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zstedc()

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 (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
June 2017
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 214 of file zstedc.f.

214 *
215 * -- LAPACK computational routine (version 3.7.1) --
216 * -- LAPACK is a software package provided by Univ. of Tennessee, --
217 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218 * June 2017
219 *
220 * .. Scalar Arguments ..
221  CHARACTER compz
222  INTEGER info, ldz, liwork, lrwork, lwork, n
223 * ..
224 * .. Array Arguments ..
225  INTEGER iwork( * )
226  DOUBLE PRECISION d( * ), e( * ), rwork( * )
227  COMPLEX*16 work( * ), z( ldz, * )
228 * ..
229 *
230 * =====================================================================
231 *
232 * .. Parameters ..
233  DOUBLE PRECISION zero, one, two
234  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
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  DOUBLE PRECISION eps, orgnrm, p, tiny
241 * ..
242 * .. External Functions ..
243  LOGICAL lsame
244  INTEGER ilaenv
245  DOUBLE PRECISION dlamch, dlanst
246  EXTERNAL lsame, ilaenv, dlamch, dlanst
247 * ..
248 * .. External Subroutines ..
249  EXTERNAL dlascl, dlaset, dstedc, dsteqr, dsterf, xerbla,
251 * ..
252 * .. Intrinsic Functions ..
253  INTRINSIC abs, dble, int, log, max, mod, 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, 'ZSTEDC', ' ', 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( dble( 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( 'ZSTEDC', -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 DSTERF 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 DSTERF to compute the eigenvalues.
347 *
348  IF( icompz.EQ.0 ) THEN
349  CALL dsterf( 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 zsteqr( compz, n, d, e, z, ldz, rwork, info )
359 *
360  ELSE
361 *
362 * If COMPZ = 'I', we simply call DSTEDC instead.
363 *
364  IF( icompz.EQ.2 ) THEN
365  CALL dlaset( 'Full', n, n, zero, one, rwork, n )
366  ll = n*n + 1
367  CALL dstedc( '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 = dlanst( 'M', n, d, e )
383  IF( orgnrm.EQ.zero )
384  $ GO TO 70
385 *
386  eps = dlamch( '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 = dlanst( 'M', m, d( start ), e( start ) )
420  CALL dlascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
421  $ info )
422  CALL dlascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
423  $ m-1, info )
424 *
425  CALL zlaed0( 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 dlascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
436  $ info )
437 *
438  ELSE
439  CALL dsteqr( 'I', m, d( start ), e( start ), rwork, m,
440  $ rwork( m*m+1 ), info )
441  CALL zlacrm( n, m, z( 1, start ), ldz, rwork, m, work, n,
442  $ rwork( m*m+1 ) )
443  CALL zlacpy( '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 *
457 * Use Selection Sort to minimize swaps of eigenvectors
458 *
459  DO 60 ii = 2, n
460  i = ii - 1
461  k = i
462  p = d( i )
463  DO 50 j = ii, n
464  IF( d( j ).LT.p ) THEN
465  k = j
466  p = d( j )
467  END IF
468  50 CONTINUE
469  IF( k.NE.i ) THEN
470  d( k ) = d( i )
471  d( i ) = p
472  CALL zswap( n, z( 1, i ), 1, z( 1, k ), 1 )
473  END IF
474  60 CONTINUE
475  END IF
476 *
477  70 CONTINUE
478  work( 1 ) = lwmin
479  rwork( 1 ) = lrwmin
480  iwork( 1 ) = liwmin
481 *
482  RETURN
483 *
484 * End of ZSTEDC
485 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:133
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 zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
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
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:134
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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
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
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:190
Here is the call graph for this function:
Here is the caller graph for this function: