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

◆ dsyev_2stage()

subroutine dsyev_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  info 
)

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

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

Purpose:
 DSYEV_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
 real symmetric matrix A using the 2stage technique for
 the reduction to tridiagonal.
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 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 = max(stage1,stage2) + (KD+1)*N + 2*N
                                             = N*KD + N*max(KD+1,FACTOPTNB)
                                               + max(2*KD*KD, KD*NTHREADS)
                                               + (KD+1)*N + 2*N
                                   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 queried. Not yet available

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK 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 181 of file dsyev_2stage.f.

183*
184 IMPLICIT NONE
185*
186* -- LAPACK driver routine --
187* -- LAPACK is a software package provided by Univ. of Tennessee, --
188* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190* .. Scalar Arguments ..
191 CHARACTER JOBZ, UPLO
192 INTEGER INFO, LDA, LWORK, N
193* ..
194* .. Array Arguments ..
195 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
196* ..
197*
198* =====================================================================
199*
200* .. Parameters ..
201 DOUBLE PRECISION ZERO, ONE
202 parameter( zero = 0.0d0, one = 1.0d0 )
203* ..
204* .. Local Scalars ..
205 LOGICAL LOWER, LQUERY, WANTZ
206 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
207 $ LLWORK, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
208 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
209 $ SMLNUM
210* ..
211* .. External Functions ..
212 LOGICAL LSAME
213 INTEGER ILAENV2STAGE
214 DOUBLE PRECISION DLAMCH, DLANSY
215 EXTERNAL lsame, dlamch, dlansy, ilaenv2stage
216* ..
217* .. External Subroutines ..
218 EXTERNAL dlascl, dorgtr, dscal, dsteqr, dsterf,
219 $ xerbla, dsytrd_2stage
220* ..
221* .. Intrinsic Functions ..
222 INTRINSIC max, sqrt
223* ..
224* .. Executable Statements ..
225*
226* Test the input parameters.
227*
228 wantz = lsame( jobz, 'V' )
229 lower = lsame( uplo, 'L' )
230 lquery = ( lwork.EQ.-1 )
231*
232 info = 0
233 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
234 info = -1
235 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
236 info = -2
237 ELSE IF( n.LT.0 ) THEN
238 info = -3
239 ELSE IF( lda.LT.max( 1, n ) ) THEN
240 info = -5
241 END IF
242*
243 IF( info.EQ.0 ) THEN
244 kd = ilaenv2stage( 1, 'DSYTRD_2STAGE', jobz, n, -1, -1, -1 )
245 ib = ilaenv2stage( 2, 'DSYTRD_2STAGE', jobz, n, kd, -1, -1 )
246 lhtrd = ilaenv2stage( 3, 'DSYTRD_2STAGE', jobz, n, kd, ib, -1 )
247 lwtrd = ilaenv2stage( 4, 'DSYTRD_2STAGE', jobz, n, kd, ib, -1 )
248 lwmin = 2*n + lhtrd + lwtrd
249 work( 1 ) = lwmin
250*
251 IF( lwork.LT.lwmin .AND. .NOT.lquery )
252 $ info = -8
253 END IF
254*
255 IF( info.NE.0 ) THEN
256 CALL xerbla( 'DSYEV_2STAGE ', -info )
257 RETURN
258 ELSE IF( lquery ) THEN
259 RETURN
260 END IF
261*
262* Quick return if possible
263*
264 IF( n.EQ.0 ) THEN
265 RETURN
266 END IF
267*
268 IF( n.EQ.1 ) THEN
269 w( 1 ) = a( 1, 1 )
270 work( 1 ) = 2
271 IF( wantz )
272 $ a( 1, 1 ) = one
273 RETURN
274 END IF
275*
276* Get machine constants.
277*
278 safmin = dlamch( 'Safe minimum' )
279 eps = dlamch( 'Precision' )
280 smlnum = safmin / eps
281 bignum = one / smlnum
282 rmin = sqrt( smlnum )
283 rmax = sqrt( bignum )
284*
285* Scale matrix to allowable range, if necessary.
286*
287 anrm = dlansy( 'M', uplo, n, a, lda, work )
288 iscale = 0
289 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
290 iscale = 1
291 sigma = rmin / anrm
292 ELSE IF( anrm.GT.rmax ) THEN
293 iscale = 1
294 sigma = rmax / anrm
295 END IF
296 IF( iscale.EQ.1 )
297 $ CALL dlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
298*
299* Call DSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
300*
301 inde = 1
302 indtau = inde + n
303 indhous = indtau + n
304 indwrk = indhous + lhtrd
305 llwork = lwork - indwrk + 1
306*
307 CALL dsytrd_2stage( jobz, uplo, n, a, lda, w, work( inde ),
308 $ work( indtau ), work( indhous ), lhtrd,
309 $ work( indwrk ), llwork, iinfo )
310*
311* For eigenvalues only, call DSTERF. For eigenvectors, first call
312* DORGTR to generate the orthogonal matrix, then call DSTEQR.
313*
314 IF( .NOT.wantz ) THEN
315 CALL dsterf( n, w, work( inde ), info )
316 ELSE
317* Not available in this release, and argument checking should not
318* let it getting here
319 RETURN
320 CALL dorgtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
321 $ llwork, iinfo )
322 CALL dsteqr( jobz, n, w, work( inde ), a, lda, work( indtau ),
323 $ info )
324 END IF
325*
326* If matrix was scaled, then rescale eigenvalues appropriately.
327*
328 IF( iscale.EQ.1 ) THEN
329 IF( info.EQ.0 ) THEN
330 imax = n
331 ELSE
332 imax = info - 1
333 END IF
334 CALL dscal( imax, one / sigma, w, 1 )
335 END IF
336*
337* Set WORK(1) to optimal workspace size.
338*
339 work( 1 ) = lwmin
340*
341 RETURN
342*
343* End of DSYEV_2STAGE
344*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dsytrd_2stage
subroutine dsytrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
DSYTRD_2STAGE
Definition dsytrd_2stage.f:224
ilaenv2stage
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
Definition ilaenv2stage.f:149
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlansy
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
dlascl
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
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dsteqr
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:131
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
dorgtr
subroutine dorgtr(uplo, n, a, lda, tau, work, lwork, info)
DORGTR
Definition dorgtr.f:123
Here is the call graph for this function:
Here is the caller graph for this function: