OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
double
LAPACK » Symmetric Matrix » Eigenvalue

Functions

subroutine dsyev (jobz, uplo, n, a, lda, w, work, lwork, info)
  DSYEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine dsyev_2stage (jobz, uplo, n, a, lda, w, work, lwork, info)
  DSYEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine dsyevd (jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
  DSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine dsyevd_2stage (jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
  DSYEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine dsyevr (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, iwork, liwork, info)
  DSYEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine dsyevr_2stage (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, iwork, liwork, info)
  DSYEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine dsyevx (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
  DSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine dsyevx_2stage (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
  DSYEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine dsygv (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, info)
 DSYGV
subroutine dsygv_2stage (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, info)
 DSYGV_2STAGE
subroutine dsygvd (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, iwork, liwork, info)
 DSYGVD
subroutine dsygvx (itype, jobz, range, uplo, n, a, lda, b, ldb, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
 DSYGVX

Detailed Description

This is the group of double eigenvalue driver functions for SY matrices

Function Documentation

◆ dsyev()

subroutine dsyev ( 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 computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> DSYEV computes all eigenvalues and, optionally, eigenvectors of a
!> real symmetric matrix A.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[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 (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK.  LWORK >= max(1,3*N-1).
!>          For optimal efficiency, LWORK >= (NB+2)*N,
!>          where NB is the blocksize for DSYTRD returned by ILAENV.
!>
!>          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.

Definition at line 131 of file dsyev.f.

132*
133* -- LAPACK driver routine --
134* -- LAPACK is a software package provided by Univ. of Tennessee, --
135* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136*
137* .. Scalar Arguments ..
138 CHARACTER JOBZ, UPLO
139 INTEGER INFO, LDA, LWORK, N
140* ..
141* .. Array Arguments ..
142 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
143* ..
144*
145* =====================================================================
146*
147* .. Parameters ..
148 DOUBLE PRECISION ZERO, ONE
149 parameter( zero = 0.0d0, one = 1.0d0 )
150* ..
151* .. Local Scalars ..
152 LOGICAL LOWER, LQUERY, WANTZ
153 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
154 $ LLWORK, LWKOPT, NB
155 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
156 $ SMLNUM
157* ..
158* .. External Functions ..
159 LOGICAL LSAME
160 INTEGER ILAENV
161 DOUBLE PRECISION DLAMCH, DLANSY
162 EXTERNAL lsame, ilaenv, dlamch, dlansy
163* ..
164* .. External Subroutines ..
165 EXTERNAL dlascl, dorgtr, dscal, dsteqr, dsterf, dsytrd,
166 $ xerbla
167* ..
168* .. Intrinsic Functions ..
169 INTRINSIC max, sqrt
170* ..
171* .. Executable Statements ..
172*
173* Test the input parameters.
174*
175 wantz = lsame( jobz, 'V' )
176 lower = lsame( uplo, 'L' )
177 lquery = ( lwork.EQ.-1 )
178*
179 info = 0
180 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
181 info = -1
182 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
183 info = -2
184 ELSE IF( n.LT.0 ) THEN
185 info = -3
186 ELSE IF( lda.LT.max( 1, n ) ) THEN
187 info = -5
188 END IF
189*
190 IF( info.EQ.0 ) THEN
191 nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
192 lwkopt = max( 1, ( nb+2 )*n )
193 work( 1 ) = lwkopt
194*
195 IF( lwork.LT.max( 1, 3*n-1 ) .AND. .NOT.lquery )
196 $ info = -8
197 END IF
198*
199 IF( info.NE.0 ) THEN
200 CALL xerbla( 'DSYEV ', -info )
201 RETURN
202 ELSE IF( lquery ) THEN
203 RETURN
204 END IF
205*
206* Quick return if possible
207*
208 IF( n.EQ.0 ) THEN
209 RETURN
210 END IF
211*
212 IF( n.EQ.1 ) THEN
213 w( 1 ) = a( 1, 1 )
214 work( 1 ) = 2
215 IF( wantz )
216 $ a( 1, 1 ) = one
217 RETURN
218 END IF
219*
220* Get machine constants.
221*
222 safmin = dlamch( 'Safe minimum' )
223 eps = dlamch( 'Precision' )
224 smlnum = safmin / eps
225 bignum = one / smlnum
226 rmin = sqrt( smlnum )
227 rmax = sqrt( bignum )
228*
229* Scale matrix to allowable range, if necessary.
230*
231 anrm = dlansy( 'M', uplo, n, a, lda, work )
232 iscale = 0
233 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
234 iscale = 1
235 sigma = rmin / anrm
236 ELSE IF( anrm.GT.rmax ) THEN
237 iscale = 1
238 sigma = rmax / anrm
239 END IF
240 IF( iscale.EQ.1 )
241 $ CALL dlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
242*
243* Call DSYTRD to reduce symmetric matrix to tridiagonal form.
244*
245 inde = 1
246 indtau = inde + n
247 indwrk = indtau + n
248 llwork = lwork - indwrk + 1
249 CALL dsytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
250 $ work( indwrk ), llwork, iinfo )
251*
252* For eigenvalues only, call DSTERF. For eigenvectors, first call
253* DORGTR to generate the orthogonal matrix, then call DSTEQR.
254*
255 IF( .NOT.wantz ) THEN
256 CALL dsterf( n, w, work( inde ), info )
257 ELSE
258 CALL dorgtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
259 $ llwork, iinfo )
260 CALL dsteqr( jobz, n, w, work( inde ), a, lda, work( indtau ),
261 $ info )
262 END IF
263*
264* If matrix was scaled, then rescale eigenvalues appropriately.
265*
266 IF( iscale.EQ.1 ) THEN
267 IF( info.EQ.0 ) THEN
268 imax = n
269 ELSE
270 imax = info - 1
271 END IF
272 CALL dscal( imax, one / sigma, w, 1 )
273 END IF
274*
275* Set WORK(1) to optimal workspace size.
276*
277 work( 1 ) = lwkopt
278*
279 RETURN
280*
281* End of DSYEV
282*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
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
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
dorgtr
subroutine dorgtr(uplo, n, a, lda, tau, work, lwork, info)
DORGTR
Definition dorgtr.f:123
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
dsytrd
subroutine dsytrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
DSYTRD
Definition dsytrd.f:192
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
max
#define max(a, b)
Definition macros.h:21

◆ 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*
ilaenv2stage
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
Definition ilaenv2stage.f:149
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

◆ dsyevd()

subroutine dsyevd ( 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 computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> DSYEVD computes all eigenvalues and, optionally, eigenvectors of a
!> real symmetric matrix A. 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.
!>
!> Because of large use of BLAS of level 3, DSYEVD needs N**2 more
!> workspace than DSYEVX.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[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 at least 2*N+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.

Definition at line 183 of file dsyevd.f.

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, LIWORK, LWORK, N
193* ..
194* .. Array Arguments ..
195 INTEGER IWORK( * )
196 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
197* ..
198*
199* =====================================================================
200*
201* .. Parameters ..
202 DOUBLE PRECISION ZERO, ONE
203 parameter( zero = 0.0d+0, one = 1.0d+0 )
204* ..
205* .. Local Scalars ..
206*
207 LOGICAL LOWER, LQUERY, WANTZ
208 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
209 $ LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN
210 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
211 $ SMLNUM
212* ..
213* .. External Functions ..
214 LOGICAL LSAME
215 INTEGER ILAENV
216 DOUBLE PRECISION DLAMCH, DLANSY
217 EXTERNAL lsame, dlamch, dlansy, ilaenv
218* ..
219* .. External Subroutines ..
220 EXTERNAL dlacpy, dlascl, dormtr, dscal, dstedc, dsterf,
221 $ dsytrd, xerbla
222* ..
223* .. Intrinsic Functions ..
224 INTRINSIC max, sqrt
225* ..
226* .. Executable Statements ..
227*
228* Test the input parameters.
229*
230 wantz = lsame( jobz, 'V' )
231 lower = lsame( uplo, 'L' )
232 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
233*
234 info = 0
235 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
236 info = -1
237 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
238 info = -2
239 ELSE IF( n.LT.0 ) THEN
240 info = -3
241 ELSE IF( lda.LT.max( 1, n ) ) THEN
242 info = -5
243 END IF
244*
245 IF( info.EQ.0 ) THEN
246 IF( n.LE.1 ) THEN
247 liwmin = 1
248 lwmin = 1
249 lopt = lwmin
250 liopt = liwmin
251 ELSE
252 IF( wantz ) THEN
253 liwmin = 3 + 5*n
254 lwmin = 1 + 6*n + 2*n**2
255 ELSE
256 liwmin = 1
257 lwmin = 2*n + 1
258 END IF
259 lopt = max( lwmin, 2*n +
260 $ ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 ) )
261 liopt = liwmin
262 END IF
263 work( 1 ) = lopt
264 iwork( 1 ) = liopt
265*
266 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
267 info = -8
268 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
269 info = -10
270 END IF
271 END IF
272*
273 IF( info.NE.0 ) THEN
274 CALL xerbla( 'DSYEVD', -info )
275 RETURN
276 ELSE IF( lquery ) THEN
277 RETURN
278 END IF
279*
280* Quick return if possible
281*
282 IF( n.EQ.0 )
283 $ RETURN
284*
285 IF( n.EQ.1 ) THEN
286 w( 1 ) = a( 1, 1 )
287 IF( wantz )
288 $ a( 1, 1 ) = one
289 RETURN
290 END IF
291*
292* Get machine constants.
293*
294 safmin = dlamch( 'Safe minimum' )
295 eps = dlamch( 'Precision' )
296 smlnum = safmin / eps
297 bignum = one / smlnum
298 rmin = sqrt( smlnum )
299 rmax = sqrt( bignum )
300*
301* Scale matrix to allowable range, if necessary.
302*
303 anrm = dlansy( 'M', uplo, n, a, lda, work )
304 iscale = 0
305 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
306 iscale = 1
307 sigma = rmin / anrm
308 ELSE IF( anrm.GT.rmax ) THEN
309 iscale = 1
310 sigma = rmax / anrm
311 END IF
312 IF( iscale.EQ.1 )
313 $ CALL dlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
314*
315* Call DSYTRD to reduce symmetric matrix to tridiagonal form.
316*
317 inde = 1
318 indtau = inde + n
319 indwrk = indtau + n
320 llwork = lwork - indwrk + 1
321 indwk2 = indwrk + n*n
322 llwrk2 = lwork - indwk2 + 1
323*
324 CALL dsytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
325 $ work( indwrk ), llwork, iinfo )
326*
327* For eigenvalues only, call DSTERF. For eigenvectors, first call
328* DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
329* tridiagonal matrix, then call DORMTR to multiply it by the
330* Householder transformations stored in A.
331*
332 IF( .NOT.wantz ) THEN
333 CALL dsterf( n, w, work( inde ), info )
334 ELSE
335 CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
336 $ work( indwk2 ), llwrk2, iwork, liwork, info )
337 CALL dormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
338 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
339 CALL dlacpy( 'A', n, n, work( indwrk ), n, a, lda )
340 END IF
341*
342* If matrix was scaled, then rescale eigenvalues appropriately.
343*
344 IF( iscale.EQ.1 )
345 $ CALL dscal( n, one / sigma, w, 1 )
346*
347 work( 1 ) = lopt
348 iwork( 1 ) = liopt
349*
350 RETURN
351*
352* End of DSYEVD
353*
dlacpy
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
dstedc
subroutine dstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
DSTEDC
Definition dstedc.f:188
dormtr
subroutine dormtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
DORMTR
Definition dormtr.f:171

◆ 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.
!>
!> 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,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 225 of file dsyevd_2stage.f.

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

◆ dsyevr()

subroutine dsyevr ( character jobz,
character range,
character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
integer, dimension( * ) isuppz,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> DSYEVR computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of
!> indices for the desired eigenvalues.
!>
!> DSYEVR first reduces the matrix A to tridiagonal form T with a call
!> to DSYTRD.  Then, whenever possible, DSYEVR calls DSTEMR to compute
!> the eigenspectrum using Relatively Robust Representations.  DSTEMR
!> computes eigenvalues by the dqds algorithm, while orthogonal
!> eigenvectors are computed from various  L D L^T representations
!> (also known as Relatively Robust Representations). Gram-Schmidt
!> orthogonalization is avoided as far as possible. More specifically,
!> the various steps of the algorithm are as follows.
!>
!> For each unreduced block (submatrix) of T,
!>    (a) Compute T - sigma I  = L D L^T, so that L and D
!>        define all the wanted eigenvalues to high relative accuracy.
!>        This means that small relative changes in the entries of D and L
!>        cause only small relative changes in the eigenvalues and
!>        eigenvectors. The standard (unfactored) representation of the
!>        tridiagonal matrix T does not have this property in general.
!>    (b) Compute the eigenvalues to suitable accuracy.
!>        If the eigenvectors are desired, the algorithm attains full
!>        accuracy of the computed eigenvalues only right before
!>        the corresponding vectors have to be computed, see steps c) and d).
!>    (c) For each cluster of close eigenvalues, select a new
!>        shift close to the cluster, find a new factorization, and refine
!>        the shifted eigenvalues to suitable accuracy.
!>    (d) For each eigenvalue with a large enough relative separation compute
!>        the corresponding eigenvector by forming a rank revealing twisted
!>        factorization. Go back to (c) for any clusters that remain.
!>
!> The desired accuracy of the output can be specified by the input
!> parameter ABSTOL.
!>
!> For more details, see DSTEMR's documentation and:
!> - Inderjit S. Dhillon and Beresford N. Parlett: 
!>   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
!> - Inderjit Dhillon and Beresford Parlett:  SIAM Journal on Matrix Analysis and Applications, Vol. 25,
!>   2004.  Also LAPACK Working Note 154.
!> - Inderjit Dhillon: ,
!>   Computer Science Division Technical Report No. UCB/CSD-97-971,
!>   UC Berkeley, May 1997.
!>
!>
!> Note 1 : DSYEVR calls DSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> DSYEVR calls DSTEBZ and DSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of DSTEMR may create NaNs and infinities and
!> hence may abort due to a floating point exception in environments
!> which do not handle NaNs and infinities in the ieee standard default
!> manner.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>          For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
!>          DSTEIN are called
!>
[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, 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).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
!>          If high relative accuracy is important, set ABSTOL to
!>          DLAMCH( 'Safe minimum' ).  Doing so will guarantee that
!>          eigenvalues are computed to high relative accuracy when
!>          possible in future releases.  The current code does not
!>          make any guarantees about high relative accuracy, but
!>          future releases will. See J. Barlow and J. Demmel,
!>          , LAPACK Working Note #7, for a discussion
!>          of which matrices define their eigenvalues to high relative
!>          accuracy.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>          Supplying N columns is always safe.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]ISUPPZ
!>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
!>          The support of the eigenvectors in Z, i.e., the indices
!>          indicating the nonzero elements in Z. The i-th eigenvector
!>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
!>          ISUPPZ( 2*i ). This is an output of DSTEMR (tridiagonal
!>          matrix). The support of the eigenvectors of A is typically
!>          1:N because of the orthogonal transformations applied by DORMTR.
!>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,26*N).
!>          For optimal efficiency, LWORK >= (NB+6)*N,
!>          where NB is the max of the blocksize for DSYTRD and DORMTR
!>          returned by ILAENV.
!>
!>          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]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.  LIWORK >= max(1,10*N).
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal size of the IWORK array,
!>          returns this value as the first entry of the IWORK array, and
!>          no error message related to 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:  Internal error
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Inderjit Dhillon, IBM Almaden, USA
Osni Marques, LBNL/NERSC, USA
Ken Stanley, Computer Science Division, University of California at Berkeley, USA
Jason Riedy, Computer Science Division, University of California at Berkeley, USA

Definition at line 331 of file dsyevr.f.

334*
335* -- LAPACK driver routine --
336* -- LAPACK is a software package provided by Univ. of Tennessee, --
337* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
338*
339* .. Scalar Arguments ..
340 CHARACTER JOBZ, RANGE, UPLO
341 INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
342 DOUBLE PRECISION ABSTOL, VL, VU
343* ..
344* .. Array Arguments ..
345 INTEGER ISUPPZ( * ), IWORK( * )
346 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
347* ..
348*
349* =====================================================================
350*
351* .. Parameters ..
352 DOUBLE PRECISION ZERO, ONE, TWO
353 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
354* ..
355* .. Local Scalars ..
356 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, VALEIG, WANTZ,
357 $ TRYRAC
358 CHARACTER ORDER
359 INTEGER I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
360 $ INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
361 $ INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
362 $ LLWORK, LLWRKN, LWKOPT, LWMIN, NB, NSPLIT
363 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
364 $ SIGMA, SMLNUM, TMP1, VLL, VUU
365* ..
366* .. External Functions ..
367 LOGICAL LSAME
368 INTEGER ILAENV
369 DOUBLE PRECISION DLAMCH, DLANSY
370 EXTERNAL lsame, ilaenv, dlamch, dlansy
371* ..
372* .. External Subroutines ..
373 EXTERNAL dcopy, dormtr, dscal, dstebz, dstemr, dstein,
374 $ dsterf, dswap, dsytrd, xerbla
375* ..
376* .. Intrinsic Functions ..
377 INTRINSIC max, min, sqrt
378* ..
379* .. Executable Statements ..
380*
381* Test the input parameters.
382*
383 ieeeok = ilaenv( 10, 'DSYEVR', 'N', 1, 2, 3, 4 )
384*
385 lower = lsame( uplo, 'L' )
386 wantz = lsame( jobz, 'V' )
387 alleig = lsame( range, 'A' )
388 valeig = lsame( range, 'V' )
389 indeig = lsame( range, 'I' )
390*
391 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
392*
393 lwmin = max( 1, 26*n )
394 liwmin = max( 1, 10*n )
395*
396 info = 0
397 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
398 info = -1
399 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
400 info = -2
401 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
402 info = -3
403 ELSE IF( n.LT.0 ) THEN
404 info = -4
405 ELSE IF( lda.LT.max( 1, n ) ) THEN
406 info = -6
407 ELSE
408 IF( valeig ) THEN
409 IF( n.GT.0 .AND. vu.LE.vl )
410 $ info = -8
411 ELSE IF( indeig ) THEN
412 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
413 info = -9
414 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
415 info = -10
416 END IF
417 END IF
418 END IF
419 IF( info.EQ.0 ) THEN
420 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
421 info = -15
422 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
423 info = -18
424 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
425 info = -20
426 END IF
427 END IF
428*
429 IF( info.EQ.0 ) THEN
430 nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
431 nb = max( nb, ilaenv( 1, 'DORMTR', uplo, n, -1, -1, -1 ) )
432 lwkopt = max( ( nb+1 )*n, lwmin )
433 work( 1 ) = lwkopt
434 iwork( 1 ) = liwmin
435 END IF
436*
437 IF( info.NE.0 ) THEN
438 CALL xerbla( 'DSYEVR', -info )
439 RETURN
440 ELSE IF( lquery ) THEN
441 RETURN
442 END IF
443*
444* Quick return if possible
445*
446 m = 0
447 IF( n.EQ.0 ) THEN
448 work( 1 ) = 1
449 RETURN
450 END IF
451*
452 IF( n.EQ.1 ) THEN
453 work( 1 ) = 7
454 IF( alleig .OR. indeig ) THEN
455 m = 1
456 w( 1 ) = a( 1, 1 )
457 ELSE
458 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
459 m = 1
460 w( 1 ) = a( 1, 1 )
461 END IF
462 END IF
463 IF( wantz ) THEN
464 z( 1, 1 ) = one
465 isuppz( 1 ) = 1
466 isuppz( 2 ) = 1
467 END IF
468 RETURN
469 END IF
470*
471* Get machine constants.
472*
473 safmin = dlamch( 'Safe minimum' )
474 eps = dlamch( 'Precision' )
475 smlnum = safmin / eps
476 bignum = one / smlnum
477 rmin = sqrt( smlnum )
478 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
479*
480* Scale matrix to allowable range, if necessary.
481*
482 iscale = 0
483 abstll = abstol
484 IF (valeig) THEN
485 vll = vl
486 vuu = vu
487 END IF
488 anrm = dlansy( 'M', uplo, n, a, lda, work )
489 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
490 iscale = 1
491 sigma = rmin / anrm
492 ELSE IF( anrm.GT.rmax ) THEN
493 iscale = 1
494 sigma = rmax / anrm
495 END IF
496 IF( iscale.EQ.1 ) THEN
497 IF( lower ) THEN
498 DO 10 j = 1, n
499 CALL dscal( n-j+1, sigma, a( j, j ), 1 )
500 10 CONTINUE
501 ELSE
502 DO 20 j = 1, n
503 CALL dscal( j, sigma, a( 1, j ), 1 )
504 20 CONTINUE
505 END IF
506 IF( abstol.GT.0 )
507 $ abstll = abstol*sigma
508 IF( valeig ) THEN
509 vll = vl*sigma
510 vuu = vu*sigma
511 END IF
512 END IF
513
514* Initialize indices into workspaces. Note: The IWORK indices are
515* used only if DSTERF or DSTEMR fail.
516
517* WORK(INDTAU:INDTAU+N-1) stores the scalar factors of the
518* elementary reflectors used in DSYTRD.
519 indtau = 1
520* WORK(INDD:INDD+N-1) stores the tridiagonal's diagonal entries.
521 indd = indtau + n
522* WORK(INDE:INDE+N-1) stores the off-diagonal entries of the
523* tridiagonal matrix from DSYTRD.
524 inde = indd + n
525* WORK(INDDD:INDDD+N-1) is a copy of the diagonal entries over
526* -written by DSTEMR (the DSTERF path copies the diagonal to W).
527 inddd = inde + n
528* WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal entries over
529* -written while computing the eigenvalues in DSTERF and DSTEMR.
530 indee = inddd + n
531* INDWK is the starting offset of the left-over workspace, and
532* LLWORK is the remaining workspace size.
533 indwk = indee + n
534 llwork = lwork - indwk + 1
535
536* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
537* stores the block indices of each of the M<=N eigenvalues.
538 indibl = 1
539* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
540* stores the starting and finishing indices of each block.
541 indisp = indibl + n
542* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
543* that corresponding to eigenvectors that fail to converge in
544* DSTEIN. This information is discarded; if any fail, the driver
545* returns INFO > 0.
546 indifl = indisp + n
547* INDIWO is the offset of the remaining integer workspace.
548 indiwo = indifl + n
549
550*
551* Call DSYTRD to reduce symmetric matrix to tridiagonal form.
552*
553 CALL dsytrd( uplo, n, a, lda, work( indd ), work( inde ),
554 $ work( indtau ), work( indwk ), llwork, iinfo )
555*
556* If all eigenvalues are desired
557* then call DSTERF or DSTEMR and DORMTR.
558*
559 IF( ( alleig .OR. ( indeig .AND. il.EQ.1 .AND. iu.EQ.n ) ) .AND.
560 $ ieeeok.EQ.1 ) THEN
561 IF( .NOT.wantz ) THEN
562 CALL dcopy( n, work( indd ), 1, w, 1 )
563 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
564 CALL dsterf( n, w, work( indee ), info )
565 ELSE
566 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
567 CALL dcopy( n, work( indd ), 1, work( inddd ), 1 )
568*
569 IF (abstol .LE. two*n*eps) THEN
570 tryrac = .true.
571 ELSE
572 tryrac = .false.
573 END IF
574 CALL dstemr( jobz, 'A', n, work( inddd ), work( indee ),
575 $ vl, vu, il, iu, m, w, z, ldz, n, isuppz,
576 $ tryrac, work( indwk ), lwork, iwork, liwork,
577 $ info )
578*
579*
580*
581* Apply orthogonal matrix used in reduction to tridiagonal
582* form to eigenvectors returned by DSTEMR.
583*
584 IF( wantz .AND. info.EQ.0 ) THEN
585 indwkn = inde
586 llwrkn = lwork - indwkn + 1
587 CALL dormtr( 'L', uplo, 'N', n, m, a, lda,
588 $ work( indtau ), z, ldz, work( indwkn ),
589 $ llwrkn, iinfo )
590 END IF
591 END IF
592*
593*
594 IF( info.EQ.0 ) THEN
595* Everything worked. Skip DSTEBZ/DSTEIN. IWORK(:) are
596* undefined.
597 m = n
598 GO TO 30
599 END IF
600 info = 0
601 END IF
602*
603* Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
604* Also call DSTEBZ and DSTEIN if DSTEMR fails.
605*
606 IF( wantz ) THEN
607 order = 'B'
608 ELSE
609 order = 'E'
610 END IF
611
612 CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
613 $ work( indd ), work( inde ), m, nsplit, w,
614 $ iwork( indibl ), iwork( indisp ), work( indwk ),
615 $ iwork( indiwo ), info )
616*
617 IF( wantz ) THEN
618 CALL dstein( n, work( indd ), work( inde ), m, w,
619 $ iwork( indibl ), iwork( indisp ), z, ldz,
620 $ work( indwk ), iwork( indiwo ), iwork( indifl ),
621 $ info )
622*
623* Apply orthogonal matrix used in reduction to tridiagonal
624* form to eigenvectors returned by DSTEIN.
625*
626 indwkn = inde
627 llwrkn = lwork - indwkn + 1
628 CALL dormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
629 $ ldz, work( indwkn ), llwrkn, iinfo )
630 END IF
631*
632* If matrix was scaled, then rescale eigenvalues appropriately.
633*
634* Jump here if DSTEMR/DSTEIN succeeded.
635 30 CONTINUE
636 IF( iscale.EQ.1 ) THEN
637 IF( info.EQ.0 ) THEN
638 imax = m
639 ELSE
640 imax = info - 1
641 END IF
642 CALL dscal( imax, one / sigma, w, 1 )
643 END IF
644*
645* If eigenvalues are not in order, then sort them, along with
646* eigenvectors. Note: We do not sort the IFAIL portion of IWORK.
647* It may not be initialized (if DSTEMR/DSTEIN succeeded), and we do
648* not return this detailed information to the user.
649*
650 IF( wantz ) THEN
651 DO 50 j = 1, m - 1
652 i = 0
653 tmp1 = w( j )
654 DO 40 jj = j + 1, m
655 IF( w( jj ).LT.tmp1 ) THEN
656 i = jj
657 tmp1 = w( jj )
658 END IF
659 40 CONTINUE
660*
661 IF( i.NE.0 ) THEN
662 w( i ) = w( j )
663 w( j ) = tmp1
664 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
665 END IF
666 50 CONTINUE
667 END IF
668*
669* Set WORK(1) to optimal workspace size.
670*
671 work( 1 ) = lwkopt
672 iwork( 1 ) = liwmin
673*
674 RETURN
675*
676* End of DSYEVR
677*
dstebz
subroutine dstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
DSTEBZ
Definition dstebz.f:273
dstein
subroutine dstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
DSTEIN
Definition dstein.f:174
dstemr
subroutine dstemr(jobz, range, n, d, e, vl, vu, il, iu, m, w, z, ldz, nzc, isuppz, tryrac, work, lwork, iwork, liwork, info)
DSTEMR
Definition dstemr.f:321
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
min
#define min(a, b)
Definition macros.h:20

◆ dsyevr_2stage()

subroutine dsyevr_2stage ( character jobz,
character range,
character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
integer, dimension( * ) isuppz,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> DSYEVR_2STAGE computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A using the 2stage technique for
!> the reduction to tridiagonal.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of
!> indices for the desired eigenvalues.
!>
!> DSYEVR_2STAGE first reduces the matrix A to tridiagonal form T with a call
!> to DSYTRD.  Then, whenever possible, DSYEVR_2STAGE calls DSTEMR to compute
!> the eigenspectrum using Relatively Robust Representations.  DSTEMR
!> computes eigenvalues by the dqds algorithm, while orthogonal
!> eigenvectors are computed from various  L D L^T representations
!> (also known as Relatively Robust Representations). Gram-Schmidt
!> orthogonalization is avoided as far as possible. More specifically,
!> the various steps of the algorithm are as follows.
!>
!> For each unreduced block (submatrix) of T,
!>    (a) Compute T - sigma I  = L D L^T, so that L and D
!>        define all the wanted eigenvalues to high relative accuracy.
!>        This means that small relative changes in the entries of D and L
!>        cause only small relative changes in the eigenvalues and
!>        eigenvectors. The standard (unfactored) representation of the
!>        tridiagonal matrix T does not have this property in general.
!>    (b) Compute the eigenvalues to suitable accuracy.
!>        If the eigenvectors are desired, the algorithm attains full
!>        accuracy of the computed eigenvalues only right before
!>        the corresponding vectors have to be computed, see steps c) and d).
!>    (c) For each cluster of close eigenvalues, select a new
!>        shift close to the cluster, find a new factorization, and refine
!>        the shifted eigenvalues to suitable accuracy.
!>    (d) For each eigenvalue with a large enough relative separation compute
!>        the corresponding eigenvector by forming a rank revealing twisted
!>        factorization. Go back to (c) for any clusters that remain.
!>
!> The desired accuracy of the output can be specified by the input
!> parameter ABSTOL.
!>
!> For more details, see DSTEMR's documentation and:
!> - Inderjit S. Dhillon and Beresford N. Parlett: 
!>   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
!> - Inderjit Dhillon and Beresford Parlett:  SIAM Journal on Matrix Analysis and Applications, Vol. 25,
!>   2004.  Also LAPACK Working Note 154.
!> - Inderjit Dhillon: ,
!>   Computer Science Division Technical Report No. UCB/CSD-97-971,
!>   UC Berkeley, May 1997.
!>
!>
!> Note 1 : DSYEVR_2STAGE calls DSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> DSYEVR_2STAGE calls DSTEBZ and SSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of DSTEMR may create NaNs and infinities and
!> hence may abort due to a floating point exception in environments
!> which do not handle NaNs and infinities in the ieee standard default
!> manner.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>                  Not available in this release.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>          For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
!>          DSTEIN are called
!>
[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, 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).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
!>          If high relative accuracy is important, set ABSTOL to
!>          DLAMCH( 'Safe minimum' ).  Doing so will guarantee that
!>          eigenvalues are computed to high relative accuracy when
!>          possible in future releases.  The current code does not
!>          make any guarantees about high relative accuracy, but
!>          future releases will. See J. Barlow and J. Demmel,
!>          , LAPACK Working Note #7, for a discussion
!>          of which matrices define their eigenvalues to high relative
!>          accuracy.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>          Supplying N columns is always safe.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]ISUPPZ
!>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
!>          The support of the eigenvectors in Z, i.e., the indices
!>          indicating the nonzero elements in Z. The i-th eigenvector
!>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
!>          ISUPPZ( 2*i ). This is an output of DSTEMR (tridiagonal
!>          matrix). The support of the eigenvectors of A is typically 
!>          1:N because of the orthogonal transformations applied by DORMTR.
!>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  
!>          If JOBZ = 'N' and N > 1, LWORK must be queried.
!>                                   LWORK = MAX(1, 26*N, dimension) where
!>                                   dimension = max(stage1,stage2) + (KD+1)*N + 5*N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 5*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]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.  LIWORK >= max(1,10*N).
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal size of the IWORK array,
!>          returns this value as the first entry of the IWORK array, and
!>          no error message related to 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:  Internal error
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Inderjit Dhillon, IBM Almaden, USA
Osni Marques, LBNL/NERSC, USA
Ken Stanley, Computer Science Division, University of California at Berkeley, USA
Jason Riedy, Computer Science Division, University of California at Berkeley, USA
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 378 of file dsyevr_2stage.f.

381*
382 IMPLICIT NONE
383*
384* -- LAPACK driver routine --
385* -- LAPACK is a software package provided by Univ. of Tennessee, --
386* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
387*
388* .. Scalar Arguments ..
389 CHARACTER JOBZ, RANGE, UPLO
390 INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
391 DOUBLE PRECISION ABSTOL, VL, VU
392* ..
393* .. Array Arguments ..
394 INTEGER ISUPPZ( * ), IWORK( * )
395 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
396* ..
397*
398* =====================================================================
399*
400* .. Parameters ..
401 DOUBLE PRECISION ZERO, ONE, TWO
402 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
403* ..
404* .. Local Scalars ..
405 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, VALEIG, WANTZ,
406 $ TRYRAC
407 CHARACTER ORDER
408 INTEGER I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
409 $ INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
410 $ INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
411 $ LLWORK, LLWRKN, LWMIN, NSPLIT,
412 $ LHTRD, LWTRD, KD, IB, INDHOUS
413 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
414 $ SIGMA, SMLNUM, TMP1, VLL, VUU
415* ..
416* .. External Functions ..
417 LOGICAL LSAME
418 INTEGER ILAENV, ILAENV2STAGE
419 DOUBLE PRECISION DLAMCH, DLANSY
420 EXTERNAL lsame, dlamch, dlansy, ilaenv, ilaenv2stage
421* ..
422* .. External Subroutines ..
423 EXTERNAL dcopy, dormtr, dscal, dstebz, dstemr, dstein,
424 $ dsterf, dswap, dsytrd_2stage, xerbla
425* ..
426* .. Intrinsic Functions ..
427 INTRINSIC max, min, sqrt
428* ..
429* .. Executable Statements ..
430*
431* Test the input parameters.
432*
433 ieeeok = ilaenv( 10, 'DSYEVR', 'N', 1, 2, 3, 4 )
434*
435 lower = lsame( uplo, 'L' )
436 wantz = lsame( jobz, 'V' )
437 alleig = lsame( range, 'A' )
438 valeig = lsame( range, 'V' )
439 indeig = lsame( range, 'I' )
440*
441 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
442*
443 kd = ilaenv2stage( 1, 'DSYTRD_2STAGE', jobz, n, -1, -1, -1 )
444 ib = ilaenv2stage( 2, 'DSYTRD_2STAGE', jobz, n, kd, -1, -1 )
445 lhtrd = ilaenv2stage( 3, 'DSYTRD_2STAGE', jobz, n, kd, ib, -1 )
446 lwtrd = ilaenv2stage( 4, 'DSYTRD_2STAGE', jobz, n, kd, ib, -1 )
447 lwmin = max( 26*n, 5*n + lhtrd + lwtrd )
448 liwmin = max( 1, 10*n )
449*
450 info = 0
451 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
452 info = -1
453 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
454 info = -2
455 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
456 info = -3
457 ELSE IF( n.LT.0 ) THEN
458 info = -4
459 ELSE IF( lda.LT.max( 1, n ) ) THEN
460 info = -6
461 ELSE
462 IF( valeig ) THEN
463 IF( n.GT.0 .AND. vu.LE.vl )
464 $ info = -8
465 ELSE IF( indeig ) THEN
466 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
467 info = -9
468 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
469 info = -10
470 END IF
471 END IF
472 END IF
473 IF( info.EQ.0 ) THEN
474 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
475 info = -15
476 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
477 info = -18
478 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
479 info = -20
480 END IF
481 END IF
482*
483 IF( info.EQ.0 ) THEN
484* NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
485* NB = MAX( NB, ILAENV( 1, 'DORMTR', UPLO, N, -1, -1, -1 ) )
486* LWKOPT = MAX( ( NB+1 )*N, LWMIN )
487 work( 1 ) = lwmin
488 iwork( 1 ) = liwmin
489 END IF
490*
491 IF( info.NE.0 ) THEN
492 CALL xerbla( 'DSYEVR_2STAGE', -info )
493 RETURN
494 ELSE IF( lquery ) THEN
495 RETURN
496 END IF
497*
498* Quick return if possible
499*
500 m = 0
501 IF( n.EQ.0 ) THEN
502 work( 1 ) = 1
503 RETURN
504 END IF
505*
506 IF( n.EQ.1 ) THEN
507 work( 1 ) = 7
508 IF( alleig .OR. indeig ) THEN
509 m = 1
510 w( 1 ) = a( 1, 1 )
511 ELSE
512 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
513 m = 1
514 w( 1 ) = a( 1, 1 )
515 END IF
516 END IF
517 IF( wantz ) THEN
518 z( 1, 1 ) = one
519 isuppz( 1 ) = 1
520 isuppz( 2 ) = 1
521 END IF
522 RETURN
523 END IF
524*
525* Get machine constants.
526*
527 safmin = dlamch( 'Safe minimum' )
528 eps = dlamch( 'Precision' )
529 smlnum = safmin / eps
530 bignum = one / smlnum
531 rmin = sqrt( smlnum )
532 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
533*
534* Scale matrix to allowable range, if necessary.
535*
536 iscale = 0
537 abstll = abstol
538 IF (valeig) THEN
539 vll = vl
540 vuu = vu
541 END IF
542 anrm = dlansy( 'M', uplo, n, a, lda, work )
543 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
544 iscale = 1
545 sigma = rmin / anrm
546 ELSE IF( anrm.GT.rmax ) THEN
547 iscale = 1
548 sigma = rmax / anrm
549 END IF
550 IF( iscale.EQ.1 ) THEN
551 IF( lower ) THEN
552 DO 10 j = 1, n
553 CALL dscal( n-j+1, sigma, a( j, j ), 1 )
554 10 CONTINUE
555 ELSE
556 DO 20 j = 1, n
557 CALL dscal( j, sigma, a( 1, j ), 1 )
558 20 CONTINUE
559 END IF
560 IF( abstol.GT.0 )
561 $ abstll = abstol*sigma
562 IF( valeig ) THEN
563 vll = vl*sigma
564 vuu = vu*sigma
565 END IF
566 END IF
567
568* Initialize indices into workspaces. Note: The IWORK indices are
569* used only if DSTERF or DSTEMR fail.
570
571* WORK(INDTAU:INDTAU+N-1) stores the scalar factors of the
572* elementary reflectors used in DSYTRD.
573 indtau = 1
574* WORK(INDD:INDD+N-1) stores the tridiagonal's diagonal entries.
575 indd = indtau + n
576* WORK(INDE:INDE+N-1) stores the off-diagonal entries of the
577* tridiagonal matrix from DSYTRD.
578 inde = indd + n
579* WORK(INDDD:INDDD+N-1) is a copy of the diagonal entries over
580* -written by DSTEMR (the DSTERF path copies the diagonal to W).
581 inddd = inde + n
582* WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal entries over
583* -written while computing the eigenvalues in DSTERF and DSTEMR.
584 indee = inddd + n
585* INDHOUS is the starting offset Householder storage of stage 2
586 indhous = indee + n
587* INDWK is the starting offset of the left-over workspace, and
588* LLWORK is the remaining workspace size.
589 indwk = indhous + lhtrd
590 llwork = lwork - indwk + 1
591
592
593* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
594* stores the block indices of each of the M<=N eigenvalues.
595 indibl = 1
596* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
597* stores the starting and finishing indices of each block.
598 indisp = indibl + n
599* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
600* that corresponding to eigenvectors that fail to converge in
601* DSTEIN. This information is discarded; if any fail, the driver
602* returns INFO > 0.
603 indifl = indisp + n
604* INDIWO is the offset of the remaining integer workspace.
605 indiwo = indifl + n
606
607*
608* Call DSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
609*
610*
611 CALL dsytrd_2stage( jobz, uplo, n, a, lda, work( indd ),
612 $ work( inde ), work( indtau ), work( indhous ),
613 $ lhtrd, work( indwk ), llwork, iinfo )
614*
615* If all eigenvalues are desired
616* then call DSTERF or DSTEMR and DORMTR.
617*
618 IF( ( alleig .OR. ( indeig .AND. il.EQ.1 .AND. iu.EQ.n ) ) .AND.
619 $ ieeeok.EQ.1 ) THEN
620 IF( .NOT.wantz ) THEN
621 CALL dcopy( n, work( indd ), 1, w, 1 )
622 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
623 CALL dsterf( n, w, work( indee ), info )
624 ELSE
625 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
626 CALL dcopy( n, work( indd ), 1, work( inddd ), 1 )
627*
628 IF (abstol .LE. two*n*eps) THEN
629 tryrac = .true.
630 ELSE
631 tryrac = .false.
632 END IF
633 CALL dstemr( jobz, 'A', n, work( inddd ), work( indee ),
634 $ vl, vu, il, iu, m, w, z, ldz, n, isuppz,
635 $ tryrac, work( indwk ), lwork, iwork, liwork,
636 $ info )
637*
638*
639*
640* Apply orthogonal matrix used in reduction to tridiagonal
641* form to eigenvectors returned by DSTEMR.
642*
643 IF( wantz .AND. info.EQ.0 ) THEN
644 indwkn = inde
645 llwrkn = lwork - indwkn + 1
646 CALL dormtr( 'L', uplo, 'N', n, m, a, lda,
647 $ work( indtau ), z, ldz, work( indwkn ),
648 $ llwrkn, iinfo )
649 END IF
650 END IF
651*
652*
653 IF( info.EQ.0 ) THEN
654* Everything worked. Skip DSTEBZ/DSTEIN. IWORK(:) are
655* undefined.
656 m = n
657 GO TO 30
658 END IF
659 info = 0
660 END IF
661*
662* Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
663* Also call DSTEBZ and DSTEIN if DSTEMR fails.
664*
665 IF( wantz ) THEN
666 order = 'B'
667 ELSE
668 order = 'E'
669 END IF
670
671 CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
672 $ work( indd ), work( inde ), m, nsplit, w,
673 $ iwork( indibl ), iwork( indisp ), work( indwk ),
674 $ iwork( indiwo ), info )
675*
676 IF( wantz ) THEN
677 CALL dstein( n, work( indd ), work( inde ), m, w,
678 $ iwork( indibl ), iwork( indisp ), z, ldz,
679 $ work( indwk ), iwork( indiwo ), iwork( indifl ),
680 $ info )
681*
682* Apply orthogonal matrix used in reduction to tridiagonal
683* form to eigenvectors returned by DSTEIN.
684*
685 indwkn = inde
686 llwrkn = lwork - indwkn + 1
687 CALL dormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
688 $ ldz, work( indwkn ), llwrkn, iinfo )
689 END IF
690*
691* If matrix was scaled, then rescale eigenvalues appropriately.
692*
693* Jump here if DSTEMR/DSTEIN succeeded.
694 30 CONTINUE
695 IF( iscale.EQ.1 ) THEN
696 IF( info.EQ.0 ) THEN
697 imax = m
698 ELSE
699 imax = info - 1
700 END IF
701 CALL dscal( imax, one / sigma, w, 1 )
702 END IF
703*
704* If eigenvalues are not in order, then sort them, along with
705* eigenvectors. Note: We do not sort the IFAIL portion of IWORK.
706* It may not be initialized (if DSTEMR/DSTEIN succeeded), and we do
707* not return this detailed information to the user.
708*
709 IF( wantz ) THEN
710 DO 50 j = 1, m - 1
711 i = 0
712 tmp1 = w( j )
713 DO 40 jj = j + 1, m
714 IF( w( jj ).LT.tmp1 ) THEN
715 i = jj
716 tmp1 = w( jj )
717 END IF
718 40 CONTINUE
719*
720 IF( i.NE.0 ) THEN
721 w( i ) = w( j )
722 w( j ) = tmp1
723 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
724 END IF
725 50 CONTINUE
726 END IF
727*
728* Set WORK(1) to optimal workspace size.
729*
730 work( 1 ) = lwmin
731 iwork( 1 ) = liwmin
732*
733 RETURN
734*
735* End of DSYEVR_2STAGE
736*

◆ dsyevx()

subroutine dsyevx ( character jobz,
character range,
character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

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

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

Purpose:
!>
!> DSYEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of indices
!> for the desired eigenvalues.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[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, 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).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[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 (MAX(1,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 8*N.
!>          For optimal efficiency, LWORK >= (NB+3)*N,
!>          where NB is the max of the blocksize for DSYTRD and DORMTR
!>          returned by ILAENV.
!>
!>          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]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, then i eigenvectors failed to converge.
!>                Their indices are stored in array IFAIL.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 250 of file dsyevx.f.

253*
254* -- LAPACK driver routine --
255* -- LAPACK is a software package provided by Univ. of Tennessee, --
256* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
257*
258* .. Scalar Arguments ..
259 CHARACTER JOBZ, RANGE, UPLO
260 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
261 DOUBLE PRECISION ABSTOL, VL, VU
262* ..
263* .. Array Arguments ..
264 INTEGER IFAIL( * ), IWORK( * )
265 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
266* ..
267*
268* =====================================================================
269*
270* .. Parameters ..
271 DOUBLE PRECISION ZERO, ONE
272 parameter( zero = 0.0d+0, one = 1.0d+0 )
273* ..
274* .. Local Scalars ..
275 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
276 $ WANTZ
277 CHARACTER ORDER
278 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
279 $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
280 $ ITMP1, J, JJ, LLWORK, LLWRKN, LWKMIN,
281 $ LWKOPT, NB, NSPLIT
282 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
283 $ SIGMA, SMLNUM, TMP1, VLL, VUU
284* ..
285* .. External Functions ..
286 LOGICAL LSAME
287 INTEGER ILAENV
288 DOUBLE PRECISION DLAMCH, DLANSY
289 EXTERNAL lsame, ilaenv, dlamch, dlansy
290* ..
291* .. External Subroutines ..
292 EXTERNAL dcopy, dlacpy, dorgtr, dormtr, dscal, dstebz,
293 $ dstein, dsteqr, dsterf, dswap, dsytrd, xerbla
294* ..
295* .. Intrinsic Functions ..
296 INTRINSIC max, min, sqrt
297* ..
298* .. Executable Statements ..
299*
300* Test the input parameters.
301*
302 lower = lsame( uplo, 'L' )
303 wantz = lsame( jobz, 'V' )
304 alleig = lsame( range, 'A' )
305 valeig = lsame( range, 'V' )
306 indeig = lsame( range, 'I' )
307 lquery = ( lwork.EQ.-1 )
308*
309 info = 0
310 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
311 info = -1
312 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
313 info = -2
314 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
315 info = -3
316 ELSE IF( n.LT.0 ) THEN
317 info = -4
318 ELSE IF( lda.LT.max( 1, n ) ) THEN
319 info = -6
320 ELSE
321 IF( valeig ) THEN
322 IF( n.GT.0 .AND. vu.LE.vl )
323 $ info = -8
324 ELSE IF( indeig ) THEN
325 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
326 info = -9
327 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
328 info = -10
329 END IF
330 END IF
331 END IF
332 IF( info.EQ.0 ) THEN
333 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
334 info = -15
335 END IF
336 END IF
337*
338 IF( info.EQ.0 ) THEN
339 IF( n.LE.1 ) THEN
340 lwkmin = 1
341 work( 1 ) = lwkmin
342 ELSE
343 lwkmin = 8*n
344 nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
345 nb = max( nb, ilaenv( 1, 'DORMTR', uplo, n, -1, -1, -1 ) )
346 lwkopt = max( lwkmin, ( nb + 3 )*n )
347 work( 1 ) = lwkopt
348 END IF
349*
350 IF( lwork.LT.lwkmin .AND. .NOT.lquery )
351 $ info = -17
352 END IF
353*
354 IF( info.NE.0 ) THEN
355 CALL xerbla( 'DSYEVX', -info )
356 RETURN
357 ELSE IF( lquery ) THEN
358 RETURN
359 END IF
360*
361* Quick return if possible
362*
363 m = 0
364 IF( n.EQ.0 ) THEN
365 RETURN
366 END IF
367*
368 IF( n.EQ.1 ) THEN
369 IF( alleig .OR. indeig ) THEN
370 m = 1
371 w( 1 ) = a( 1, 1 )
372 ELSE
373 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
374 m = 1
375 w( 1 ) = a( 1, 1 )
376 END IF
377 END IF
378 IF( wantz )
379 $ z( 1, 1 ) = one
380 RETURN
381 END IF
382*
383* Get machine constants.
384*
385 safmin = dlamch( 'Safe minimum' )
386 eps = dlamch( 'Precision' )
387 smlnum = safmin / eps
388 bignum = one / smlnum
389 rmin = sqrt( smlnum )
390 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
391*
392* Scale matrix to allowable range, if necessary.
393*
394 iscale = 0
395 abstll = abstol
396 IF( valeig ) THEN
397 vll = vl
398 vuu = vu
399 END IF
400 anrm = dlansy( 'M', uplo, n, a, lda, work )
401 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
402 iscale = 1
403 sigma = rmin / anrm
404 ELSE IF( anrm.GT.rmax ) THEN
405 iscale = 1
406 sigma = rmax / anrm
407 END IF
408 IF( iscale.EQ.1 ) THEN
409 IF( lower ) THEN
410 DO 10 j = 1, n
411 CALL dscal( n-j+1, sigma, a( j, j ), 1 )
412 10 CONTINUE
413 ELSE
414 DO 20 j = 1, n
415 CALL dscal( j, sigma, a( 1, j ), 1 )
416 20 CONTINUE
417 END IF
418 IF( abstol.GT.0 )
419 $ abstll = abstol*sigma
420 IF( valeig ) THEN
421 vll = vl*sigma
422 vuu = vu*sigma
423 END IF
424 END IF
425*
426* Call DSYTRD to reduce symmetric matrix to tridiagonal form.
427*
428 indtau = 1
429 inde = indtau + n
430 indd = inde + n
431 indwrk = indd + n
432 llwork = lwork - indwrk + 1
433 CALL dsytrd( uplo, n, a, lda, work( indd ), work( inde ),
434 $ work( indtau ), work( indwrk ), llwork, iinfo )
435*
436* If all eigenvalues are desired and ABSTOL is less than or equal to
437* zero, then call DSTERF or DORGTR and SSTEQR. If this fails for
438* some eigenvalue, then try DSTEBZ.
439*
440 test = .false.
441 IF( indeig ) THEN
442 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
443 test = .true.
444 END IF
445 END IF
446 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
447 CALL dcopy( n, work( indd ), 1, w, 1 )
448 indee = indwrk + 2*n
449 IF( .NOT.wantz ) THEN
450 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
451 CALL dsterf( n, w, work( indee ), info )
452 ELSE
453 CALL dlacpy( 'A', n, n, a, lda, z, ldz )
454 CALL dorgtr( uplo, n, z, ldz, work( indtau ),
455 $ work( indwrk ), llwork, iinfo )
456 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
457 CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
458 $ work( indwrk ), info )
459 IF( info.EQ.0 ) THEN
460 DO 30 i = 1, n
461 ifail( i ) = 0
462 30 CONTINUE
463 END IF
464 END IF
465 IF( info.EQ.0 ) THEN
466 m = n
467 GO TO 40
468 END IF
469 info = 0
470 END IF
471*
472* Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
473*
474 IF( wantz ) THEN
475 order = 'B'
476 ELSE
477 order = 'E'
478 END IF
479 indibl = 1
480 indisp = indibl + n
481 indiwo = indisp + n
482 CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
483 $ work( indd ), work( inde ), m, nsplit, w,
484 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
485 $ iwork( indiwo ), info )
486*
487 IF( wantz ) THEN
488 CALL dstein( n, work( indd ), work( inde ), m, w,
489 $ iwork( indibl ), iwork( indisp ), z, ldz,
490 $ work( indwrk ), iwork( indiwo ), ifail, info )
491*
492* Apply orthogonal matrix used in reduction to tridiagonal
493* form to eigenvectors returned by DSTEIN.
494*
495 indwkn = inde
496 llwrkn = lwork - indwkn + 1
497 CALL dormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
498 $ ldz, work( indwkn ), llwrkn, iinfo )
499 END IF
500*
501* If matrix was scaled, then rescale eigenvalues appropriately.
502*
503 40 CONTINUE
504 IF( iscale.EQ.1 ) THEN
505 IF( info.EQ.0 ) THEN
506 imax = m
507 ELSE
508 imax = info - 1
509 END IF
510 CALL dscal( imax, one / sigma, w, 1 )
511 END IF
512*
513* If eigenvalues are not in order, then sort them, along with
514* eigenvectors.
515*
516 IF( wantz ) THEN
517 DO 60 j = 1, m - 1
518 i = 0
519 tmp1 = w( j )
520 DO 50 jj = j + 1, m
521 IF( w( jj ).LT.tmp1 ) THEN
522 i = jj
523 tmp1 = w( jj )
524 END IF
525 50 CONTINUE
526*
527 IF( i.NE.0 ) THEN
528 itmp1 = iwork( indibl+i-1 )
529 w( i ) = w( j )
530 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
531 w( j ) = tmp1
532 iwork( indibl+j-1 ) = itmp1
533 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
534 IF( info.NE.0 ) THEN
535 itmp1 = ifail( i )
536 ifail( i ) = ifail( j )
537 ifail( j ) = itmp1
538 END IF
539 END IF
540 60 CONTINUE
541 END IF
542*
543* Set WORK(1) to optimal workspace size.
544*
545 work( 1 ) = lwkopt
546*
547 RETURN
548*
549* End of DSYEVX
550*

◆ dsyevx_2stage()

subroutine dsyevx_2stage ( character jobz,
character range,
character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

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

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

Purpose:
!>
!> DSYEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A using the 2stage technique for
!> the reduction to tridiagonal.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of indices
!> for the desired eigenvalues.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>                  Not available in this release.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[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, 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).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[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 (MAX(1,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, 8*N, dimension) where
!>                                   dimension = max(stage1,stage2) + (KD+1)*N + 3*N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 3*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]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, then i eigenvectors failed to converge.
!>                Their indices are stored in array IFAIL.
!>
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 297 of file dsyevx_2stage.f.

300*
301 IMPLICIT NONE
302*
303* -- LAPACK driver routine --
304* -- LAPACK is a software package provided by Univ. of Tennessee, --
305* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
306*
307* .. Scalar Arguments ..
308 CHARACTER JOBZ, RANGE, UPLO
309 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
310 DOUBLE PRECISION ABSTOL, VL, VU
311* ..
312* .. Array Arguments ..
313 INTEGER IFAIL( * ), IWORK( * )
314 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
315* ..
316*
317* =====================================================================
318*
319* .. Parameters ..
320 DOUBLE PRECISION ZERO, ONE
321 parameter( zero = 0.0d+0, one = 1.0d+0 )
322* ..
323* .. Local Scalars ..
324 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
325 $ WANTZ
326 CHARACTER ORDER
327 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
328 $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
329 $ ITMP1, J, JJ, LLWORK, LLWRKN,
330 $ NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
331 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
332 $ SIGMA, SMLNUM, TMP1, VLL, VUU
333* ..
334* .. External Functions ..
335 LOGICAL LSAME
336 INTEGER ILAENV2STAGE
337 DOUBLE PRECISION DLAMCH, DLANSY
338 EXTERNAL lsame, dlamch, dlansy, ilaenv2stage
339* ..
340* .. External Subroutines ..
341 EXTERNAL dcopy, dlacpy, dorgtr, dormtr, dscal, dstebz,
342 $ dstein, dsteqr, dsterf, dswap, xerbla,
343 $ dsytrd_2stage
344* ..
345* .. Intrinsic Functions ..
346 INTRINSIC max, min, sqrt
347* ..
348* .. Executable Statements ..
349*
350* Test the input parameters.
351*
352 lower = lsame( uplo, 'L' )
353 wantz = lsame( jobz, 'V' )
354 alleig = lsame( range, 'A' )
355 valeig = lsame( range, 'V' )
356 indeig = lsame( range, 'I' )
357 lquery = ( lwork.EQ.-1 )
358*
359 info = 0
360 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
361 info = -1
362 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
363 info = -2
364 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
365 info = -3
366 ELSE IF( n.LT.0 ) THEN
367 info = -4
368 ELSE IF( lda.LT.max( 1, n ) ) THEN
369 info = -6
370 ELSE
371 IF( valeig ) THEN
372 IF( n.GT.0 .AND. vu.LE.vl )
373 $ info = -8
374 ELSE IF( indeig ) THEN
375 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
376 info = -9
377 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
378 info = -10
379 END IF
380 END IF
381 END IF
382 IF( info.EQ.0 ) THEN
383 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
384 info = -15
385 END IF
386 END IF
387*
388 IF( info.EQ.0 ) THEN
389 IF( n.LE.1 ) THEN
390 lwmin = 1
391 work( 1 ) = lwmin
392 ELSE
393 kd = ilaenv2stage( 1, 'DSYTRD_2STAGE', jobz,
394 $ n, -1, -1, -1 )
395 ib = ilaenv2stage( 2, 'DSYTRD_2STAGE', jobz,
396 $ n, kd, -1, -1 )
397 lhtrd = ilaenv2stage( 3, 'DSYTRD_2STAGE', jobz,
398 $ n, kd, ib, -1 )
399 lwtrd = ilaenv2stage( 4, 'DSYTRD_2STAGE', jobz,
400 $ n, kd, ib, -1 )
401 lwmin = max( 8*n, 3*n + lhtrd + lwtrd )
402 work( 1 ) = lwmin
403 END IF
404*
405 IF( lwork.LT.lwmin .AND. .NOT.lquery )
406 $ info = -17
407 END IF
408*
409 IF( info.NE.0 ) THEN
410 CALL xerbla( 'DSYEVX_2STAGE', -info )
411 RETURN
412 ELSE IF( lquery ) THEN
413 RETURN
414 END IF
415*
416* Quick return if possible
417*
418 m = 0
419 IF( n.EQ.0 ) THEN
420 RETURN
421 END IF
422*
423 IF( n.EQ.1 ) THEN
424 IF( alleig .OR. indeig ) THEN
425 m = 1
426 w( 1 ) = a( 1, 1 )
427 ELSE
428 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
429 m = 1
430 w( 1 ) = a( 1, 1 )
431 END IF
432 END IF
433 IF( wantz )
434 $ z( 1, 1 ) = one
435 RETURN
436 END IF
437*
438* Get machine constants.
439*
440 safmin = dlamch( 'Safe minimum' )
441 eps = dlamch( 'Precision' )
442 smlnum = safmin / eps
443 bignum = one / smlnum
444 rmin = sqrt( smlnum )
445 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
446*
447* Scale matrix to allowable range, if necessary.
448*
449 iscale = 0
450 abstll = abstol
451 IF( valeig ) THEN
452 vll = vl
453 vuu = vu
454 END IF
455 anrm = dlansy( 'M', uplo, n, a, lda, work )
456 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
457 iscale = 1
458 sigma = rmin / anrm
459 ELSE IF( anrm.GT.rmax ) THEN
460 iscale = 1
461 sigma = rmax / anrm
462 END IF
463 IF( iscale.EQ.1 ) THEN
464 IF( lower ) THEN
465 DO 10 j = 1, n
466 CALL dscal( n-j+1, sigma, a( j, j ), 1 )
467 10 CONTINUE
468 ELSE
469 DO 20 j = 1, n
470 CALL dscal( j, sigma, a( 1, j ), 1 )
471 20 CONTINUE
472 END IF
473 IF( abstol.GT.0 )
474 $ abstll = abstol*sigma
475 IF( valeig ) THEN
476 vll = vl*sigma
477 vuu = vu*sigma
478 END IF
479 END IF
480*
481* Call DSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
482*
483 indtau = 1
484 inde = indtau + n
485 indd = inde + n
486 indhous = indd + n
487 indwrk = indhous + lhtrd
488 llwork = lwork - indwrk + 1
489*
490 CALL dsytrd_2stage( jobz, uplo, n, a, lda, work( indd ),
491 $ work( inde ), work( indtau ), work( indhous ),
492 $ lhtrd, work( indwrk ), llwork, iinfo )
493*
494* If all eigenvalues are desired and ABSTOL is less than or equal to
495* zero, then call DSTERF or DORGTR and SSTEQR. If this fails for
496* some eigenvalue, then try DSTEBZ.
497*
498 test = .false.
499 IF( indeig ) THEN
500 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
501 test = .true.
502 END IF
503 END IF
504 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
505 CALL dcopy( n, work( indd ), 1, w, 1 )
506 indee = indwrk + 2*n
507 IF( .NOT.wantz ) THEN
508 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
509 CALL dsterf( n, w, work( indee ), info )
510 ELSE
511 CALL dlacpy( 'A', n, n, a, lda, z, ldz )
512 CALL dorgtr( uplo, n, z, ldz, work( indtau ),
513 $ work( indwrk ), llwork, iinfo )
514 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
515 CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
516 $ work( indwrk ), info )
517 IF( info.EQ.0 ) THEN
518 DO 30 i = 1, n
519 ifail( i ) = 0
520 30 CONTINUE
521 END IF
522 END IF
523 IF( info.EQ.0 ) THEN
524 m = n
525 GO TO 40
526 END IF
527 info = 0
528 END IF
529*
530* Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
531*
532 IF( wantz ) THEN
533 order = 'B'
534 ELSE
535 order = 'E'
536 END IF
537 indibl = 1
538 indisp = indibl + n
539 indiwo = indisp + n
540 CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
541 $ work( indd ), work( inde ), m, nsplit, w,
542 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
543 $ iwork( indiwo ), info )
544*
545 IF( wantz ) THEN
546 CALL dstein( n, work( indd ), work( inde ), m, w,
547 $ iwork( indibl ), iwork( indisp ), z, ldz,
548 $ work( indwrk ), iwork( indiwo ), ifail, info )
549*
550* Apply orthogonal matrix used in reduction to tridiagonal
551* form to eigenvectors returned by DSTEIN.
552*
553 indwkn = inde
554 llwrkn = lwork - indwkn + 1
555 CALL dormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
556 $ ldz, work( indwkn ), llwrkn, iinfo )
557 END IF
558*
559* If matrix was scaled, then rescale eigenvalues appropriately.
560*
561 40 CONTINUE
562 IF( iscale.EQ.1 ) THEN
563 IF( info.EQ.0 ) THEN
564 imax = m
565 ELSE
566 imax = info - 1
567 END IF
568 CALL dscal( imax, one / sigma, w, 1 )
569 END IF
570*
571* If eigenvalues are not in order, then sort them, along with
572* eigenvectors.
573*
574 IF( wantz ) THEN
575 DO 60 j = 1, m - 1
576 i = 0
577 tmp1 = w( j )
578 DO 50 jj = j + 1, m
579 IF( w( jj ).LT.tmp1 ) THEN
580 i = jj
581 tmp1 = w( jj )
582 END IF
583 50 CONTINUE
584*
585 IF( i.NE.0 ) THEN
586 itmp1 = iwork( indibl+i-1 )
587 w( i ) = w( j )
588 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
589 w( j ) = tmp1
590 iwork( indibl+j-1 ) = itmp1
591 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
592 IF( info.NE.0 ) THEN
593 itmp1 = ifail( i )
594 ifail( i ) = ifail( j )
595 ifail( j ) = itmp1
596 END IF
597 END IF
598 60 CONTINUE
599 END IF
600*
601* Set WORK(1) to optimal workspace size.
602*
603 work( 1 ) = lwmin
604*
605 RETURN
606*
607* End of DSYEVX_2STAGE
608*

◆ dsygv()

subroutine dsygv ( integer itype,
character jobz,
character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) w,
double precision, dimension( * ) work,
integer lwork,
integer info )

DSYGV

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

Purpose:
!>
!> DSYGV computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
!> Here A and B are assumed to be symmetric and B is also
!> positive definite.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  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
!>          matrix Z of eigenvectors.  The eigenvectors are normalized
!>          as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
!>          or the lower triangle (if UPLO='L') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB, N)
!>          On entry, the symmetric positive definite matrix B.
!>          If UPLO = 'U', the leading N-by-N upper triangular part of B
!>          contains the upper triangular part of the matrix B.
!>          If UPLO = 'L', the leading N-by-N lower triangular part of B
!>          contains the lower triangular part of the matrix B.
!>
!>          On exit, if INFO <= N, the part of B containing the matrix is
!>          overwritten by the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= 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 (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK.  LWORK >= max(1,3*N-1).
!>          For optimal efficiency, LWORK >= (NB+2)*N,
!>          where NB is the blocksize for DSYTRD returned by ILAENV.
!>
!>          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:  DPOTRF or DSYEV returned an error code:
!>             <= N:  if INFO = i, DSYEV failed to converge;
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not converge to zero;
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 173 of file dsygv.f.

175*
176* -- LAPACK driver routine --
177* -- LAPACK is a software package provided by Univ. of Tennessee, --
178* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
179*
180* .. Scalar Arguments ..
181 CHARACTER JOBZ, UPLO
182 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
183* ..
184* .. Array Arguments ..
185 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 DOUBLE PRECISION ONE
192 parameter( one = 1.0d+0 )
193* ..
194* .. Local Scalars ..
195 LOGICAL LQUERY, UPPER, WANTZ
196 CHARACTER TRANS
197 INTEGER LWKMIN, LWKOPT, NB, NEIG
198* ..
199* .. External Functions ..
200 LOGICAL LSAME
201 INTEGER ILAENV
202 EXTERNAL lsame, ilaenv
203* ..
204* .. External Subroutines ..
205 EXTERNAL dpotrf, dsyev, dsygst, dtrmm, dtrsm, xerbla
206* ..
207* .. Intrinsic Functions ..
208 INTRINSIC max
209* ..
210* .. Executable Statements ..
211*
212* Test the input parameters.
213*
214 wantz = lsame( jobz, 'V' )
215 upper = lsame( uplo, 'U' )
216 lquery = ( lwork.EQ.-1 )
217*
218 info = 0
219 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
220 info = -1
221 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
222 info = -2
223 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
224 info = -3
225 ELSE IF( n.LT.0 ) THEN
226 info = -4
227 ELSE IF( lda.LT.max( 1, n ) ) THEN
228 info = -6
229 ELSE IF( ldb.LT.max( 1, n ) ) THEN
230 info = -8
231 END IF
232*
233 IF( info.EQ.0 ) THEN
234 lwkmin = max( 1, 3*n - 1 )
235 nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
236 lwkopt = max( lwkmin, ( nb + 2 )*n )
237 work( 1 ) = lwkopt
238*
239 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
240 info = -11
241 END IF
242 END IF
243*
244 IF( info.NE.0 ) THEN
245 CALL xerbla( 'DSYGV ', -info )
246 RETURN
247 ELSE IF( lquery ) THEN
248 RETURN
249 END IF
250*
251* Quick return if possible
252*
253 IF( n.EQ.0 )
254 $ RETURN
255*
256* Form a Cholesky factorization of B.
257*
258 CALL dpotrf( uplo, n, b, ldb, info )
259 IF( info.NE.0 ) THEN
260 info = n + info
261 RETURN
262 END IF
263*
264* Transform problem to standard eigenvalue problem and solve.
265*
266 CALL dsygst( itype, uplo, n, a, lda, b, ldb, info )
267 CALL dsyev( jobz, uplo, n, a, lda, w, work, lwork, info )
268*
269 IF( wantz ) THEN
270*
271* Backtransform eigenvectors to the original problem.
272*
273 neig = n
274 IF( info.GT.0 )
275 $ neig = info - 1
276 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
277*
278* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
279* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
280*
281 IF( upper ) THEN
282 trans = 'N'
283 ELSE
284 trans = 'T'
285 END IF
286*
287 CALL dtrsm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
288 $ b, ldb, a, lda )
289*
290 ELSE IF( itype.EQ.3 ) THEN
291*
292* For B*A*x=(lambda)*x;
293* backtransform eigenvectors: x = L*y or U**T*y
294*
295 IF( upper ) THEN
296 trans = 'T'
297 ELSE
298 trans = 'N'
299 END IF
300*
301 CALL dtrmm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
302 $ b, ldb, a, lda )
303 END IF
304 END IF
305*
306 work( 1 ) = lwkopt
307 RETURN
308*
309* End of DSYGV
310*
dpotrf
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
Definition dpotrf.f:107
dsygst
subroutine dsygst(itype, uplo, n, a, lda, b, ldb, info)
DSYGST
Definition dsygst.f:127
dsyev
subroutine dsyev(jobz, uplo, n, a, lda, w, work, lwork, info)
DSYEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition dsyev.f:132
dtrmm
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181

◆ dsygv_2stage()

subroutine dsygv_2stage ( integer itype,
character jobz,
character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) w,
double precision, dimension( * ) work,
integer lwork,
integer info )

DSYGV_2STAGE

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

Purpose:
!>
!> DSYGV_2STAGE computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
!> Here A and B are assumed to be symmetric and B is also
!> positive definite.
!> This routine use the 2stage technique for the reduction to tridiagonal
!> which showed higher performance on recent architecture and for large
!> sizes N>2000.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[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 triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  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
!>          matrix Z of eigenvectors.  The eigenvectors are normalized
!>          as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
!>          or the lower triangle (if UPLO='L') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB, N)
!>          On entry, the symmetric positive definite matrix B.
!>          If UPLO = 'U', the leading N-by-N upper triangular part of B
!>          contains the upper triangular part of the matrix B.
!>          If UPLO = 'L', the leading N-by-N lower triangular part of B
!>          contains the lower triangular part of the matrix B.
!>
!>          On exit, if INFO <= N, the part of B containing the matrix is
!>          overwritten by the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= 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 (MAX(1,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:  DPOTRF or DSYEV returned an error code:
!>             <= N:  if INFO = i, DSYEV failed to converge;
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not converge to zero;
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
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 224 of file dsygv_2stage.f.

226*
227 IMPLICIT NONE
228*
229* -- LAPACK driver routine --
230* -- LAPACK is a software package provided by Univ. of Tennessee, --
231* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232*
233* .. Scalar Arguments ..
234 CHARACTER JOBZ, UPLO
235 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
236* ..
237* .. Array Arguments ..
238 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 DOUBLE PRECISION ONE
245 parameter( one = 1.0d+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL LQUERY, UPPER, WANTZ
249 CHARACTER TRANS
250 INTEGER NEIG, LWMIN, LHTRD, LWTRD, KD, IB
251* ..
252* .. External Functions ..
253 LOGICAL LSAME
254 INTEGER ILAENV2STAGE
255 EXTERNAL lsame, ilaenv2stage
256* ..
257* .. External Subroutines ..
258 EXTERNAL dpotrf, dsygst, dtrmm, dtrsm, xerbla,
259 $ dsyev_2stage
260* ..
261* .. Intrinsic Functions ..
262 INTRINSIC max
263* ..
264* .. Executable Statements ..
265*
266* Test the input parameters.
267*
268 wantz = lsame( jobz, 'V' )
269 upper = lsame( uplo, 'U' )
270 lquery = ( lwork.EQ.-1 )
271*
272 info = 0
273 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
274 info = -1
275 ELSE IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
276 info = -2
277 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
278 info = -3
279 ELSE IF( n.LT.0 ) THEN
280 info = -4
281 ELSE IF( lda.LT.max( 1, n ) ) THEN
282 info = -6
283 ELSE IF( ldb.LT.max( 1, n ) ) THEN
284 info = -8
285 END IF
286*
287 IF( info.EQ.0 ) THEN
288 kd = ilaenv2stage( 1, 'DSYTRD_2STAGE', jobz, n, -1, -1, -1 )
289 ib = ilaenv2stage( 2, 'DSYTRD_2STAGE', jobz, n, kd, -1, -1 )
290 lhtrd = ilaenv2stage( 3, 'DSYTRD_2STAGE', jobz, n, kd, ib, -1 )
291 lwtrd = ilaenv2stage( 4, 'DSYTRD_2STAGE', jobz, n, kd, ib, -1 )
292 lwmin = 2*n + lhtrd + lwtrd
293 work( 1 ) = lwmin
294*
295 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
296 info = -11
297 END IF
298 END IF
299*
300 IF( info.NE.0 ) THEN
301 CALL xerbla( 'DSYGV_2STAGE ', -info )
302 RETURN
303 ELSE IF( lquery ) THEN
304 RETURN
305 END IF
306*
307* Quick return if possible
308*
309 IF( n.EQ.0 )
310 $ RETURN
311*
312* Form a Cholesky factorization of B.
313*
314 CALL dpotrf( uplo, n, b, ldb, info )
315 IF( info.NE.0 ) THEN
316 info = n + info
317 RETURN
318 END IF
319*
320* Transform problem to standard eigenvalue problem and solve.
321*
322 CALL dsygst( itype, uplo, n, a, lda, b, ldb, info )
323 CALL dsyev_2stage( jobz, uplo, n, a, lda, w, work, lwork, info )
324*
325 IF( wantz ) THEN
326*
327* Backtransform eigenvectors to the original problem.
328*
329 neig = n
330 IF( info.GT.0 )
331 $ neig = info - 1
332 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
333*
334* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
335* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
336*
337 IF( upper ) THEN
338 trans = 'N'
339 ELSE
340 trans = 'T'
341 END IF
342*
343 CALL dtrsm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
344 $ b, ldb, a, lda )
345*
346 ELSE IF( itype.EQ.3 ) THEN
347*
348* For B*A*x=(lambda)*x;
349* backtransform eigenvectors: x = L*y or U**T*y
350*
351 IF( upper ) THEN
352 trans = 'T'
353 ELSE
354 trans = 'N'
355 END IF
356*
357 CALL dtrmm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
358 $ b, ldb, a, lda )
359 END IF
360 END IF
361*
362 work( 1 ) = lwmin
363 RETURN
364*
365* End of DSYGV_2STAGE
366*
dsyev_2stage
subroutine dsyev_2stage(jobz, uplo, n, a, lda, w, work, lwork, info)
DSYEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matr...
Definition dsyev_2stage.f:183

◆ dsygvd()

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

DSYGVD

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

Purpose:
!>
!> DSYGVD computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
!> B are assumed to be symmetric and B is also positive definite.
!> 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]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  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
!>          matrix Z of eigenvectors.  The eigenvectors are normalized
!>          as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
!>          or the lower triangle (if UPLO='L') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB, N)
!>          On entry, the symmetric matrix B.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of B contains the
!>          upper triangular part of the matrix B.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of B contains
!>          the lower triangular part of the matrix B.
!>
!>          On exit, if INFO <= N, the part of B containing the matrix is
!>          overwritten by the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= 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 (MAX(1,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 >= 1.
!>          If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
!>          If JOBZ = 'V' and N > 1, LWORK >= 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 >= 1.
!>          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
!>          If JOBZ  = 'V' and N > 1, LIWORK >= 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:  DPOTRF or DSYEVD returned an error code:
!>             <= N:  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);
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Modified so that no backsubstitution is performed if DSYEVD fails to
!>  converge (NEIG in old code could be greater than N causing out of
!>  bounds reference to A - reported by Ralf Meyer).  Also corrected the
!>  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
!>
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 225 of file dsygvd.f.

227*
228* -- LAPACK driver routine --
229* -- LAPACK is a software package provided by Univ. of Tennessee, --
230* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231*
232* .. Scalar Arguments ..
233 CHARACTER JOBZ, UPLO
234 INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
235* ..
236* .. Array Arguments ..
237 INTEGER IWORK( * )
238 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 DOUBLE PRECISION ONE
245 parameter( one = 1.0d+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL LQUERY, UPPER, WANTZ
249 CHARACTER TRANS
250 INTEGER LIOPT, LIWMIN, LOPT, LWMIN
251* ..
252* .. External Functions ..
253 LOGICAL LSAME
254 EXTERNAL lsame
255* ..
256* .. External Subroutines ..
257 EXTERNAL dpotrf, dsyevd, dsygst, dtrmm, dtrsm, xerbla
258* ..
259* .. Intrinsic Functions ..
260 INTRINSIC dble, max
261* ..
262* .. Executable Statements ..
263*
264* Test the input parameters.
265*
266 wantz = lsame( jobz, 'V' )
267 upper = lsame( uplo, 'U' )
268 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
269*
270 info = 0
271 IF( n.LE.1 ) THEN
272 liwmin = 1
273 lwmin = 1
274 ELSE IF( wantz ) THEN
275 liwmin = 3 + 5*n
276 lwmin = 1 + 6*n + 2*n**2
277 ELSE
278 liwmin = 1
279 lwmin = 2*n + 1
280 END IF
281 lopt = lwmin
282 liopt = liwmin
283 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
284 info = -1
285 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
286 info = -2
287 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
288 info = -3
289 ELSE IF( n.LT.0 ) THEN
290 info = -4
291 ELSE IF( lda.LT.max( 1, n ) ) THEN
292 info = -6
293 ELSE IF( ldb.LT.max( 1, n ) ) THEN
294 info = -8
295 END IF
296*
297 IF( info.EQ.0 ) THEN
298 work( 1 ) = lopt
299 iwork( 1 ) = liopt
300*
301 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
302 info = -11
303 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
304 info = -13
305 END IF
306 END IF
307*
308 IF( info.NE.0 ) THEN
309 CALL xerbla( 'DSYGVD', -info )
310 RETURN
311 ELSE IF( lquery ) THEN
312 RETURN
313 END IF
314*
315* Quick return if possible
316*
317 IF( n.EQ.0 )
318 $ RETURN
319*
320* Form a Cholesky factorization of B.
321*
322 CALL dpotrf( uplo, n, b, ldb, info )
323 IF( info.NE.0 ) THEN
324 info = n + info
325 RETURN
326 END IF
327*
328* Transform problem to standard eigenvalue problem and solve.
329*
330 CALL dsygst( itype, uplo, n, a, lda, b, ldb, info )
331 CALL dsyevd( jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork,
332 $ info )
333 lopt = max( dble( lopt ), dble( work( 1 ) ) )
334 liopt = max( dble( liopt ), dble( iwork( 1 ) ) )
335*
336 IF( wantz .AND. info.EQ.0 ) THEN
337*
338* Backtransform eigenvectors to the original problem.
339*
340 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
341*
342* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
343* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
344*
345 IF( upper ) THEN
346 trans = 'N'
347 ELSE
348 trans = 'T'
349 END IF
350*
351 CALL dtrsm( 'Left', uplo, trans, 'Non-unit', n, n, one,
352 $ b, ldb, a, lda )
353*
354 ELSE IF( itype.EQ.3 ) THEN
355*
356* For B*A*x=(lambda)*x;
357* backtransform eigenvectors: x = L*y or U**T*y
358*
359 IF( upper ) THEN
360 trans = 'T'
361 ELSE
362 trans = 'N'
363 END IF
364*
365 CALL dtrmm( 'Left', uplo, trans, 'Non-unit', n, n, one,
366 $ b, ldb, a, lda )
367 END IF
368 END IF
369*
370 work( 1 ) = lopt
371 iwork( 1 ) = liopt
372*
373 RETURN
374*
375* End of DSYGVD
376*
dsyevd
subroutine dsyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
DSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition dsyevd.f:185

◆ dsygvx()

subroutine dsygvx ( integer itype,
character jobz,
character range,
character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

DSYGVX

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

Purpose:
!>
!> DSYGVX computes selected eigenvalues, and optionally, eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
!> and B are assumed to be symmetric and B is also positive definite.
!> Eigenvalues and eigenvectors can be selected by specifying either a
!> range of values or a range of indices for the desired eigenvalues.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A and B are stored;
!>          = 'L':  Lower triangle of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix pencil (A,B).  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, 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).
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB, N)
!>          On entry, the symmetric matrix B.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of B contains the
!>          upper triangular part of the matrix B.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of B contains
!>          the lower triangular part of the matrix B.
!>
!>          On exit, if INFO <= N, the part of B containing the matrix is
!>          overwritten by the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing C to tridiagonal form, where C is the symmetric
!>          matrix of the standard symmetric problem to which the
!>          generalized problem is transformed.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'N', then Z is not referenced.
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          The eigenvectors are normalized as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[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 (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK.  LWORK >= max(1,8*N).
!>          For optimal efficiency, LWORK >= (NB+3)*N,
!>          where NB is the blocksize for DSYTRD returned by ILAENV.
!>
!>          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]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  DPOTRF or DSYEVX returned an error code:
!>             <= N:  if INFO = i, DSYEVX failed to converge;
!>                    i eigenvectors failed to converge.  Their indices
!>                    are stored in array IFAIL.
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 294 of file dsygvx.f.

297*
298* -- LAPACK driver routine --
299* -- LAPACK is a software package provided by Univ. of Tennessee, --
300* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
301*
302* .. Scalar Arguments ..
303 CHARACTER JOBZ, RANGE, UPLO
304 INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
305 DOUBLE PRECISION ABSTOL, VL, VU
306* ..
307* .. Array Arguments ..
308 INTEGER IFAIL( * ), IWORK( * )
309 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ),
310 $ Z( LDZ, * )
311* ..
312*
313* =====================================================================
314*
315* .. Parameters ..
316 DOUBLE PRECISION ONE
317 parameter( one = 1.0d+0 )
318* ..
319* .. Local Scalars ..
320 LOGICAL ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ
321 CHARACTER TRANS
322 INTEGER LWKMIN, LWKOPT, NB
323* ..
324* .. External Functions ..
325 LOGICAL LSAME
326 INTEGER ILAENV
327 EXTERNAL lsame, ilaenv
328* ..
329* .. External Subroutines ..
330 EXTERNAL dpotrf, dsyevx, dsygst, dtrmm, dtrsm, xerbla
331* ..
332* .. Intrinsic Functions ..
333 INTRINSIC max, min
334* ..
335* .. Executable Statements ..
336*
337* Test the input parameters.
338*
339 upper = lsame( uplo, 'U' )
340 wantz = lsame( jobz, 'V' )
341 alleig = lsame( range, 'A' )
342 valeig = lsame( range, 'V' )
343 indeig = lsame( range, 'I' )
344 lquery = ( lwork.EQ.-1 )
345*
346 info = 0
347 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
348 info = -1
349 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
350 info = -2
351 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
352 info = -3
353 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
354 info = -4
355 ELSE IF( n.LT.0 ) THEN
356 info = -5
357 ELSE IF( lda.LT.max( 1, n ) ) THEN
358 info = -7
359 ELSE IF( ldb.LT.max( 1, n ) ) THEN
360 info = -9
361 ELSE
362 IF( valeig ) THEN
363 IF( n.GT.0 .AND. vu.LE.vl )
364 $ info = -11
365 ELSE IF( indeig ) THEN
366 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
367 info = -12
368 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
369 info = -13
370 END IF
371 END IF
372 END IF
373 IF (info.EQ.0) THEN
374 IF (ldz.LT.1 .OR. (wantz .AND. ldz.LT.n)) THEN
375 info = -18
376 END IF
377 END IF
378*
379 IF( info.EQ.0 ) THEN
380 lwkmin = max( 1, 8*n )
381 nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
382 lwkopt = max( lwkmin, ( nb + 3 )*n )
383 work( 1 ) = lwkopt
384*
385 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
386 info = -20
387 END IF
388 END IF
389*
390 IF( info.NE.0 ) THEN
391 CALL xerbla( 'DSYGVX', -info )
392 RETURN
393 ELSE IF( lquery ) THEN
394 RETURN
395 END IF
396*
397* Quick return if possible
398*
399 m = 0
400 IF( n.EQ.0 ) THEN
401 RETURN
402 END IF
403*
404* Form a Cholesky factorization of B.
405*
406 CALL dpotrf( uplo, n, b, ldb, info )
407 IF( info.NE.0 ) THEN
408 info = n + info
409 RETURN
410 END IF
411*
412* Transform problem to standard eigenvalue problem and solve.
413*
414 CALL dsygst( itype, uplo, n, a, lda, b, ldb, info )
415 CALL dsyevx( jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol,
416 $ m, w, z, ldz, work, lwork, iwork, ifail, info )
417*
418 IF( wantz ) THEN
419*
420* Backtransform eigenvectors to the original problem.
421*
422 IF( info.GT.0 )
423 $ m = info - 1
424 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
425*
426* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
427* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
428*
429 IF( upper ) THEN
430 trans = 'N'
431 ELSE
432 trans = 'T'
433 END IF
434*
435 CALL dtrsm( 'Left', uplo, trans, 'Non-unit', n, m, one, b,
436 $ ldb, z, ldz )
437*
438 ELSE IF( itype.EQ.3 ) THEN
439*
440* For B*A*x=(lambda)*x;
441* backtransform eigenvectors: x = L*y or U**T*y
442*
443 IF( upper ) THEN
444 trans = 'T'
445 ELSE
446 trans = 'N'
447 END IF
448*
449 CALL dtrmm( 'Left', uplo, trans, 'Non-unit', n, m, one, b,
450 $ ldb, z, ldz )
451 END IF
452 END IF
453*
454* Set WORK(1) to optimal workspace size.
455*
456 work( 1 ) = lwkopt
457*
458 RETURN
459*
460* End of DSYGVX
461*
dsyevx
subroutine dsyevx(jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
DSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition dsyevx.f:253