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

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.```
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).
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: