LAPACK  3.10.1
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.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 210 of file zstedc.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  DOUBLE PRECISION D( * ), E( * ), RWORK( * )
224  COMPLEX*16 WORK( * ), Z( LDZ, * )
225 * ..
226 *
227 * =====================================================================
228 *
229 * .. Parameters ..
230  DOUBLE PRECISION ZERO, ONE, TWO
231  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
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  DOUBLE PRECISION EPS, ORGNRM, P, TINY
238 * ..
239 * .. External Functions ..
240  LOGICAL LSAME
241  INTEGER ILAENV
242  DOUBLE PRECISION DLAMCH, DLANST
243  EXTERNAL lsame, ilaenv, dlamch, dlanst
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL dlascl, dlaset, dstedc, dsteqr, dsterf, xerbla,
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, dble, int, log, max, mod, 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, 'ZSTEDC', ' ', 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( dble( 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( 'ZSTEDC', -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 DSTERF 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 DSTERF to compute the eigenvalues.
344 *
345  IF( icompz.EQ.0 ) THEN
346  CALL dsterf( 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 zsteqr( compz, n, d, e, z, ldz, rwork, info )
356 *
357  ELSE
358 *
359 * If COMPZ = 'I', we simply call DSTEDC instead.
360 *
361  IF( icompz.EQ.2 ) THEN
362  CALL dlaset( 'Full', n, n, zero, one, rwork, n )
363  ll = n*n + 1
364  CALL dstedc( '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 = dlanst( 'M', n, d, e )
380  IF( orgnrm.EQ.zero )
381  $ GO TO 70
382 *
383  eps = dlamch( '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 = dlanst( 'M', m, d( start ), e( start ) )
417  CALL dlascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
418  $ info )
419  CALL dlascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
420  $ m-1, info )
421 *
422  CALL zlaed0( 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 dlascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
433  $ info )
434 *
435  ELSE
436  CALL dsteqr( 'I', m, d( start ), e( start ), rwork, m,
437  $ rwork( m*m+1 ), info )
438  CALL zlacrm( n, m, z( 1, start ), ldz, rwork, m, work, n,
439  $ rwork( m*m+1 ) )
440  CALL zlacpy( '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 zswap( 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 ZSTEDC
482 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlanst.f:100
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:143
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: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 dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:188
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
ZLACRM multiplies a complex matrix by a square real matrix.
Definition: zlacrm.f:114
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:132
subroutine zlaed0(QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, IWORK, INFO)
ZLAED0 used by ZSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: zlaed0.f:145
Here is the call graph for this function:
Here is the caller graph for this function: