OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
complex
LAPACK » Hermitian Matrix » Eigenvalue

Functions

subroutine cheev (jobz, uplo, n, a, lda, w, work, lwork, rwork, info)
  CHEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
subroutine cheev_2stage (jobz, uplo, n, a, lda, w, work, lwork, rwork, info)
  CHEEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
subroutine cheevd (jobz, uplo, n, a, lda, w, work, lwork, rwork, lrwork, iwork, liwork, info)
  CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
subroutine cheevd_2stage (jobz, uplo, n, a, lda, w, work, lwork, rwork, lrwork, iwork, liwork, info)
  CHEEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
subroutine cheevr (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, rwork, lrwork, iwork, liwork, info)
  CHEEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
subroutine cheevr_2stage (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, rwork, lrwork, iwork, liwork, info)
  CHEEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
subroutine cheevx (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, rwork, iwork, ifail, info)
  CHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
subroutine cheevx_2stage (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, rwork, iwork, ifail, info)
  CHEEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
subroutine chegv (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, rwork, info)
 CHEGV
subroutine chegv_2stage (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, rwork, info)
 CHEGV_2STAGE
subroutine chegvd (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, rwork, lrwork, iwork, liwork, info)
 CHEGVD
subroutine chegvx (itype, jobz, range, uplo, n, a, lda, b, ldb, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, rwork, iwork, ifail, info)
 CHEGVX

Detailed Description

This is the group of complex eigenvalue driver functions for HE matrices

Function Documentation

◆ cheev()

subroutine cheev ( character jobz,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real, dimension( * ) w,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

CHEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
!>
!> CHEEV computes all eigenvalues and, optionally, eigenvectors of a
!> complex Hermitian 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is COMPLEX 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,2*N-1).
!>          For optimal efficiency, LWORK >= (NB+1)*N,
!>          where NB is the blocksize for CHETRD 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]RWORK
!>          RWORK is REAL array, dimension (max(1, 3*N-2))
!>
[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 138 of file cheev.f.

140*
141* -- LAPACK driver routine --
142* -- LAPACK is a software package provided by Univ. of Tennessee, --
143* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144*
145* .. Scalar Arguments ..
146 CHARACTER JOBZ, UPLO
147 INTEGER INFO, LDA, LWORK, N
148* ..
149* .. Array Arguments ..
150 REAL RWORK( * ), W( * )
151 COMPLEX A( LDA, * ), WORK( * )
152* ..
153*
154* =====================================================================
155*
156* .. Parameters ..
157 REAL ZERO, ONE
158 parameter( zero = 0.0e0, one = 1.0e0 )
159 COMPLEX CONE
160 parameter( cone = ( 1.0e0, 0.0e0 ) )
161* ..
162* .. Local Scalars ..
163 LOGICAL LOWER, LQUERY, WANTZ
164 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
165 $ LLWORK, LWKOPT, NB
166 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
167 $ SMLNUM
168* ..
169* .. External Functions ..
170 LOGICAL LSAME
171 INTEGER ILAENV
172 REAL CLANHE, SLAMCH
173 EXTERNAL ilaenv, lsame, clanhe, slamch
174* ..
175* .. External Subroutines ..
176 EXTERNAL chetrd, clascl, csteqr, cungtr, sscal, ssterf,
177 $ xerbla
178* ..
179* .. Intrinsic Functions ..
180 INTRINSIC max, sqrt
181* ..
182* .. Executable Statements ..
183*
184* Test the input parameters.
185*
186 wantz = lsame( jobz, 'V' )
187 lower = lsame( uplo, 'L' )
188 lquery = ( lwork.EQ.-1 )
189*
190 info = 0
191 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
192 info = -1
193 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
194 info = -2
195 ELSE IF( n.LT.0 ) THEN
196 info = -3
197 ELSE IF( lda.LT.max( 1, n ) ) THEN
198 info = -5
199 END IF
200*
201 IF( info.EQ.0 ) THEN
202 nb = ilaenv( 1, 'CHETRD', uplo, n, -1, -1, -1 )
203 lwkopt = max( 1, ( nb+1 )*n )
204 work( 1 ) = lwkopt
205*
206 IF( lwork.LT.max( 1, 2*n-1 ) .AND. .NOT.lquery )
207 $ info = -8
208 END IF
209*
210 IF( info.NE.0 ) THEN
211 CALL xerbla( 'CHEEV ', -info )
212 RETURN
213 ELSE IF( lquery ) THEN
214 RETURN
215 END IF
216*
217* Quick return if possible
218*
219 IF( n.EQ.0 ) THEN
220 RETURN
221 END IF
222*
223 IF( n.EQ.1 ) THEN
224 w( 1 ) = real( a( 1, 1 ) )
225 work( 1 ) = 1
226 IF( wantz )
227 $ a( 1, 1 ) = cone
228 RETURN
229 END IF
230*
231* Get machine constants.
232*
233 safmin = slamch( 'Safe minimum' )
234 eps = slamch( 'Precision' )
235 smlnum = safmin / eps
236 bignum = one / smlnum
237 rmin = sqrt( smlnum )
238 rmax = sqrt( bignum )
239*
240* Scale matrix to allowable range, if necessary.
241*
242 anrm = clanhe( 'M', uplo, n, a, lda, rwork )
243 iscale = 0
244 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
245 iscale = 1
246 sigma = rmin / anrm
247 ELSE IF( anrm.GT.rmax ) THEN
248 iscale = 1
249 sigma = rmax / anrm
250 END IF
251 IF( iscale.EQ.1 )
252 $ CALL clascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
253*
254* Call CHETRD to reduce Hermitian matrix to tridiagonal form.
255*
256 inde = 1
257 indtau = 1
258 indwrk = indtau + n
259 llwork = lwork - indwrk + 1
260 CALL chetrd( uplo, n, a, lda, w, rwork( inde ), work( indtau ),
261 $ work( indwrk ), llwork, iinfo )
262*
263* For eigenvalues only, call SSTERF. For eigenvectors, first call
264* CUNGTR to generate the unitary matrix, then call CSTEQR.
265*
266 IF( .NOT.wantz ) THEN
267 CALL ssterf( n, w, rwork( inde ), info )
268 ELSE
269 CALL cungtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
270 $ llwork, iinfo )
271 indwrk = inde + n
272 CALL csteqr( jobz, n, w, rwork( inde ), a, lda,
273 $ rwork( indwrk ), info )
274 END IF
275*
276* If matrix was scaled, then rescale eigenvalues appropriately.
277*
278 IF( iscale.EQ.1 ) THEN
279 IF( info.EQ.0 ) THEN
280 imax = n
281 ELSE
282 imax = info - 1
283 END IF
284 CALL sscal( imax, one / sigma, w, 1 )
285 END IF
286*
287* Set WORK(1) to optimal complex workspace size.
288*
289 work( 1 ) = lwkopt
290*
291 RETURN
292*
293* End of CHEEV
294*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
clanhe
real function clanhe(norm, uplo, n, a, lda, work)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhe.f:124
chetrd
subroutine chetrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
CHETRD
Definition chetrd.f:192
clascl
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
cungtr
subroutine cungtr(uplo, n, a, lda, tau, work, lwork, info)
CUNGTR
Definition cungtr.f:123
csteqr
subroutine csteqr(compz, n, d, e, z, ldz, work, info)
CSTEQR
Definition csteqr.f:132
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
max
#define max(a, b)
Definition macros.h:21

◆ cheev_2stage()

subroutine cheev_2stage ( character jobz,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real, dimension( * ) w,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

CHEEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
!>
!> CHEEV_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
!> complex Hermitian 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is COMPLEX 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 + N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 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]RWORK
!>          RWORK is REAL array, dimension (max(1, 3*N-2))
!>
[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 187 of file cheev_2stage.f.

189*
190 IMPLICIT NONE
191*
192* -- LAPACK driver routine --
193* -- LAPACK is a software package provided by Univ. of Tennessee, --
194* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195*
196* .. Scalar Arguments ..
197 CHARACTER JOBZ, UPLO
198 INTEGER INFO, LDA, LWORK, N
199* ..
200* .. Array Arguments ..
201 REAL RWORK( * ), W( * )
202 COMPLEX A( LDA, * ), WORK( * )
203* ..
204*
205* =====================================================================
206*
207* .. Parameters ..
208 REAL ZERO, ONE
209 parameter( zero = 0.0e0, one = 1.0e0 )
210 COMPLEX CONE
211 parameter( cone = ( 1.0e0, 0.0e0 ) )
212* ..
213* .. Local Scalars ..
214 LOGICAL LOWER, LQUERY, WANTZ
215 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
216 $ LLWORK, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
217 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
218 $ SMLNUM
219* ..
220* .. External Functions ..
221 LOGICAL LSAME
222 INTEGER ILAENV2STAGE
223 REAL SLAMCH, CLANHE
224 EXTERNAL lsame, slamch, clanhe, ilaenv2stage
225* ..
226* .. External Subroutines ..
227 EXTERNAL sscal, ssterf, xerbla, clascl, csteqr,
228 $ cungtr, chetrd_2stage
229* ..
230* .. Intrinsic Functions ..
231 INTRINSIC real, max, sqrt
232* ..
233* .. Executable Statements ..
234*
235* Test the input parameters.
236*
237 wantz = lsame( jobz, 'V' )
238 lower = lsame( uplo, 'L' )
239 lquery = ( lwork.EQ.-1 )
240*
241 info = 0
242 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
243 info = -1
244 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
245 info = -2
246 ELSE IF( n.LT.0 ) THEN
247 info = -3
248 ELSE IF( lda.LT.max( 1, n ) ) THEN
249 info = -5
250 END IF
251*
252 IF( info.EQ.0 ) THEN
253 kd = ilaenv2stage( 1, 'CHETRD_2STAGE', jobz, n, -1, -1, -1 )
254 ib = ilaenv2stage( 2, 'CHETRD_2STAGE', jobz, n, kd, -1, -1 )
255 lhtrd = ilaenv2stage( 3, 'CHETRD_2STAGE', jobz, n, kd, ib, -1 )
256 lwtrd = ilaenv2stage( 4, 'CHETRD_2STAGE', jobz, n, kd, ib, -1 )
257 lwmin = n + lhtrd + lwtrd
258 work( 1 ) = lwmin
259*
260 IF( lwork.LT.lwmin .AND. .NOT.lquery )
261 $ info = -8
262 END IF
263*
264 IF( info.NE.0 ) THEN
265 CALL xerbla( 'CHEEV_2STAGE ', -info )
266 RETURN
267 ELSE IF( lquery ) THEN
268 RETURN
269 END IF
270*
271* Quick return if possible
272*
273 IF( n.EQ.0 ) THEN
274 RETURN
275 END IF
276*
277 IF( n.EQ.1 ) THEN
278 w( 1 ) = real( a( 1, 1 ) )
279 work( 1 ) = 1
280 IF( wantz )
281 $ a( 1, 1 ) = cone
282 RETURN
283 END IF
284*
285* Get machine constants.
286*
287 safmin = slamch( 'Safe minimum' )
288 eps = slamch( 'Precision' )
289 smlnum = safmin / eps
290 bignum = one / smlnum
291 rmin = sqrt( smlnum )
292 rmax = sqrt( bignum )
293*
294* Scale matrix to allowable range, if necessary.
295*
296 anrm = clanhe( 'M', uplo, n, a, lda, rwork )
297 iscale = 0
298 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
299 iscale = 1
300 sigma = rmin / anrm
301 ELSE IF( anrm.GT.rmax ) THEN
302 iscale = 1
303 sigma = rmax / anrm
304 END IF
305 IF( iscale.EQ.1 )
306 $ CALL clascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
307*
308* Call CHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
309*
310 inde = 1
311 indtau = 1
312 indhous = indtau + n
313 indwrk = indhous + lhtrd
314 llwork = lwork - indwrk + 1
315*
316 CALL chetrd_2stage( jobz, uplo, n, a, lda, w, rwork( inde ),
317 $ work( indtau ), work( indhous ), lhtrd,
318 $ work( indwrk ), llwork, iinfo )
319*
320* For eigenvalues only, call SSTERF. For eigenvectors, first call
321* CUNGTR to generate the unitary matrix, then call CSTEQR.
322*
323 IF( .NOT.wantz ) THEN
324 CALL ssterf( n, w, rwork( inde ), info )
325 ELSE
326 CALL cungtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
327 $ llwork, iinfo )
328 indwrk = inde + n
329 CALL csteqr( jobz, n, w, rwork( inde ), a, lda,
330 $ rwork( indwrk ), info )
331 END IF
332*
333* If matrix was scaled, then rescale eigenvalues appropriately.
334*
335 IF( iscale.EQ.1 ) THEN
336 IF( info.EQ.0 ) THEN
337 imax = n
338 ELSE
339 imax = info - 1
340 END IF
341 CALL sscal( imax, one / sigma, w, 1 )
342 END IF
343*
344* Set WORK(1) to optimal complex workspace size.
345*
346 work( 1 ) = lwmin
347*
348 RETURN
349*
350* End of CHEEV_2STAGE
351*
ilaenv2stage
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
Definition ilaenv2stage.f:149
chetrd_2stage
subroutine chetrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
CHETRD_2STAGE
Definition chetrd_2stage.f:224

◆ cheevd()

subroutine cheevd ( character jobz,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real, dimension( * ) w,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
!>
!> CHEEVD computes all eigenvalues and, optionally, eigenvectors of a
!> complex Hermitian 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.
!>
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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is COMPLEX 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.
!>          If N <= 1,                LWORK must be at least 1.
!>          If JOBZ  = 'N' and N > 1, LWORK must be at least N + 1.
!>          If JOBZ  = 'V' and N > 1, LWORK must be at least 2*N + N**2.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK, RWORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array,
!>                                         dimension (LRWORK)
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of the array RWORK.
!>          If N <= 1,                LRWORK must be at least 1.
!>          If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
!>          If JOBZ  = 'V' and N > 1, LRWORK must be at least
!>                         1 + 5*N + 2*N**2.
!>
!>          If LRWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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.
Further Details:
Modified description of INFO. Sven, 16 Feb 05.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 203 of file cheevd.f.

205*
206* -- LAPACK driver routine --
207* -- LAPACK is a software package provided by Univ. of Tennessee, --
208* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209*
210* .. Scalar Arguments ..
211 CHARACTER JOBZ, UPLO
212 INTEGER INFO, LDA, LIWORK, LRWORK, LWORK, N
213* ..
214* .. Array Arguments ..
215 INTEGER IWORK( * )
216 REAL RWORK( * ), W( * )
217 COMPLEX A( LDA, * ), WORK( * )
218* ..
219*
220* =====================================================================
221*
222* .. Parameters ..
223 REAL ZERO, ONE
224 parameter( zero = 0.0e0, one = 1.0e0 )
225 COMPLEX CONE
226 parameter( cone = ( 1.0e0, 0.0e0 ) )
227* ..
228* .. Local Scalars ..
229 LOGICAL LOWER, LQUERY, WANTZ
230 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
231 $ INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK,
232 $ LLWRK2, LOPT, LROPT, LRWMIN, LWMIN
233 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
234 $ SMLNUM
235* ..
236* .. External Functions ..
237 LOGICAL LSAME
238 INTEGER ILAENV
239 REAL CLANHE, SLAMCH
240 EXTERNAL ilaenv, lsame, clanhe, slamch
241* ..
242* .. External Subroutines ..
243 EXTERNAL chetrd, clacpy, clascl, cstedc, cunmtr, sscal,
244 $ ssterf, xerbla
245* ..
246* .. Intrinsic Functions ..
247 INTRINSIC max, sqrt
248* ..
249* .. Executable Statements ..
250*
251* Test the input parameters.
252*
253 wantz = lsame( jobz, 'V' )
254 lower = lsame( uplo, 'L' )
255 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
256*
257 info = 0
258 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
259 info = -1
260 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
261 info = -2
262 ELSE IF( n.LT.0 ) THEN
263 info = -3
264 ELSE IF( lda.LT.max( 1, n ) ) THEN
265 info = -5
266 END IF
267*
268 IF( info.EQ.0 ) THEN
269 IF( n.LE.1 ) THEN
270 lwmin = 1
271 lrwmin = 1
272 liwmin = 1
273 lopt = lwmin
274 lropt = lrwmin
275 liopt = liwmin
276 ELSE
277 IF( wantz ) THEN
278 lwmin = 2*n + n*n
279 lrwmin = 1 + 5*n + 2*n**2
280 liwmin = 3 + 5*n
281 ELSE
282 lwmin = n + 1
283 lrwmin = n
284 liwmin = 1
285 END IF
286 lopt = max( lwmin, n +
287 $ ilaenv( 1, 'CHETRD', uplo, n, -1, -1, -1 ) )
288 lropt = lrwmin
289 liopt = liwmin
290 END IF
291 work( 1 ) = lopt
292 rwork( 1 ) = lropt
293 iwork( 1 ) = liopt
294*
295 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
296 info = -8
297 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
298 info = -10
299 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
300 info = -12
301 END IF
302 END IF
303*
304 IF( info.NE.0 ) THEN
305 CALL xerbla( 'CHEEVD', -info )
306 RETURN
307 ELSE IF( lquery ) THEN
308 RETURN
309 END IF
310*
311* Quick return if possible
312*
313 IF( n.EQ.0 )
314 $ RETURN
315*
316 IF( n.EQ.1 ) THEN
317 w( 1 ) = real( a( 1, 1 ) )
318 IF( wantz )
319 $ a( 1, 1 ) = cone
320 RETURN
321 END IF
322*
323* Get machine constants.
324*
325 safmin = slamch( 'Safe minimum' )
326 eps = slamch( 'Precision' )
327 smlnum = safmin / eps
328 bignum = one / smlnum
329 rmin = sqrt( smlnum )
330 rmax = sqrt( bignum )
331*
332* Scale matrix to allowable range, if necessary.
333*
334 anrm = clanhe( 'M', uplo, n, a, lda, rwork )
335 iscale = 0
336 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
337 iscale = 1
338 sigma = rmin / anrm
339 ELSE IF( anrm.GT.rmax ) THEN
340 iscale = 1
341 sigma = rmax / anrm
342 END IF
343 IF( iscale.EQ.1 )
344 $ CALL clascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
345*
346* Call CHETRD to reduce Hermitian matrix to tridiagonal form.
347*
348 inde = 1
349 indtau = 1
350 indwrk = indtau + n
351 indrwk = inde + n
352 indwk2 = indwrk + n*n
353 llwork = lwork - indwrk + 1
354 llwrk2 = lwork - indwk2 + 1
355 llrwk = lrwork - indrwk + 1
356 CALL chetrd( uplo, n, a, lda, w, rwork( inde ), work( indtau ),
357 $ work( indwrk ), llwork, iinfo )
358*
359* For eigenvalues only, call SSTERF. For eigenvectors, first call
360* CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
361* tridiagonal matrix, then call CUNMTR to multiply it to the
362* Householder transformations represented as Householder vectors in
363* A.
364*
365 IF( .NOT.wantz ) THEN
366 CALL ssterf( n, w, rwork( inde ), info )
367 ELSE
368 CALL cstedc( 'I', n, w, rwork( inde ), work( indwrk ), n,
369 $ work( indwk2 ), llwrk2, rwork( indrwk ), llrwk,
370 $ iwork, liwork, info )
371 CALL cunmtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
372 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
373 CALL clacpy( 'A', n, n, work( indwrk ), n, a, lda )
374 END IF
375*
376* If matrix was scaled, then rescale eigenvalues appropriately.
377*
378 IF( iscale.EQ.1 ) THEN
379 IF( info.EQ.0 ) THEN
380 imax = n
381 ELSE
382 imax = info - 1
383 END IF
384 CALL sscal( imax, one / sigma, w, 1 )
385 END IF
386*
387 work( 1 ) = lopt
388 rwork( 1 ) = lropt
389 iwork( 1 ) = liopt
390*
391 RETURN
392*
393* End of CHEEVD
394*
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
cunmtr
subroutine cunmtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
CUNMTR
Definition cunmtr.f:172
cstedc
subroutine cstedc(compz, n, d, e, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
CSTEDC
Definition cstedc.f:212

◆ cheevd_2stage()

subroutine cheevd_2stage ( character jobz,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real, dimension( * ) w,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

CHEEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
!>
!> CHEEVD_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
!> complex Hermitian 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is COMPLEX 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 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 + N+1
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 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 2*N + N**2
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK, RWORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array,
!>                                         dimension (LRWORK)
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of the array RWORK.
!>          If N <= 1,                LRWORK must be at least 1.
!>          If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
!>          If JOBZ  = 'V' and N > 1, LRWORK must be at least
!>                         1 + 5*N + 2*N**2.
!>
!>          If LRWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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.
Further Details:
Modified description of INFO. Sven, 16 Feb 05.
Contributors:
Jeff Rutter, 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 251 of file cheevd_2stage.f.

253*
254 IMPLICIT NONE
255*
256* -- LAPACK driver routine --
257* -- LAPACK is a software package provided by Univ. of Tennessee, --
258* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
259*
260* .. Scalar Arguments ..
261 CHARACTER JOBZ, UPLO
262 INTEGER INFO, LDA, LIWORK, LRWORK, LWORK, N
263* ..
264* .. Array Arguments ..
265 INTEGER IWORK( * )
266 REAL RWORK( * ), W( * )
267 COMPLEX A( LDA, * ), WORK( * )
268* ..
269*
270* =====================================================================
271*
272* .. Parameters ..
273 REAL ZERO, ONE
274 parameter( zero = 0.0e0, one = 1.0e0 )
275 COMPLEX CONE
276 parameter( cone = ( 1.0e0, 0.0e0 ) )
277* ..
278* .. Local Scalars ..
279 LOGICAL LOWER, LQUERY, WANTZ
280 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
281 $ INDWRK, ISCALE, LIWMIN, LLRWK, LLWORK,
282 $ LLWRK2, LRWMIN, LWMIN,
283 $ LHTRD, LWTRD, KD, IB, INDHOUS
284
285
286 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
287 $ SMLNUM
288* ..
289* .. External Functions ..
290 LOGICAL LSAME
291 INTEGER ILAENV2STAGE
292 REAL SLAMCH, CLANHE
293 EXTERNAL lsame, slamch, clanhe, ilaenv2stage
294* ..
295* .. External Subroutines ..
296 EXTERNAL sscal, ssterf, xerbla, clacpy, clascl,
297 $ cstedc, cunmtr, chetrd_2stage
298* ..
299* .. Intrinsic Functions ..
300 INTRINSIC real, max, sqrt
301* ..
302* .. Executable Statements ..
303*
304* Test the input parameters.
305*
306 wantz = lsame( jobz, 'V' )
307 lower = lsame( uplo, 'L' )
308 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
309*
310 info = 0
311 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
312 info = -1
313 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
314 info = -2
315 ELSE IF( n.LT.0 ) THEN
316 info = -3
317 ELSE IF( lda.LT.max( 1, n ) ) THEN
318 info = -5
319 END IF
320*
321 IF( info.EQ.0 ) THEN
322 IF( n.LE.1 ) THEN
323 lwmin = 1
324 lrwmin = 1
325 liwmin = 1
326 ELSE
327 kd = ilaenv2stage( 1, 'CHETRD_2STAGE', jobz,
328 $ n, -1, -1, -1 )
329 ib = ilaenv2stage( 2, 'CHETRD_2STAGE', jobz,
330 $ n, kd, -1, -1 )
331 lhtrd = ilaenv2stage( 3, 'CHETRD_2STAGE', jobz,
332 $ n, kd, ib, -1 )
333 lwtrd = ilaenv2stage( 4, 'CHETRD_2STAGE', jobz,
334 $ n, kd, ib, -1 )
335 IF( wantz ) THEN
336 lwmin = 2*n + n*n
337 lrwmin = 1 + 5*n + 2*n**2
338 liwmin = 3 + 5*n
339 ELSE
340 lwmin = n + 1 + lhtrd + lwtrd
341 lrwmin = n
342 liwmin = 1
343 END IF
344 END IF
345 work( 1 ) = lwmin
346 rwork( 1 ) = lrwmin
347 iwork( 1 ) = liwmin
348*
349 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
350 info = -8
351 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
352 info = -10
353 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
354 info = -12
355 END IF
356 END IF
357*
358 IF( info.NE.0 ) THEN
359 CALL xerbla( 'CHEEVD_2STAGE', -info )
360 RETURN
361 ELSE IF( lquery ) THEN
362 RETURN
363 END IF
364*
365* Quick return if possible
366*
367 IF( n.EQ.0 )
368 $ RETURN
369*
370 IF( n.EQ.1 ) THEN
371 w( 1 ) = real( a( 1, 1 ) )
372 IF( wantz )
373 $ a( 1, 1 ) = cone
374 RETURN
375 END IF
376*
377* Get machine constants.
378*
379 safmin = slamch( 'Safe minimum' )
380 eps = slamch( 'Precision' )
381 smlnum = safmin / eps
382 bignum = one / smlnum
383 rmin = sqrt( smlnum )
384 rmax = sqrt( bignum )
385*
386* Scale matrix to allowable range, if necessary.
387*
388 anrm = clanhe( 'M', uplo, n, a, lda, rwork )
389 iscale = 0
390 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
391 iscale = 1
392 sigma = rmin / anrm
393 ELSE IF( anrm.GT.rmax ) THEN
394 iscale = 1
395 sigma = rmax / anrm
396 END IF
397 IF( iscale.EQ.1 )
398 $ CALL clascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
399*
400* Call CHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
401*
402 inde = 1
403 indrwk = inde + n
404 llrwk = lrwork - indrwk + 1
405 indtau = 1
406 indhous = indtau + n
407 indwrk = indhous + lhtrd
408 llwork = lwork - indwrk + 1
409 indwk2 = indwrk + n*n
410 llwrk2 = lwork - indwk2 + 1
411*
412 CALL chetrd_2stage( jobz, uplo, n, a, lda, w, rwork( inde ),
413 $ work( indtau ), work( indhous ), lhtrd,
414 $ work( indwrk ), llwork, iinfo )
415*
416* For eigenvalues only, call SSTERF. For eigenvectors, first call
417* CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
418* tridiagonal matrix, then call CUNMTR to multiply it to the
419* Householder transformations represented as Householder vectors in
420* A.
421*
422 IF( .NOT.wantz ) THEN
423 CALL ssterf( n, w, rwork( inde ), info )
424 ELSE
425 CALL cstedc( 'I', n, w, rwork( inde ), work( indwrk ), n,
426 $ work( indwk2 ), llwrk2, rwork( indrwk ), llrwk,
427 $ iwork, liwork, info )
428 CALL cunmtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
429 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
430 CALL clacpy( 'A', n, n, work( indwrk ), n, a, lda )
431 END IF
432*
433* If matrix was scaled, then rescale eigenvalues appropriately.
434*
435 IF( iscale.EQ.1 ) THEN
436 IF( info.EQ.0 ) THEN
437 imax = n
438 ELSE
439 imax = info - 1
440 END IF
441 CALL sscal( imax, one / sigma, w, 1 )
442 END IF
443*
444 work( 1 ) = lwmin
445 rwork( 1 ) = lrwmin
446 iwork( 1 ) = liwmin
447*
448 RETURN
449*
450* End of CHEEVD_2STAGE
451*

◆ cheevr()

subroutine cheevr ( character jobz,
character range,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
complex, dimension( ldz, * ) z,
integer ldz,
integer, dimension( * ) isuppz,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

CHEEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
!>
!> CHEEVR computes selected eigenvalues and, optionally, eigenvectors
!> of a complex Hermitian matrix A.  Eigenvalues and eigenvectors can
!> be selected by specifying either a range of values or a range of
!> indices for the desired eigenvalues.
!>
!> CHEEVR first reduces the matrix A to tridiagonal form T with a call
!> to CHETRD.  Then, whenever possible, CHEEVR calls CSTEMR to compute
!> the eigenspectrum using Relatively Robust Representations.  CSTEMR
!> 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 CSTEMR'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 : CHEEVR calls CSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> CHEEVR calls SSTEBZ and CSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of CSTEMR 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, SSTEBZ and
!>          CSTEIN 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 REAL
!>          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 REAL
!>          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 REAL
!>          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
!>          SLAMCH( '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 REAL array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is COMPLEX 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.
!>
[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 CSTEMR (tridiagonal
!>          matrix). The support of the eigenvectors of A is typically
!>          1:N because of the unitary transformations applied by CUNMTR.
!>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
!>
[out]WORK
!>          WORK is COMPLEX 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,2*N).
!>          For optimal efficiency, LWORK >= (NB+1)*N,
!>          where NB is the max of the blocksize for CHETRD and for
!>          CUNMTR as returned by ILAENV.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK, RWORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the optimal
!>          (and minimal) LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The length of the array RWORK.  LRWORK >= max(1,24*N).
!>
!>          If LRWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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
!>          (and minimal) LIWORK.
!>
[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 sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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:  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 354 of file cheevr.f.

357*
358* -- LAPACK driver routine --
359* -- LAPACK is a software package provided by Univ. of Tennessee, --
360* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
361*
362* .. Scalar Arguments ..
363 CHARACTER JOBZ, RANGE, UPLO
364 INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK,
365 $ M, N
366 REAL ABSTOL, VL, VU
367* ..
368* .. Array Arguments ..
369 INTEGER ISUPPZ( * ), IWORK( * )
370 REAL RWORK( * ), W( * )
371 COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * )
372* ..
373*
374* =====================================================================
375*
376* .. Parameters ..
377 REAL ZERO, ONE, TWO
378 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
379* ..
380* .. Local Scalars ..
381 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
382 $ WANTZ, TRYRAC
383 CHARACTER ORDER
384 INTEGER I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
385 $ INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
386 $ INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
387 $ LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
388 $ LWKOPT, LWMIN, NB, NSPLIT
389 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
390 $ SIGMA, SMLNUM, TMP1, VLL, VUU
391* ..
392* .. External Functions ..
393 LOGICAL LSAME
394 INTEGER ILAENV
395 REAL CLANSY, SLAMCH
396 EXTERNAL lsame, ilaenv, clansy, slamch
397* ..
398* .. External Subroutines ..
399 EXTERNAL chetrd, csscal, cstemr, cstein, cswap, cunmtr,
400 $ scopy, sscal, sstebz, ssterf, xerbla
401* ..
402* .. Intrinsic Functions ..
403 INTRINSIC max, min, real, sqrt
404* ..
405* .. Executable Statements ..
406*
407* Test the input parameters.
408*
409 ieeeok = ilaenv( 10, 'CHEEVR', 'N', 1, 2, 3, 4 )
410*
411 lower = lsame( uplo, 'L' )
412 wantz = lsame( jobz, 'V' )
413 alleig = lsame( range, 'A' )
414 valeig = lsame( range, 'V' )
415 indeig = lsame( range, 'I' )
416*
417 lquery = ( ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 ) .OR.
418 $ ( liwork.EQ.-1 ) )
419*
420 lrwmin = max( 1, 24*n )
421 liwmin = max( 1, 10*n )
422 lwmin = max( 1, 2*n )
423*
424 info = 0
425 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
426 info = -1
427 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
428 info = -2
429 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
430 info = -3
431 ELSE IF( n.LT.0 ) THEN
432 info = -4
433 ELSE IF( lda.LT.max( 1, n ) ) THEN
434 info = -6
435 ELSE
436 IF( valeig ) THEN
437 IF( n.GT.0 .AND. vu.LE.vl )
438 $ info = -8
439 ELSE IF( indeig ) THEN
440 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
441 info = -9
442 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
443 info = -10
444 END IF
445 END IF
446 END IF
447 IF( info.EQ.0 ) THEN
448 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
449 info = -15
450 END IF
451 END IF
452*
453 IF( info.EQ.0 ) THEN
454 nb = ilaenv( 1, 'CHETRD', uplo, n, -1, -1, -1 )
455 nb = max( nb, ilaenv( 1, 'CUNMTR', uplo, n, -1, -1, -1 ) )
456 lwkopt = max( ( nb+1 )*n, lwmin )
457 work( 1 ) = lwkopt
458 rwork( 1 ) = lrwmin
459 iwork( 1 ) = liwmin
460*
461 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
462 info = -18
463 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
464 info = -20
465 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
466 info = -22
467 END IF
468 END IF
469*
470 IF( info.NE.0 ) THEN
471 CALL xerbla( 'CHEEVR', -info )
472 RETURN
473 ELSE IF( lquery ) THEN
474 RETURN
475 END IF
476*
477* Quick return if possible
478*
479 m = 0
480 IF( n.EQ.0 ) THEN
481 work( 1 ) = 1
482 RETURN
483 END IF
484*
485 IF( n.EQ.1 ) THEN
486 work( 1 ) = 2
487 IF( alleig .OR. indeig ) THEN
488 m = 1
489 w( 1 ) = real( a( 1, 1 ) )
490 ELSE
491 IF( vl.LT.real( a( 1, 1 ) ) .AND. vu.GE.real( a( 1, 1 ) ) )
492 $ THEN
493 m = 1
494 w( 1 ) = real( a( 1, 1 ) )
495 END IF
496 END IF
497 IF( wantz ) THEN
498 z( 1, 1 ) = one
499 isuppz( 1 ) = 1
500 isuppz( 2 ) = 1
501 END IF
502 RETURN
503 END IF
504*
505* Get machine constants.
506*
507 safmin = slamch( 'Safe minimum' )
508 eps = slamch( 'Precision' )
509 smlnum = safmin / eps
510 bignum = one / smlnum
511 rmin = sqrt( smlnum )
512 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
513*
514* Scale matrix to allowable range, if necessary.
515*
516 iscale = 0
517 abstll = abstol
518 IF (valeig) THEN
519 vll = vl
520 vuu = vu
521 END IF
522 anrm = clansy( 'M', uplo, n, a, lda, rwork )
523 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
524 iscale = 1
525 sigma = rmin / anrm
526 ELSE IF( anrm.GT.rmax ) THEN
527 iscale = 1
528 sigma = rmax / anrm
529 END IF
530 IF( iscale.EQ.1 ) THEN
531 IF( lower ) THEN
532 DO 10 j = 1, n
533 CALL csscal( n-j+1, sigma, a( j, j ), 1 )
534 10 CONTINUE
535 ELSE
536 DO 20 j = 1, n
537 CALL csscal( j, sigma, a( 1, j ), 1 )
538 20 CONTINUE
539 END IF
540 IF( abstol.GT.0 )
541 $ abstll = abstol*sigma
542 IF( valeig ) THEN
543 vll = vl*sigma
544 vuu = vu*sigma
545 END IF
546 END IF
547
548* Initialize indices into workspaces. Note: The IWORK indices are
549* used only if SSTERF or CSTEMR fail.
550
551* WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
552* elementary reflectors used in CHETRD.
553 indtau = 1
554* INDWK is the starting offset of the remaining complex workspace,
555* and LLWORK is the remaining complex workspace size.
556 indwk = indtau + n
557 llwork = lwork - indwk + 1
558
559* RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal
560* entries.
561 indrd = 1
562* RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the
563* tridiagonal matrix from CHETRD.
564 indre = indrd + n
565* RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
566* -written by CSTEMR (the SSTERF path copies the diagonal to W).
567 indrdd = indre + n
568* RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over
569* -written while computing the eigenvalues in SSTERF and CSTEMR.
570 indree = indrdd + n
571* INDRWK is the starting offset of the left-over real workspace, and
572* LLRWORK is the remaining workspace size.
573 indrwk = indree + n
574 llrwork = lrwork - indrwk + 1
575
576* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and
577* stores the block indices of each of the M<=N eigenvalues.
578 indibl = 1
579* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and
580* stores the starting and finishing indices of each block.
581 indisp = indibl + n
582* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
583* that corresponding to eigenvectors that fail to converge in
584* SSTEIN. This information is discarded; if any fail, the driver
585* returns INFO > 0.
586 indifl = indisp + n
587* INDIWO is the offset of the remaining integer workspace.
588 indiwo = indifl + n
589
590*
591* Call CHETRD to reduce Hermitian matrix to tridiagonal form.
592*
593 CALL chetrd( uplo, n, a, lda, rwork( indrd ), rwork( indre ),
594 $ work( indtau ), work( indwk ), llwork, iinfo )
595*
596* If all eigenvalues are desired
597* then call SSTERF or CSTEMR and CUNMTR.
598*
599 test = .false.
600 IF( indeig ) THEN
601 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
602 test = .true.
603 END IF
604 END IF
605 IF( ( alleig.OR.test ) .AND. ( ieeeok.EQ.1 ) ) THEN
606 IF( .NOT.wantz ) THEN
607 CALL scopy( n, rwork( indrd ), 1, w, 1 )
608 CALL scopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
609 CALL ssterf( n, w, rwork( indree ), info )
610 ELSE
611 CALL scopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
612 CALL scopy( n, rwork( indrd ), 1, rwork( indrdd ), 1 )
613*
614 IF (abstol .LE. two*n*eps) THEN
615 tryrac = .true.
616 ELSE
617 tryrac = .false.
618 END IF
619 CALL cstemr( jobz, 'A', n, rwork( indrdd ),
620 $ rwork( indree ), vl, vu, il, iu, m, w,
621 $ z, ldz, n, isuppz, tryrac,
622 $ rwork( indrwk ), llrwork,
623 $ iwork, liwork, info )
624*
625* Apply unitary matrix used in reduction to tridiagonal
626* form to eigenvectors returned by CSTEMR.
627*
628 IF( wantz .AND. info.EQ.0 ) THEN
629 indwkn = indwk
630 llwrkn = lwork - indwkn + 1
631 CALL cunmtr( 'L', uplo, 'N', n, m, a, lda,
632 $ work( indtau ), z, ldz, work( indwkn ),
633 $ llwrkn, iinfo )
634 END IF
635 END IF
636*
637*
638 IF( info.EQ.0 ) THEN
639 m = n
640 GO TO 30
641 END IF
642 info = 0
643 END IF
644*
645* Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
646* Also call SSTEBZ and CSTEIN if CSTEMR fails.
647*
648 IF( wantz ) THEN
649 order = 'B'
650 ELSE
651 order = 'E'
652 END IF
653
654 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
655 $ rwork( indrd ), rwork( indre ), m, nsplit, w,
656 $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
657 $ iwork( indiwo ), info )
658*
659 IF( wantz ) THEN
660 CALL cstein( n, rwork( indrd ), rwork( indre ), m, w,
661 $ iwork( indibl ), iwork( indisp ), z, ldz,
662 $ rwork( indrwk ), iwork( indiwo ), iwork( indifl ),
663 $ info )
664*
665* Apply unitary matrix used in reduction to tridiagonal
666* form to eigenvectors returned by CSTEIN.
667*
668 indwkn = indwk
669 llwrkn = lwork - indwkn + 1
670 CALL cunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
671 $ ldz, work( indwkn ), llwrkn, iinfo )
672 END IF
673*
674* If matrix was scaled, then rescale eigenvalues appropriately.
675*
676 30 CONTINUE
677 IF( iscale.EQ.1 ) THEN
678 IF( info.EQ.0 ) THEN
679 imax = m
680 ELSE
681 imax = info - 1
682 END IF
683 CALL sscal( imax, one / sigma, w, 1 )
684 END IF
685*
686* If eigenvalues are not in order, then sort them, along with
687* eigenvectors.
688*
689 IF( wantz ) THEN
690 DO 50 j = 1, m - 1
691 i = 0
692 tmp1 = w( j )
693 DO 40 jj = j + 1, m
694 IF( w( jj ).LT.tmp1 ) THEN
695 i = jj
696 tmp1 = w( jj )
697 END IF
698 40 CONTINUE
699*
700 IF( i.NE.0 ) THEN
701 itmp1 = iwork( indibl+i-1 )
702 w( i ) = w( j )
703 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
704 w( j ) = tmp1
705 iwork( indibl+j-1 ) = itmp1
706 CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
707 END IF
708 50 CONTINUE
709 END IF
710*
711* Set WORK(1) to optimal workspace size.
712*
713 work( 1 ) = lwkopt
714 rwork( 1 ) = lrwmin
715 iwork( 1 ) = liwmin
716*
717 RETURN
718*
719* End of CHEEVR
720*
sstebz
subroutine sstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
SSTEBZ
Definition sstebz.f:273
cstein
subroutine cstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
CSTEIN
Definition cstein.f:182
cstemr
subroutine cstemr(jobz, range, n, d, e, vl, vu, il, iu, m, w, z, ldz, nzc, isuppz, tryrac, work, lwork, iwork, liwork, info)
CSTEMR
Definition cstemr.f:338
clansy
real function clansy(norm, uplo, n, a, lda, work)
CLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clansy.f:123
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
min
#define min(a, b)
Definition macros.h:20

◆ cheevr_2stage()

subroutine cheevr_2stage ( character jobz,
character range,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
complex, dimension( ldz, * ) z,
integer ldz,
integer, dimension( * ) isuppz,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

CHEEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
!>
!> CHEEVR_2STAGE computes selected eigenvalues and, optionally, eigenvectors
!> of a complex Hermitian 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.
!>
!> CHEEVR_2STAGE first reduces the matrix A to tridiagonal form T with a call
!> to CHETRD.  Then, whenever possible, CHEEVR_2STAGE calls CSTEMR to compute
!> eigenspectrum using Relatively Robust Representations.  CSTEMR
!> 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 CSTEMR'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 : CHEEVR_2STAGE calls CSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> CHEEVR_2STAGE calls SSTEBZ and CSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of CSTEMR 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, SSTEBZ and
!>          CSTEIN 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 REAL
!>          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 REAL
!>          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 REAL
!>          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
!>          SLAMCH( '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 REAL array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is COMPLEX 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.
!>
[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 CSTEMR (tridiagonal
!>          matrix). The support of the eigenvectors of A is typically 
!>          1:N because of the unitary transformations applied by CUNMTR.
!>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
!>
[out]WORK
!>          WORK is COMPLEX 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 + N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 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 sizes of the WORK, RWORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the optimal
!>          (and minimal) LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The length of the array RWORK.  LRWORK >= max(1,24*N).
!>
!>          If LRWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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
!>          (and minimal) LIWORK.
!>
[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 sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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:  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 402 of file cheevr_2stage.f.

406*
407 IMPLICIT NONE
408*
409* -- LAPACK driver routine --
410* -- LAPACK is a software package provided by Univ. of Tennessee, --
411* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
412*
413* .. Scalar Arguments ..
414 CHARACTER JOBZ, RANGE, UPLO
415 INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK,
416 $ M, N
417 REAL ABSTOL, VL, VU
418* ..
419* .. Array Arguments ..
420 INTEGER ISUPPZ( * ), IWORK( * )
421 REAL RWORK( * ), W( * )
422 COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * )
423* ..
424*
425* =====================================================================
426*
427* .. Parameters ..
428 REAL ZERO, ONE, TWO
429 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
430* ..
431* .. Local Scalars ..
432 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
433 $ WANTZ, TRYRAC
434 CHARACTER ORDER
435 INTEGER I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
436 $ INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
437 $ INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
438 $ LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
439 $ LWMIN, NSPLIT, LHTRD, LWTRD, KD, IB, INDHOUS
440 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
441 $ SIGMA, SMLNUM, TMP1, VLL, VUU
442* ..
443* .. External Functions ..
444 LOGICAL LSAME
445 INTEGER ILAENV, ILAENV2STAGE
446 REAL SLAMCH, CLANSY
447 EXTERNAL lsame, slamch, clansy, ilaenv, ilaenv2stage
448* ..
449* .. External Subroutines ..
450 EXTERNAL scopy, sscal, sstebz, ssterf, xerbla, csscal,
451 $ chetrd_2stage, cstemr, cstein, cswap, cunmtr
452* ..
453* .. Intrinsic Functions ..
454 INTRINSIC real, max, min, sqrt
455* ..
456* .. Executable Statements ..
457*
458* Test the input parameters.
459*
460 ieeeok = ilaenv( 10, 'CHEEVR', 'N', 1, 2, 3, 4 )
461*
462 lower = lsame( uplo, 'L' )
463 wantz = lsame( jobz, 'V' )
464 alleig = lsame( range, 'A' )
465 valeig = lsame( range, 'V' )
466 indeig = lsame( range, 'I' )
467*
468 lquery = ( ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 ) .OR.
469 $ ( liwork.EQ.-1 ) )
470*
471 kd = ilaenv2stage( 1, 'CHETRD_2STAGE', jobz, n, -1, -1, -1 )
472 ib = ilaenv2stage( 2, 'CHETRD_2STAGE', jobz, n, kd, -1, -1 )
473 lhtrd = ilaenv2stage( 3, 'CHETRD_2STAGE', jobz, n, kd, ib, -1 )
474 lwtrd = ilaenv2stage( 4, 'CHETRD_2STAGE', jobz, n, kd, ib, -1 )
475 lwmin = n + lhtrd + lwtrd
476 lrwmin = max( 1, 24*n )
477 liwmin = max( 1, 10*n )
478*
479 info = 0
480 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
481 info = -1
482 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
483 info = -2
484 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
485 info = -3
486 ELSE IF( n.LT.0 ) THEN
487 info = -4
488 ELSE IF( lda.LT.max( 1, n ) ) THEN
489 info = -6
490 ELSE
491 IF( valeig ) THEN
492 IF( n.GT.0 .AND. vu.LE.vl )
493 $ info = -8
494 ELSE IF( indeig ) THEN
495 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
496 info = -9
497 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
498 info = -10
499 END IF
500 END IF
501 END IF
502 IF( info.EQ.0 ) THEN
503 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
504 info = -15
505 END IF
506 END IF
507*
508 IF( info.EQ.0 ) THEN
509 work( 1 ) = lwmin
510 rwork( 1 ) = lrwmin
511 iwork( 1 ) = liwmin
512*
513 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
514 info = -18
515 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
516 info = -20
517 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
518 info = -22
519 END IF
520 END IF
521*
522 IF( info.NE.0 ) THEN
523 CALL xerbla( 'CHEEVR_2STAGE', -info )
524 RETURN
525 ELSE IF( lquery ) THEN
526 RETURN
527 END IF
528*
529* Quick return if possible
530*
531 m = 0
532 IF( n.EQ.0 ) THEN
533 work( 1 ) = 1
534 RETURN
535 END IF
536*
537 IF( n.EQ.1 ) THEN
538 work( 1 ) = 2
539 IF( alleig .OR. indeig ) THEN
540 m = 1
541 w( 1 ) = real( a( 1, 1 ) )
542 ELSE
543 IF( vl.LT.real( a( 1, 1 ) ) .AND. vu.GE.real( a( 1, 1 ) ) )
544 $ THEN
545 m = 1
546 w( 1 ) = real( a( 1, 1 ) )
547 END IF
548 END IF
549 IF( wantz ) THEN
550 z( 1, 1 ) = one
551 isuppz( 1 ) = 1
552 isuppz( 2 ) = 1
553 END IF
554 RETURN
555 END IF
556*
557* Get machine constants.
558*
559 safmin = slamch( 'Safe minimum' )
560 eps = slamch( 'Precision' )
561 smlnum = safmin / eps
562 bignum = one / smlnum
563 rmin = sqrt( smlnum )
564 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
565*
566* Scale matrix to allowable range, if necessary.
567*
568 iscale = 0
569 abstll = abstol
570 IF (valeig) THEN
571 vll = vl
572 vuu = vu
573 END IF
574 anrm = clansy( 'M', uplo, n, a, lda, rwork )
575 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
576 iscale = 1
577 sigma = rmin / anrm
578 ELSE IF( anrm.GT.rmax ) THEN
579 iscale = 1
580 sigma = rmax / anrm
581 END IF
582 IF( iscale.EQ.1 ) THEN
583 IF( lower ) THEN
584 DO 10 j = 1, n
585 CALL csscal( n-j+1, sigma, a( j, j ), 1 )
586 10 CONTINUE
587 ELSE
588 DO 20 j = 1, n
589 CALL csscal( j, sigma, a( 1, j ), 1 )
590 20 CONTINUE
591 END IF
592 IF( abstol.GT.0 )
593 $ abstll = abstol*sigma
594 IF( valeig ) THEN
595 vll = vl*sigma
596 vuu = vu*sigma
597 END IF
598 END IF
599
600* Initialize indices into workspaces. Note: The IWORK indices are
601* used only if SSTERF or CSTEMR fail.
602
603* WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
604* elementary reflectors used in CHETRD.
605 indtau = 1
606* INDWK is the starting offset of the remaining complex workspace,
607* and LLWORK is the remaining complex workspace size.
608 indhous = indtau + n
609 indwk = indhous + lhtrd
610 llwork = lwork - indwk + 1
611
612* RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal
613* entries.
614 indrd = 1
615* RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the
616* tridiagonal matrix from CHETRD.
617 indre = indrd + n
618* RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
619* -written by CSTEMR (the SSTERF path copies the diagonal to W).
620 indrdd = indre + n
621* RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over
622* -written while computing the eigenvalues in SSTERF and CSTEMR.
623 indree = indrdd + n
624* INDRWK is the starting offset of the left-over real workspace, and
625* LLRWORK is the remaining workspace size.
626 indrwk = indree + n
627 llrwork = lrwork - indrwk + 1
628
629* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and
630* stores the block indices of each of the M<=N eigenvalues.
631 indibl = 1
632* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and
633* stores the starting and finishing indices of each block.
634 indisp = indibl + n
635* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
636* that corresponding to eigenvectors that fail to converge in
637* CSTEIN. This information is discarded; if any fail, the driver
638* returns INFO > 0.
639 indifl = indisp + n
640* INDIWO is the offset of the remaining integer workspace.
641 indiwo = indifl + n
642
643*
644* Call CHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
645*
646 CALL chetrd_2stage( jobz, uplo, n, a, lda, rwork( indrd ),
647 $ rwork( indre ), work( indtau ),
648 $ work( indhous ), lhtrd,
649 $ work( indwk ), llwork, iinfo )
650*
651* If all eigenvalues are desired
652* then call SSTERF or CSTEMR and CUNMTR.
653*
654 test = .false.
655 IF( indeig ) THEN
656 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
657 test = .true.
658 END IF
659 END IF
660 IF( ( alleig.OR.test ) .AND. ( ieeeok.EQ.1 ) ) THEN
661 IF( .NOT.wantz ) THEN
662 CALL scopy( n, rwork( indrd ), 1, w, 1 )
663 CALL scopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
664 CALL ssterf( n, w, rwork( indree ), info )
665 ELSE
666 CALL scopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
667 CALL scopy( n, rwork( indrd ), 1, rwork( indrdd ), 1 )
668*
669 IF (abstol .LE. two*n*eps) THEN
670 tryrac = .true.
671 ELSE
672 tryrac = .false.
673 END IF
674 CALL cstemr( jobz, 'A', n, rwork( indrdd ),
675 $ rwork( indree ), vl, vu, il, iu, m, w,
676 $ z, ldz, n, isuppz, tryrac,
677 $ rwork( indrwk ), llrwork,
678 $ iwork, liwork, info )
679*
680* Apply unitary matrix used in reduction to tridiagonal
681* form to eigenvectors returned by CSTEMR.
682*
683 IF( wantz .AND. info.EQ.0 ) THEN
684 indwkn = indwk
685 llwrkn = lwork - indwkn + 1
686 CALL cunmtr( 'L', uplo, 'N', n, m, a, lda,
687 $ work( indtau ), z, ldz, work( indwkn ),
688 $ llwrkn, iinfo )
689 END IF
690 END IF
691*
692*
693 IF( info.EQ.0 ) THEN
694 m = n
695 GO TO 30
696 END IF
697 info = 0
698 END IF
699*
700* Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
701* Also call SSTEBZ and CSTEIN if CSTEMR fails.
702*
703 IF( wantz ) THEN
704 order = 'B'
705 ELSE
706 order = 'E'
707 END IF
708
709 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
710 $ rwork( indrd ), rwork( indre ), m, nsplit, w,
711 $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
712 $ iwork( indiwo ), info )
713*
714 IF( wantz ) THEN
715 CALL cstein( n, rwork( indrd ), rwork( indre ), m, w,
716 $ iwork( indibl ), iwork( indisp ), z, ldz,
717 $ rwork( indrwk ), iwork( indiwo ), iwork( indifl ),
718 $ info )
719*
720* Apply unitary matrix used in reduction to tridiagonal
721* form to eigenvectors returned by CSTEIN.
722*
723 indwkn = indwk
724 llwrkn = lwork - indwkn + 1
725 CALL cunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
726 $ ldz, work( indwkn ), llwrkn, iinfo )
727 END IF
728*
729* If matrix was scaled, then rescale eigenvalues appropriately.
730*
731 30 CONTINUE
732 IF( iscale.EQ.1 ) THEN
733 IF( info.EQ.0 ) THEN
734 imax = m
735 ELSE
736 imax = info - 1
737 END IF
738 CALL sscal( imax, one / sigma, w, 1 )
739 END IF
740*
741* If eigenvalues are not in order, then sort them, along with
742* eigenvectors.
743*
744 IF( wantz ) THEN
745 DO 50 j = 1, m - 1
746 i = 0
747 tmp1 = w( j )
748 DO 40 jj = j + 1, m
749 IF( w( jj ).LT.tmp1 ) THEN
750 i = jj
751 tmp1 = w( jj )
752 END IF
753 40 CONTINUE
754*
755 IF( i.NE.0 ) THEN
756 itmp1 = iwork( indibl+i-1 )
757 w( i ) = w( j )
758 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
759 w( j ) = tmp1
760 iwork( indibl+j-1 ) = itmp1
761 CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
762 END IF
763 50 CONTINUE
764 END IF
765*
766* Set WORK(1) to optimal workspace size.
767*
768 work( 1 ) = lwmin
769 rwork( 1 ) = lrwmin
770 iwork( 1 ) = liwmin
771*
772 RETURN
773*
774* End of CHEEVR_2STAGE
775*

◆ cheevx()

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

CHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
!>
!> CHEEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a complex Hermitian 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 REAL
!>          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 REAL
!>          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 REAL
!>          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*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('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 REAL array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 COMPLEX 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 2*N.
!>          For optimal efficiency, LWORK >= (NB+1)*N,
!>          where NB is the max of the blocksize for CHETRD and for
!>          CUNMTR as 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]RWORK
!>          RWORK is REAL array, dimension (7*N)
!>
[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 256 of file cheevx.f.

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

◆ cheevx_2stage()

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

CHEEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
!>
!> CHEEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
!> of a complex Hermitian 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 REAL
!>          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 REAL
!>          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 REAL
!>          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*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('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 REAL array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 COMPLEX 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 + N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 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]RWORK
!>          RWORK is REAL array, dimension (7*N)
!>
[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 303 of file cheevx_2stage.f.

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

◆ chegv()

subroutine chegv ( integer itype,
character jobz,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) w,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

CHEGV

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

Purpose:
!>
!> CHEGV computes all the eigenvalues, and optionally, the eigenvectors
!> of a complex generalized Hermitian-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 Hermitian 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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**H*B*Z = I;
!>          if ITYPE = 3, Z**H*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 COMPLEX array, dimension (LDB, N)
!>          On entry, the Hermitian 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**H*U or B = L*L**H.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is COMPLEX 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,2*N-1).
!>          For optimal efficiency, LWORK >= (NB+1)*N,
!>          where NB is the blocksize for CHETRD 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]RWORK
!>          RWORK is REAL array, dimension (max(1, 3*N-2))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  CPOTRF or CHEEV returned an error code:
!>             <= N:  if INFO = i, CHEEV 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 179 of file chegv.f.

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

◆ chegv_2stage()

subroutine chegv_2stage ( integer itype,
character jobz,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) w,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

CHEGV_2STAGE

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

Purpose:
!>
!> CHEGV_2STAGE computes all the eigenvalues, and optionally, the eigenvectors
!> of a complex generalized Hermitian-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 Hermitian 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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**H*B*Z = I;
!>          if ITYPE = 3, Z**H*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 COMPLEX array, dimension (LDB, N)
!>          On entry, the Hermitian 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**H*U or B = L*L**H.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is COMPLEX 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 + N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 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]RWORK
!>          RWORK is REAL array, dimension (max(1, 3*N-2))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  CPOTRF or CHEEV returned an error code:
!>             <= N:  if INFO = i, CHEEV 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 230 of file chegv_2stage.f.

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

◆ chegvd()

subroutine chegvd ( integer itype,
character jobz,
character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) w,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

CHEGVD

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

Purpose:
!>
!> CHEGVD computes all the eigenvalues, and optionally, the eigenvectors
!> of a complex generalized Hermitian-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 Hermitian 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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**H*B*Z = I;
!>          if ITYPE = 3, Z**H*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 COMPLEX array, dimension (LDB, N)
!>          On entry, the Hermitian 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**H*U or B = L*L**H.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is COMPLEX 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.
!>          If N <= 1,                LWORK >= 1.
!>          If JOBZ  = 'N' and N > 1, LWORK >= N + 1.
!>          If JOBZ  = 'V' and N > 1, LWORK >= 2*N + N**2.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK, RWORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of the array RWORK.
!>          If N <= 1,                LRWORK >= 1.
!>          If JOBZ  = 'N' and N > 1, LRWORK >= N.
!>          If JOBZ  = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
!>
!>          If LRWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK 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:  CPOTRF or CHEEVD 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 CHEEVD 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 247 of file chegvd.f.

249*
250* -- LAPACK driver routine --
251* -- LAPACK is a software package provided by Univ. of Tennessee, --
252* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253*
254* .. Scalar Arguments ..
255 CHARACTER JOBZ, UPLO
256 INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N
257* ..
258* .. Array Arguments ..
259 INTEGER IWORK( * )
260 REAL RWORK( * ), W( * )
261 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
262* ..
263*
264* =====================================================================
265*
266* .. Parameters ..
267 COMPLEX CONE
268 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
269* ..
270* .. Local Scalars ..
271 LOGICAL LQUERY, UPPER, WANTZ
272 CHARACTER TRANS
273 INTEGER LIOPT, LIWMIN, LOPT, LROPT, LRWMIN, LWMIN
274* ..
275* .. External Functions ..
276 LOGICAL LSAME
277 EXTERNAL lsame
278* ..
279* .. External Subroutines ..
280 EXTERNAL cheevd, chegst, cpotrf, ctrmm, ctrsm, xerbla
281* ..
282* .. Intrinsic Functions ..
283 INTRINSIC max, real
284* ..
285* .. Executable Statements ..
286*
287* Test the input parameters.
288*
289 wantz = lsame( jobz, 'V' )
290 upper = lsame( uplo, 'U' )
291 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
292*
293 info = 0
294 IF( n.LE.1 ) THEN
295 lwmin = 1
296 lrwmin = 1
297 liwmin = 1
298 ELSE IF( wantz ) THEN
299 lwmin = 2*n + n*n
300 lrwmin = 1 + 5*n + 2*n*n
301 liwmin = 3 + 5*n
302 ELSE
303 lwmin = n + 1
304 lrwmin = n
305 liwmin = 1
306 END IF
307 lopt = lwmin
308 lropt = lrwmin
309 liopt = liwmin
310 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
311 info = -1
312 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
313 info = -2
314 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) 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 IF( ldb.LT.max( 1, n ) ) THEN
321 info = -8
322 END IF
323*
324 IF( info.EQ.0 ) THEN
325 work( 1 ) = lopt
326 rwork( 1 ) = lropt
327 iwork( 1 ) = liopt
328*
329 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
330 info = -11
331 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
332 info = -13
333 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
334 info = -15
335 END IF
336 END IF
337*
338 IF( info.NE.0 ) THEN
339 CALL xerbla( 'CHEGVD', -info )
340 RETURN
341 ELSE IF( lquery ) THEN
342 RETURN
343 END IF
344*
345* Quick return if possible
346*
347 IF( n.EQ.0 )
348 $ RETURN
349*
350* Form a Cholesky factorization of B.
351*
352 CALL cpotrf( uplo, n, b, ldb, info )
353 IF( info.NE.0 ) THEN
354 info = n + info
355 RETURN
356 END IF
357*
358* Transform problem to standard eigenvalue problem and solve.
359*
360 CALL chegst( itype, uplo, n, a, lda, b, ldb, info )
361 CALL cheevd( jobz, uplo, n, a, lda, w, work, lwork, rwork, lrwork,
362 $ iwork, liwork, info )
363 lopt = max( real( lopt ), real( work( 1 ) ) )
364 lropt = max( real( lropt ), real( rwork( 1 ) ) )
365 liopt = max( real( liopt ), real( iwork( 1 ) ) )
366*
367 IF( wantz .AND. info.EQ.0 ) THEN
368*
369* Backtransform eigenvectors to the original problem.
370*
371 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
372*
373* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
374* backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
375*
376 IF( upper ) THEN
377 trans = 'N'
378 ELSE
379 trans = 'C'
380 END IF
381*
382 CALL ctrsm( 'Left', uplo, trans, 'Non-unit', n, n, cone,
383 $ b, ldb, a, lda )
384*
385 ELSE IF( itype.EQ.3 ) THEN
386*
387* For B*A*x=(lambda)*x;
388* backtransform eigenvectors: x = L*y or U**H *y
389*
390 IF( upper ) THEN
391 trans = 'C'
392 ELSE
393 trans = 'N'
394 END IF
395*
396 CALL ctrmm( 'Left', uplo, trans, 'Non-unit', n, n, cone,
397 $ b, ldb, a, lda )
398 END IF
399 END IF
400*
401 work( 1 ) = lopt
402 rwork( 1 ) = lropt
403 iwork( 1 ) = liopt
404*
405 RETURN
406*
407* End of CHEGVD
408*
cheevd
subroutine cheevd(jobz, uplo, n, a, lda, w, work, lwork, rwork, lrwork, iwork, liwork, info)
CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
Definition cheevd.f:205

◆ chegvx()

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

CHEGVX

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

Purpose:
!>
!> CHEGVX computes selected eigenvalues, and optionally, eigenvectors
!> of a complex generalized Hermitian-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 Hermitian 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 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 COMPLEX array, dimension (LDA, N)
!>          On entry, the Hermitian 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 COMPLEX array, dimension (LDB, N)
!>          On entry, the Hermitian 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**H*U or B = L*L**H.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[in]VL
!>          VL is REAL
!>
!>          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 REAL
!>
!>          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 REAL
!>          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*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('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 REAL array, dimension (N)
!>          The first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 COMPLEX 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,2*N).
!>          For optimal efficiency, LWORK >= (NB+1)*N,
!>          where NB is the blocksize for CHETRD 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]RWORK
!>          RWORK is REAL array, dimension (7*N)
!>
[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:  CPOTRF or CHEEVX returned an error code:
!>             <= N:  if INFO = i, CHEEVX 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 304 of file chegvx.f.

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