LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine dsbevd_2stage ( character  JOBZ,
character  UPLO,
integer  N,
integer  KD,
double precision, dimension( ldab, * )  AB,
integer  LDAB,
double precision, dimension( * )  W,
double precision, dimension( ldz, * )  Z,
integer  LDZ,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
integer  INFO 
)

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

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

Purpose:
 DSBEVD_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
 a real symmetric 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 DOUBLE PRECISION array, dimension (LDAB, N)
          On entry, the upper or lower triangle of the symmetric 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension 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 + N
                                   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 and IWORK
          arrays, returns these values as the first entries of the WORK
          and IWORK arrays, and no error message related to LWORK 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 JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
          If JOBZ  = 'V' and N > 2, 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 and
          IWORK arrays, returns these values as the first entries of
          the WORK and IWORK arrays, and no error message related to
          LWORK 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 236 of file dsbevd_2stage.f.

236 *
237  IMPLICIT NONE
238 *
239 * -- LAPACK driver routine (version 3.7.0) --
240 * -- LAPACK is a software package provided by Univ. of Tennessee, --
241 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
242 * December 2016
243 *
244 * .. Scalar Arguments ..
245  CHARACTER jobz, uplo
246  INTEGER info, kd, ldab, ldz, liwork, lwork, n
247 * ..
248 * .. Array Arguments ..
249  INTEGER iwork( * )
250  DOUBLE PRECISION ab( ldab, * ), w( * ), work( * ), z( ldz, * )
251 * ..
252 *
253 * =====================================================================
254 *
255 * .. Parameters ..
256  DOUBLE PRECISION zero, one
257  parameter ( zero = 0.0d+0, one = 1.0d+0 )
258 * ..
259 * .. Local Scalars ..
260  LOGICAL lower, lquery, wantz
261  INTEGER iinfo, inde, indwk2, indwrk, iscale, liwmin,
262  $ llwork, lwmin, lhtrd, lwtrd, ib, indhous,
263  $ llwrk2
264  DOUBLE PRECISION anrm, bignum, eps, rmax, rmin, safmin, sigma,
265  $ smlnum
266 * ..
267 * .. External Functions ..
268  LOGICAL lsame
269  INTEGER ilaenv
270  DOUBLE PRECISION dlamch, dlansb
271  EXTERNAL lsame, dlamch, dlansb, ilaenv
272 * ..
273 * .. External Subroutines ..
274  EXTERNAL dgemm, dlacpy, dlascl, dscal, dstedc,
275  $ dsterf, xerbla, dsytrd_sb2st
276 * ..
277 * .. Intrinsic Functions ..
278  INTRINSIC sqrt
279 * ..
280 * .. Executable Statements ..
281 *
282 * Test the input parameters.
283 *
284  wantz = lsame( jobz, 'V' )
285  lower = lsame( uplo, 'L' )
286  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
287 *
288  info = 0
289  IF( n.LE.1 ) THEN
290  liwmin = 1
291  lwmin = 1
292  ELSE
293  ib = ilaenv( 18, 'DSYTRD_SB2ST', jobz, n, kd, -1, -1 )
294  lhtrd = ilaenv( 19, 'DSYTRD_SB2ST', jobz, n, kd, ib, -1 )
295  lwtrd = ilaenv( 20, 'DSYTRD_SB2ST', jobz, n, kd, ib, -1 )
296  IF( wantz ) THEN
297  liwmin = 3 + 5*n
298  lwmin = 1 + 5*n + 2*n**2
299  ELSE
300  liwmin = 1
301  lwmin = max( 2*n, n+lhtrd+lwtrd )
302  END IF
303  END IF
304  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
305  info = -1
306  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
307  info = -2
308  ELSE IF( n.LT.0 ) THEN
309  info = -3
310  ELSE IF( kd.LT.0 ) THEN
311  info = -4
312  ELSE IF( ldab.LT.kd+1 ) THEN
313  info = -6
314  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
315  info = -9
316  END IF
317 *
318  IF( info.EQ.0 ) THEN
319  work( 1 ) = lwmin
320  iwork( 1 ) = liwmin
321 *
322  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
323  info = -11
324  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
325  info = -13
326  END IF
327  END IF
328 *
329  IF( info.NE.0 ) THEN
330  CALL xerbla( 'DSBEVD_2STAGE', -info )
331  RETURN
332  ELSE IF( lquery ) THEN
333  RETURN
334  END IF
335 *
336 * Quick return if possible
337 *
338  IF( n.EQ.0 )
339  $ RETURN
340 *
341  IF( n.EQ.1 ) THEN
342  w( 1 ) = ab( 1, 1 )
343  IF( wantz )
344  $ z( 1, 1 ) = one
345  RETURN
346  END IF
347 *
348 * Get machine constants.
349 *
350  safmin = dlamch( 'Safe minimum' )
351  eps = dlamch( 'Precision' )
352  smlnum = safmin / eps
353  bignum = one / smlnum
354  rmin = sqrt( smlnum )
355  rmax = sqrt( bignum )
356 *
357 * Scale matrix to allowable range, if necessary.
358 *
359  anrm = dlansb( 'M', uplo, n, kd, ab, ldab, work )
360  iscale = 0
361  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
362  iscale = 1
363  sigma = rmin / anrm
364  ELSE IF( anrm.GT.rmax ) THEN
365  iscale = 1
366  sigma = rmax / anrm
367  END IF
368  IF( iscale.EQ.1 ) THEN
369  IF( lower ) THEN
370  CALL dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
371  ELSE
372  CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
373  END IF
374  END IF
375 *
376 * Call DSYTRD_SB2ST to reduce band symmetric matrix to tridiagonal form.
377 *
378  inde = 1
379  indhous = inde + n
380  indwrk = indhous + lhtrd
381  llwork = lwork - indwrk + 1
382  indwk2 = indwrk + n*n
383  llwrk2 = lwork - indwk2 + 1
384 *
385  CALL dsytrd_sb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
386  $ work( inde ), work( indhous ), lhtrd,
387  $ work( indwrk ), llwork, iinfo )
388 *
389 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
390 *
391  IF( .NOT.wantz ) THEN
392  CALL dsterf( n, w, work( inde ), info )
393  ELSE
394  CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
395  $ work( indwk2 ), llwrk2, iwork, liwork, info )
396  CALL dgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
397  $ zero, work( indwk2 ), n )
398  CALL dlacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
399  END IF
400 *
401 * If matrix was scaled, then rescale eigenvalues appropriately.
402 *
403  IF( iscale.EQ.1 )
404  $ CALL dscal( n, one / sigma, w, 1 )
405 *
406  work( 1 ) = lwmin
407  iwork( 1 ) = liwmin
408  RETURN
409 *
410 * End of DSBEVD_2STAGE
411 *
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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 dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
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 dlansb(NORM, UPLO, N, K, AB, LDAB, WORK)
DLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric band matrix.
Definition: dlansb.f:131
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:191

Here is the call graph for this function:

Here is the caller graph for this function: