LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dsyevd_2stage()

subroutine dsyevd_2stage ( character  jobz,
character  uplo,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( * )  w,
double precision, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer  liwork,
integer  info 
)

DSYEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
 DSYEVD_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
 real symmetric matrix A using the 2stage technique for
 the reduction to tridiagonal. If eigenvectors are desired, it uses a
 divide and conquer algorithm.
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 DOUBLE PRECISION array, dimension (LDA, N)
          On entry, the symmetric 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 DOUBLE PRECISION array,
                                         dimension (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 + 2*N+1
                                             = N*KD + N*max(KD+1,FACTOPTNB)
                                               + max(2*KD*KD, KD*NTHREADS)
                                               + (KD+1)*N + 2*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
                                                1 + 6*N + 2*N**2.

          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 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 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 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.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee
Modified description of INFO. Sven, 16 Feb 05.
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 219 of file dsyevd_2stage.f.

221*
222 IMPLICIT NONE
223*
224* -- LAPACK driver routine --
225* -- LAPACK is a software package provided by Univ. of Tennessee, --
226* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227*
228* .. Scalar Arguments ..
229 CHARACTER JOBZ, UPLO
230 INTEGER INFO, LDA, LIWORK, LWORK, N
231* ..
232* .. Array Arguments ..
233 INTEGER IWORK( * )
234 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
235* ..
236*
237* =====================================================================
238*
239* .. Parameters ..
240 DOUBLE PRECISION ZERO, ONE
241 parameter( zero = 0.0d+0, one = 1.0d+0 )
242* ..
243* .. Local Scalars ..
244*
245 LOGICAL LOWER, LQUERY, WANTZ
246 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
247 $ LIWMIN, LLWORK, LLWRK2, LWMIN,
248 $ LHTRD, LWTRD, KD, IB, INDHOUS
249 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
250 $ SMLNUM
251* ..
252* .. External Functions ..
253 LOGICAL LSAME
254 INTEGER ILAENV2STAGE
255 DOUBLE PRECISION DLAMCH, DLANSY
256 EXTERNAL lsame, dlamch, dlansy, ilaenv2stage
257* ..
258* .. External Subroutines ..
259 EXTERNAL dlacpy, dlascl, dormtr, dscal, dstedc, dsterf,
261* ..
262* .. Intrinsic Functions ..
263 INTRINSIC max, sqrt
264* ..
265* .. Executable Statements ..
266*
267* Test the input parameters.
268*
269 wantz = lsame( jobz, 'V' )
270 lower = lsame( uplo, 'L' )
271 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
272*
273 info = 0
274 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
275 info = -1
276 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
277 info = -2
278 ELSE IF( n.LT.0 ) THEN
279 info = -3
280 ELSE IF( lda.LT.max( 1, n ) ) THEN
281 info = -5
282 END IF
283*
284 IF( info.EQ.0 ) THEN
285 IF( n.LE.1 ) THEN
286 liwmin = 1
287 lwmin = 1
288 ELSE
289 kd = ilaenv2stage( 1, 'DSYTRD_2STAGE', jobz,
290 $ n, -1, -1, -1 )
291 ib = ilaenv2stage( 2, 'DSYTRD_2STAGE', jobz,
292 $ n, kd, -1, -1 )
293 lhtrd = ilaenv2stage( 3, 'DSYTRD_2STAGE', jobz,
294 $ n, kd, ib, -1 )
295 lwtrd = ilaenv2stage( 4, 'DSYTRD_2STAGE', jobz,
296 $ n, kd, ib, -1 )
297 IF( wantz ) THEN
298 liwmin = 3 + 5*n
299 lwmin = 1 + 6*n + 2*n**2
300 ELSE
301 liwmin = 1
302 lwmin = 2*n + 1 + lhtrd + lwtrd
303 END IF
304 END IF
305 work( 1 ) = lwmin
306 iwork( 1 ) = liwmin
307*
308 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
309 info = -8
310 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
311 info = -10
312 END IF
313 END IF
314*
315 IF( info.NE.0 ) THEN
316 CALL xerbla( 'DSYEVD_2STAGE', -info )
317 RETURN
318 ELSE IF( lquery ) THEN
319 RETURN
320 END IF
321*
322* Quick return if possible
323*
324 IF( n.EQ.0 )
325 $ RETURN
326*
327 IF( n.EQ.1 ) THEN
328 w( 1 ) = a( 1, 1 )
329 IF( wantz )
330 $ a( 1, 1 ) = one
331 RETURN
332 END IF
333*
334* Get machine constants.
335*
336 safmin = dlamch( 'Safe minimum' )
337 eps = dlamch( 'Precision' )
338 smlnum = safmin / eps
339 bignum = one / smlnum
340 rmin = sqrt( smlnum )
341 rmax = sqrt( bignum )
342*
343* Scale matrix to allowable range, if necessary.
344*
345 anrm = dlansy( 'M', uplo, n, a, lda, work )
346 iscale = 0
347 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
348 iscale = 1
349 sigma = rmin / anrm
350 ELSE IF( anrm.GT.rmax ) THEN
351 iscale = 1
352 sigma = rmax / anrm
353 END IF
354 IF( iscale.EQ.1 )
355 $ CALL dlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
356*
357* Call DSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
358*
359 inde = 1
360 indtau = inde + n
361 indhous = indtau + n
362 indwrk = indhous + lhtrd
363 llwork = lwork - indwrk + 1
364 indwk2 = indwrk + n*n
365 llwrk2 = lwork - indwk2 + 1
366*
367 CALL dsytrd_2stage( jobz, uplo, n, a, lda, w, work( inde ),
368 $ work( indtau ), work( indhous ), lhtrd,
369 $ work( indwrk ), llwork, iinfo )
370*
371* For eigenvalues only, call DSTERF. For eigenvectors, first call
372* DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
373* tridiagonal matrix, then call DORMTR to multiply it by the
374* Householder transformations stored in A.
375*
376 IF( .NOT.wantz ) THEN
377 CALL dsterf( n, w, work( inde ), info )
378 ELSE
379* Not available in this release, and argument checking should not
380* let it getting here
381 RETURN
382 CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
383 $ work( indwk2 ), llwrk2, iwork, liwork, info )
384 CALL dormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
385 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
386 CALL dlacpy( 'A', n, n, work( indwrk ), n, a, lda )
387 END IF
388*
389* If matrix was scaled, then rescale eigenvalues appropriately.
390*
391 IF( iscale.EQ.1 )
392 $ CALL dscal( n, one / sigma, w, 1 )
393*
394 work( 1 ) = lwmin
395 iwork( 1 ) = liwmin
396*
397 RETURN
398*
399* End of DSYEVD_2STAGE
400*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsytrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
DSYTRD_2STAGE
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function dlansy(norm, uplo, n, a, lda, work)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansy.f:122
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
DSTEDC
Definition dstedc.f:182
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine dormtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
DORMTR
Definition dormtr.f:171
Here is the call graph for this function:
Here is the caller graph for this function: