LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine zheevd_2stage ( character  JOBZ,
character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  W,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  LRWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
integer  INFO 
)

ZHEEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
 ZHEEVD_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
 complex Hermitian matrix A using the 2stage technique for
 the reduction to tridiagonal.  If eigenvectors are desired, it uses a
 divide and conquer algorithm.

 The divide and conquer algorithm 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.
Parameters
[in]JOBZ
          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.
                  Not available in this release.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the Hermitian matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
          orthonormal eigenvectors of the matrix A.
          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
          or the upper triangle (if UPLO='U') of A, including the
          diagonal, is destroyed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]W
          W is DOUBLE PRECISION array, dimension (N)
          If INFO = 0, the eigenvalues in ascending order.
[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 N <= 1,               LWORK must be at least 1.
          If JOBZ = 'N' and N > 1, LWORK must be queried.
                                   LWORK = MAX(1, dimension) where
                                   dimension = max(stage1,stage2) + (KD+1)*N + N+1
                                             = N*KD + N*max(KD+1,FACTOPTNB) 
                                               + max(2*KD*KD, KD*NTHREADS) 
                                               + (KD+1)*N + N+1
                                   where KD is the blocking size of the reduction,
                                   FACTOPTNB is the blocking used by the QR or LQ
                                   algorithm, usually FACTOPTNB=128 is a good choice
                                   NTHREADS is the number of threads used when
                                   openMP compilation is enabled, otherwise =1.
          If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + N**2

          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 N <= 1,                LRWORK must be at least 1.
          If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
          If JOBZ  = 'V' and N > 1, LRWORK must be at least
                         1 + 5*N + 2*N**2.

          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 N <= 1,                LIWORK must be at least 1.
          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.

          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:  if INFO = i and JOBZ = 'N', then the algorithm failed
                to converge; i off-diagonal elements of an intermediate
                tridiagonal form did not converge to zero;
                if INFO = i and JOBZ = 'V', then 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
December 2016
Further Details:
Modified description of INFO. Sven, 16 Feb 05.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Further Details:
  All details about the 2stage techniques are available in:

  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
  Parallel reduction to condensed forms for symmetric eigenvalue problems
  using aggregated fine-grained and memory-aware kernels. In Proceedings
  of 2011 International Conference for High Performance Computing,
  Networking, Storage and Analysis (SC '11), New York, NY, USA,
  Article 8 , 11 pages.
  http://doi.acm.org/10.1145/2063384.2063394

  A. Haidar, J. Kurzak, P. Luszczek, 2013.
  An improved parallel singular value algorithm and its implementation 
  for multicore hardware, In Proceedings of 2013 International Conference
  for High Performance Computing, Networking, Storage and Analysis (SC '13).
  Denver, Colorado, USA, 2013.
  Article 90, 12 pages.
  http://doi.acm.org/10.1145/2503210.2503292

  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
  calculations based on fine-grained memory aware tasks.
  International Journal of High Performance Computing Applications.
  Volume 28 Issue 2, Pages 196-209, May 2014.
  http://hpc.sagepub.com/content/28/2/196 

Definition at line 255 of file zheevd_2stage.f.

255 *
256  IMPLICIT NONE
257 *
258 * -- LAPACK driver routine (version 3.7.0) --
259 * -- LAPACK is a software package provided by Univ. of Tennessee, --
260 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
261 * December 2016
262 *
263 * .. Scalar Arguments ..
264  CHARACTER jobz, uplo
265  INTEGER info, lda, liwork, lrwork, lwork, n
266 * ..
267 * .. Array Arguments ..
268  INTEGER iwork( * )
269  DOUBLE PRECISION rwork( * ), w( * )
270  COMPLEX*16 a( lda, * ), work( * )
271 * ..
272 *
273 * =====================================================================
274 *
275 * .. Parameters ..
276  DOUBLE PRECISION zero, one
277  parameter ( zero = 0.0d0, one = 1.0d0 )
278  COMPLEX*16 cone
279  parameter ( cone = ( 1.0d0, 0.0d0 ) )
280 * ..
281 * .. Local Scalars ..
282  LOGICAL lower, lquery, wantz
283  INTEGER iinfo, imax, inde, indrwk, indtau, indwk2,
284  $ indwrk, iscale, liwmin, llrwk, llwork,
285  $ llwrk2, lrwmin, lwmin,
286  $ lhtrd, lwtrd, kd, ib, indhous
287 
288 
289  DOUBLE PRECISION anrm, bignum, eps, rmax, rmin, safmin, sigma,
290  $ smlnum
291 * ..
292 * .. External Functions ..
293  LOGICAL lsame
294  INTEGER ilaenv
295  DOUBLE PRECISION dlamch, zlanhe
296  EXTERNAL lsame, ilaenv, dlamch, zlanhe
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL dscal, dsterf, xerbla, zlacpy, zlascl,
301 * ..
302 * .. Intrinsic Functions ..
303  INTRINSIC dble, max, sqrt
304 * ..
305 * .. Executable Statements ..
306 *
307 * Test the input parameters.
308 *
309  wantz = lsame( jobz, 'V' )
310  lower = lsame( uplo, 'L' )
311  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
312 *
313  info = 0
314  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
315  info = -1
316  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
317  info = -2
318  ELSE IF( n.LT.0 ) THEN
319  info = -3
320  ELSE IF( lda.LT.max( 1, n ) ) THEN
321  info = -5
322  END IF
323 *
324  IF( info.EQ.0 ) THEN
325  IF( n.LE.1 ) THEN
326  lwmin = 1
327  lrwmin = 1
328  liwmin = 1
329  ELSE
330  kd = ilaenv( 17, 'ZHETRD_2STAGE', jobz, n, -1, -1, -1 )
331  ib = ilaenv( 18, 'ZHETRD_2STAGE', jobz, n, kd, -1, -1 )
332  lhtrd = ilaenv( 19, 'ZHETRD_2STAGE', jobz, n, kd, ib, -1 )
333  lwtrd = ilaenv( 20, 'ZHETRD_2STAGE', jobz, n, kd, ib, -1 )
334  IF( wantz ) THEN
335  lwmin = 2*n + n*n
336  lrwmin = 1 + 5*n + 2*n**2
337  liwmin = 3 + 5*n
338  ELSE
339  lwmin = n + 1 + lhtrd + lwtrd
340  lrwmin = n
341  liwmin = 1
342  END IF
343  END IF
344  work( 1 ) = lwmin
345  rwork( 1 ) = lrwmin
346  iwork( 1 ) = liwmin
347 *
348  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
349  info = -8
350  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
351  info = -10
352  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
353  info = -12
354  END IF
355  END IF
356 *
357  IF( info.NE.0 ) THEN
358  CALL xerbla( 'ZHEEVD_2STAGE', -info )
359  RETURN
360  ELSE IF( lquery ) THEN
361  RETURN
362  END IF
363 *
364 * Quick return if possible
365 *
366  IF( n.EQ.0 )
367  $ RETURN
368 *
369  IF( n.EQ.1 ) THEN
370  w( 1 ) = dble( a( 1, 1 ) )
371  IF( wantz )
372  $ a( 1, 1 ) = cone
373  RETURN
374  END IF
375 *
376 * Get machine constants.
377 *
378  safmin = dlamch( 'Safe minimum' )
379  eps = dlamch( 'Precision' )
380  smlnum = safmin / eps
381  bignum = one / smlnum
382  rmin = sqrt( smlnum )
383  rmax = sqrt( bignum )
384 *
385 * Scale matrix to allowable range, if necessary.
386 *
387  anrm = zlanhe( 'M', uplo, n, a, lda, rwork )
388  iscale = 0
389  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
390  iscale = 1
391  sigma = rmin / anrm
392  ELSE IF( anrm.GT.rmax ) THEN
393  iscale = 1
394  sigma = rmax / anrm
395  END IF
396  IF( iscale.EQ.1 )
397  $ CALL zlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
398 *
399 * Call ZHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
400 *
401  inde = 1
402  indrwk = inde + n
403  llrwk = lrwork - indrwk + 1
404  indtau = 1
405  indhous = indtau + n
406  indwrk = indhous + lhtrd
407  llwork = lwork - indwrk + 1
408  indwk2 = indwrk + n*n
409  llwrk2 = lwork - indwk2 + 1
410 *
411  CALL zhetrd_2stage( jobz, uplo, n, a, lda, w, rwork( inde ),
412  $ work( indtau ), work( indhous ), lhtrd,
413  $ work( indwrk ), llwork, iinfo )
414 *
415 * For eigenvalues only, call DSTERF. For eigenvectors, first call
416 * ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
417 * tridiagonal matrix, then call ZUNMTR to multiply it to the
418 * Householder transformations represented as Householder vectors in
419 * A.
420 *
421  IF( .NOT.wantz ) THEN
422  CALL dsterf( n, w, rwork( inde ), info )
423  ELSE
424  CALL zstedc( 'I', n, w, rwork( inde ), work( indwrk ), n,
425  $ work( indwk2 ), llwrk2, rwork( indrwk ), llrwk,
426  $ iwork, liwork, info )
427  CALL zunmtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
428  $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
429  CALL zlacpy( 'A', n, n, work( indwrk ), n, a, lda )
430  END IF
431 *
432 * If matrix was scaled, then rescale eigenvalues appropriately.
433 *
434  IF( iscale.EQ.1 ) THEN
435  IF( info.EQ.0 ) THEN
436  imax = n
437  ELSE
438  imax = info - 1
439  END IF
440  CALL dscal( imax, one / sigma, w, 1 )
441  END IF
442 *
443  work( 1 ) = lwmin
444  rwork( 1 ) = lrwmin
445  iwork( 1 ) = liwmin
446 *
447  RETURN
448 *
449 * End of ZHEEVD_2STAGE
450 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
Definition: zlanhe.f:126
subroutine zhetrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
ZHETRD_2STAGE
subroutine zstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
ZSTEDC
Definition: zstedc.f:215
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine zunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMTR
Definition: zunmtr.f:173
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: