LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dsbevd_2stage()

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.
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 232 of file dsbevd_2stage.f.

234 *
235  IMPLICIT NONE
236 *
237 * -- LAPACK driver routine --
238 * -- LAPACK is a software package provided by Univ. of Tennessee, --
239 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240 *
241 * .. Scalar Arguments ..
242  CHARACTER JOBZ, UPLO
243  INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
244 * ..
245 * .. Array Arguments ..
246  INTEGER IWORK( * )
247  DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
248 * ..
249 *
250 * =====================================================================
251 *
252 * .. Parameters ..
253  DOUBLE PRECISION ZERO, ONE
254  parameter( zero = 0.0d+0, one = 1.0d+0 )
255 * ..
256 * .. Local Scalars ..
257  LOGICAL LOWER, LQUERY, WANTZ
258  INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
259  $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
260  $ LLWRK2
261  DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
262  $ SMLNUM
263 * ..
264 * .. External Functions ..
265  LOGICAL LSAME
266  INTEGER ILAENV2STAGE
267  DOUBLE PRECISION DLAMCH, DLANSB
268  EXTERNAL lsame, dlamch, dlansb, ilaenv2stage
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL dgemm, dlacpy, dlascl, dscal, dstedc,
273 * ..
274 * .. Intrinsic Functions ..
275  INTRINSIC sqrt
276 * ..
277 * .. Executable Statements ..
278 *
279 * Test the input parameters.
280 *
281  wantz = lsame( jobz, 'V' )
282  lower = lsame( uplo, 'L' )
283  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
284 *
285  info = 0
286  IF( n.LE.1 ) THEN
287  liwmin = 1
288  lwmin = 1
289  ELSE
290  ib = ilaenv2stage( 2, 'DSYTRD_SB2ST', jobz, n, kd, -1, -1 )
291  lhtrd = ilaenv2stage( 3, 'DSYTRD_SB2ST', jobz, n, kd, ib, -1 )
292  lwtrd = ilaenv2stage( 4, 'DSYTRD_SB2ST', jobz, n, kd, ib, -1 )
293  IF( wantz ) THEN
294  liwmin = 3 + 5*n
295  lwmin = 1 + 5*n + 2*n**2
296  ELSE
297  liwmin = 1
298  lwmin = max( 2*n, n+lhtrd+lwtrd )
299  END IF
300  END IF
301  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
302  info = -1
303  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
304  info = -2
305  ELSE IF( n.LT.0 ) THEN
306  info = -3
307  ELSE IF( kd.LT.0 ) THEN
308  info = -4
309  ELSE IF( ldab.LT.kd+1 ) THEN
310  info = -6
311  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
312  info = -9
313  END IF
314 *
315  IF( info.EQ.0 ) THEN
316  work( 1 ) = lwmin
317  iwork( 1 ) = liwmin
318 *
319  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
320  info = -11
321  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
322  info = -13
323  END IF
324  END IF
325 *
326  IF( info.NE.0 ) THEN
327  CALL xerbla( 'DSBEVD_2STAGE', -info )
328  RETURN
329  ELSE IF( lquery ) THEN
330  RETURN
331  END IF
332 *
333 * Quick return if possible
334 *
335  IF( n.EQ.0 )
336  $ RETURN
337 *
338  IF( n.EQ.1 ) THEN
339  w( 1 ) = ab( 1, 1 )
340  IF( wantz )
341  $ z( 1, 1 ) = one
342  RETURN
343  END IF
344 *
345 * Get machine constants.
346 *
347  safmin = dlamch( 'Safe minimum' )
348  eps = dlamch( 'Precision' )
349  smlnum = safmin / eps
350  bignum = one / smlnum
351  rmin = sqrt( smlnum )
352  rmax = sqrt( bignum )
353 *
354 * Scale matrix to allowable range, if necessary.
355 *
356  anrm = dlansb( 'M', uplo, n, kd, ab, ldab, work )
357  iscale = 0
358  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
359  iscale = 1
360  sigma = rmin / anrm
361  ELSE IF( anrm.GT.rmax ) THEN
362  iscale = 1
363  sigma = rmax / anrm
364  END IF
365  IF( iscale.EQ.1 ) THEN
366  IF( lower ) THEN
367  CALL dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
368  ELSE
369  CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
370  END IF
371  END IF
372 *
373 * Call DSYTRD_SB2ST to reduce band symmetric matrix to tridiagonal form.
374 *
375  inde = 1
376  indhous = inde + n
377  indwrk = indhous + lhtrd
378  llwork = lwork - indwrk + 1
379  indwk2 = indwrk + n*n
380  llwrk2 = lwork - indwk2 + 1
381 *
382  CALL dsytrd_sb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
383  $ work( inde ), work( indhous ), lhtrd,
384  $ work( indwrk ), llwork, iinfo )
385 *
386 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
387 *
388  IF( .NOT.wantz ) THEN
389  CALL dsterf( n, w, work( inde ), info )
390  ELSE
391  CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
392  $ work( indwk2 ), llwrk2, iwork, liwork, info )
393  CALL dgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
394  $ zero, work( indwk2 ), n )
395  CALL dlacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
396  END IF
397 *
398 * If matrix was scaled, then rescale eigenvalues appropriately.
399 *
400  IF( iscale.EQ.1 )
401  $ CALL dscal( n, one / sigma, w, 1 )
402 *
403  work( 1 ) = lwmin
404  iwork( 1 ) = liwmin
405  RETURN
406 *
407 * End of DSBEVD_2STAGE
408 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dsytrd_sb2st(STAGE1, VECT, UPLO, N, KD, AB, LDAB, D, E, HOUS, LHOUS, WORK, LWORK, INFO)
DSYTRD_SB2ST reduces a real symmetric band matrix A to real symmetric tridiagonal form T
Definition: dsytrd_sb2st.F:230
integer function ilaenv2stage(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV2STAGE
Definition: ilaenv2stage.f:149
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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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 dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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,...
Definition: dlansb.f:129
Here is the call graph for this function:
Here is the caller graph for this function: