LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine zhbevd_2stage ( character  JOBZ,
character  UPLO,
integer  N,
integer  KD,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
double precision, dimension( * )  W,
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 
)

ZHBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
 ZHBEVD_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
 a complex Hermitian band 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]KD
          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
[in,out]AB
          AB is COMPLEX*16 array, dimension (LDAB, N)
          On entry, the upper or lower triangle of the Hermitian band
          matrix A, stored in the first KD+1 rows of the array.  The
          j-th column of A is stored in the j-th column of the array AB
          as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

          On exit, AB is overwritten by values generated during the
          reduction to tridiagonal form.  If UPLO = 'U', the first
          superdiagonal and the diagonal of the tridiagonal matrix T
          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
          the diagonal and first subdiagonal of T are returned in the
          first two rows of AB.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD + 1.
[out]W
          W is DOUBLE PRECISION array, dimension (N)
          If INFO = 0, the eigenvalues in ascending order.
[out]Z
          Z is COMPLEX*16 array, dimension (LDZ, N)
          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
          eigenvectors of the matrix A, with the i-th column of Z
          holding the eigenvector associated with W(i).
          If JOBZ = 'N', then Z is not referenced.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', 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 length of the array WORK. LWORK >= 1, when N <= 1;
          otherwise  
          If JOBZ = 'N' and N > 1, LWORK must be queried.
                                   LWORK = MAX(1, dimension) where
                                   dimension = (2KD+1)*N + KD*NTHREADS
                                   where KD is the size of the band.
                                   NTHREADS is the number of threads used when
                                   openMP compilation is enabled, otherwise =1.
          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.

          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 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 array IWORK.
          If JOBZ = 'N' or 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, the algorithm failed to converge; i
                off-diagonal elements of an intermediate tridiagonal
                form did not converge to zero.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
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 262 of file zhbevd_2stage.f.

262 *
263  IMPLICIT NONE
264 *
265 * -- LAPACK driver routine (version 3.7.0) --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * December 2016
269 *
270 * .. Scalar Arguments ..
271  CHARACTER jobz, uplo
272  INTEGER info, kd, ldab, ldz, liwork, lrwork, lwork, n
273 * ..
274 * .. Array Arguments ..
275  INTEGER iwork( * )
276  DOUBLE PRECISION rwork( * ), w( * )
277  COMPLEX*16 ab( ldab, * ), work( * ), z( ldz, * )
278 * ..
279 *
280 * =====================================================================
281 *
282 * .. Parameters ..
283  DOUBLE PRECISION zero, one
284  parameter ( zero = 0.0d0, one = 1.0d0 )
285  COMPLEX*16 czero, cone
286  parameter ( czero = ( 0.0d0, 0.0d0 ),
287  $ cone = ( 1.0d0, 0.0d0 ) )
288 * ..
289 * .. Local Scalars ..
290  LOGICAL lower, lquery, wantz
291  INTEGER iinfo, imax, inde, indwk2, indrwk, iscale,
292  $ llwork, indwk, lhtrd, lwtrd, ib, indhous,
293  $ liwmin, llrwk, llwk2, lrwmin, lwmin
294  DOUBLE PRECISION anrm, bignum, eps, rmax, rmin, safmin, sigma,
295  $ smlnum
296 * ..
297 * .. External Functions ..
298  LOGICAL lsame
299  INTEGER ilaenv
300  DOUBLE PRECISION dlamch, zlanhb
301  EXTERNAL lsame, dlamch, zlanhb, ilaenv
302 * ..
303 * .. External Subroutines ..
304  EXTERNAL dscal, dsterf, xerbla, zgemm, zlacpy,
305  $ zlascl, zstedc, zhetrd_hb2st
306 * ..
307 * .. Intrinsic Functions ..
308  INTRINSIC dble, sqrt
309 * ..
310 * .. Executable Statements ..
311 *
312 * Test the input parameters.
313 *
314  wantz = lsame( jobz, 'V' )
315  lower = lsame( uplo, 'L' )
316  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 .OR. lrwork.EQ.-1 )
317 *
318  info = 0
319  IF( n.LE.1 ) THEN
320  lwmin = 1
321  lrwmin = 1
322  liwmin = 1
323  ELSE
324  ib = ilaenv( 18, 'ZHETRD_HB2ST', jobz, n, kd, -1, -1 )
325  lhtrd = ilaenv( 19, 'ZHETRD_HB2ST', jobz, n, kd, ib, -1 )
326  lwtrd = ilaenv( 20, 'ZHETRD_HB2ST', jobz, n, kd, ib, -1 )
327  IF( wantz ) THEN
328  lwmin = 2*n**2
329  lrwmin = 1 + 5*n + 2*n**2
330  liwmin = 3 + 5*n
331  ELSE
332  lwmin = max( n, lhtrd + lwtrd )
333  lrwmin = n
334  liwmin = 1
335  END IF
336  END IF
337  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
338  info = -1
339  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
340  info = -2
341  ELSE IF( n.LT.0 ) THEN
342  info = -3
343  ELSE IF( kd.LT.0 ) THEN
344  info = -4
345  ELSE IF( ldab.LT.kd+1 ) THEN
346  info = -6
347  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
348  info = -9
349  END IF
350 *
351  IF( info.EQ.0 ) THEN
352  work( 1 ) = lwmin
353  rwork( 1 ) = lrwmin
354  iwork( 1 ) = liwmin
355 *
356  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
357  info = -11
358  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
359  info = -13
360  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
361  info = -15
362  END IF
363  END IF
364 *
365  IF( info.NE.0 ) THEN
366  CALL xerbla( 'ZHBEVD_2STAGE', -info )
367  RETURN
368  ELSE IF( lquery ) THEN
369  RETURN
370  END IF
371 *
372 * Quick return if possible
373 *
374  IF( n.EQ.0 )
375  $ RETURN
376 *
377  IF( n.EQ.1 ) THEN
378  w( 1 ) = dble( ab( 1, 1 ) )
379  IF( wantz )
380  $ z( 1, 1 ) = cone
381  RETURN
382  END IF
383 *
384 * Get machine constants.
385 *
386  safmin = dlamch( 'Safe minimum' )
387  eps = dlamch( 'Precision' )
388  smlnum = safmin / eps
389  bignum = one / smlnum
390  rmin = sqrt( smlnum )
391  rmax = sqrt( bignum )
392 *
393 * Scale matrix to allowable range, if necessary.
394 *
395  anrm = zlanhb( 'M', uplo, n, kd, ab, ldab, rwork )
396  iscale = 0
397  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
398  iscale = 1
399  sigma = rmin / anrm
400  ELSE IF( anrm.GT.rmax ) THEN
401  iscale = 1
402  sigma = rmax / anrm
403  END IF
404  IF( iscale.EQ.1 ) THEN
405  IF( lower ) THEN
406  CALL zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
407  ELSE
408  CALL zlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
409  END IF
410  END IF
411 *
412 * Call ZHBTRD_HB2ST to reduce Hermitian band matrix to tridiagonal form.
413 *
414  inde = 1
415  indrwk = inde + n
416  llrwk = lrwork - indrwk + 1
417  indhous = 1
418  indwk = indhous + lhtrd
419  llwork = lwork - indwk + 1
420  indwk2 = indwk + n*n
421  llwk2 = lwork - indwk2 + 1
422 *
423  CALL zhetrd_hb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
424  $ rwork( inde ), work( indhous ), lhtrd,
425  $ work( indwk ), llwork, iinfo )
426 *
427 * For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEDC.
428 *
429  IF( .NOT.wantz ) THEN
430  CALL dsterf( n, w, rwork( inde ), info )
431  ELSE
432  CALL zstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
433  $ llwk2, rwork( indrwk ), llrwk, iwork, liwork,
434  $ info )
435  CALL zgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
436  $ work( indwk2 ), n )
437  CALL zlacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
438  END IF
439 *
440 * If matrix was scaled, then rescale eigenvalues appropriately.
441 *
442  IF( iscale.EQ.1 ) THEN
443  IF( info.EQ.0 ) THEN
444  imax = n
445  ELSE
446  imax = info - 1
447  END IF
448  CALL dscal( imax, one / sigma, w, 1 )
449  END IF
450 *
451  work( 1 ) = lwmin
452  rwork( 1 ) = lrwmin
453  iwork( 1 ) = liwmin
454  RETURN
455 *
456 * End of ZHBEVD_2STAGE
457 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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 zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
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 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
double precision function zlanhb(NORM, UPLO, N, K, AB, LDAB, WORK)
ZLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hermitian band matrix.
Definition: zlanhb.f:134

Here is the call graph for this function:

Here is the caller graph for this function: