OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
complex16 Other Eigenvalue routines
LAPACK » Eigenvalue

Functions

subroutine zggglm (n, m, p, a, lda, b, ldb, d, x, y, work, lwork, info)
 ZGGGLM
subroutine zhbev (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, rwork, info)
  ZHBEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhbev_2stage (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, rwork, info)
  ZHBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhbevd (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
  ZHBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhbevd_2stage (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
  ZHBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhbevx (jobz, range, uplo, n, kd, ab, ldab, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
  ZHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhbevx_2stage (jobz, range, uplo, n, kd, ab, ldab, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, rwork, iwork, ifail, info)
  ZHBEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhbgv (jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z, ldz, work, rwork, info)
 ZHBGV
subroutine zhbgvd (jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
 ZHBGVD
subroutine zhbgvx (jobz, range, uplo, n, ka, kb, ab, ldab, bb, ldbb, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
 ZHBGVX
subroutine zhpev (jobz, uplo, n, ap, w, z, ldz, work, rwork, info)
  ZHPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhpevd (jobz, uplo, n, ap, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
  ZHPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhpevx (jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
  ZHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine zhpgv (itype, jobz, uplo, n, ap, bp, w, z, ldz, work, rwork, info)
 ZHPGV
subroutine zhpgvd (itype, jobz, uplo, n, ap, bp, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
 ZHPGVD
subroutine zhpgvx (itype, jobz, range, uplo, n, ap, bp, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
 ZHPGVX

Detailed Description

This is the group of complex16 Other Eigenvalue routines

Function Documentation

◆ zggglm()

subroutine zggglm ( integer n,
integer m,
integer p,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( * ) d,
complex*16, dimension( * ) x,
complex*16, dimension( * ) y,
complex*16, dimension( * ) work,
integer lwork,
integer info )

ZGGGLM

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

Purpose:
!>
!> ZGGGLM solves a general Gauss-Markov linear model (GLM) problem:
!>
!>         minimize || y ||_2   subject to   d = A*x + B*y
!>             x
!>
!> where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
!> given N-vector. It is assumed that M <= N <= M+P, and
!>
!>            rank(A) = M    and    rank( A B ) = N.
!>
!> Under these assumptions, the constrained equation is always
!> consistent, and there is a unique solution x and a minimal 2-norm
!> solution y, which is obtained using a generalized QR factorization
!> of the matrices (A, B) given by
!>
!>    A = Q*(R),   B = Q*T*Z.
!>          (0)
!>
!> In particular, if matrix B is square nonsingular, then the problem
!> GLM is equivalent to the following weighted linear least squares
!> problem
!>
!>              minimize || inv(B)*(d-A*x) ||_2
!>                  x
!>
!> where inv(B) denotes the inverse of B.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The number of rows of the matrices A and B.  N >= 0.
!>
[in]M
!>          M is INTEGER
!>          The number of columns of the matrix A.  0 <= M <= N.
!>
[in]P
!>          P is INTEGER
!>          The number of columns of the matrix B.  P >= N-M.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA,M)
!>          On entry, the N-by-M matrix A.
!>          On exit, the upper triangular part of the array A contains
!>          the M-by-M upper triangular matrix R.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,N).
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,P)
!>          On entry, the N-by-P matrix B.
!>          On exit, if N <= P, the upper triangle of the subarray
!>          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
!>          if N > P, the elements on and above the (N-P)th subdiagonal
!>          contain the N-by-P upper trapezoidal matrix T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,N).
!>
[in,out]D
!>          D is COMPLEX*16 array, dimension (N)
!>          On entry, D is the left hand side of the GLM equation.
!>          On exit, D is destroyed.
!>
[out]X
!>          X is COMPLEX*16 array, dimension (M)
!>
[out]Y
!>          Y is COMPLEX*16 array, dimension (P)
!>
!>          On exit, X and Y are the solutions of the GLM problem.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= max(1,N+M+P).
!>          For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB,
!>          where NB is an upper bound for the optimal blocksizes for
!>          ZGEQRF, ZGERQF, ZUNMQR and ZUNMRQ.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1:  the upper triangular factor R associated with A in the
!>                generalized QR factorization of the pair (A, B) is
!>                singular, so that rank(A) < M; the least squares
!>                solution could not be computed.
!>          = 2:  the bottom (N-M) by (N-M) part of the upper trapezoidal
!>                factor T associated with B in the generalized QR
!>                factorization of the pair (A, B) is singular, so that
!>                rank( A B ) < N; the least squares solution could not
!>                be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 183 of file zggglm.f.

185*
186* -- LAPACK driver routine --
187* -- LAPACK is a software package provided by Univ. of Tennessee, --
188* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190* .. Scalar Arguments ..
191 INTEGER INFO, LDA, LDB, LWORK, M, N, P
192* ..
193* .. Array Arguments ..
194 COMPLEX*16 A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195 $ X( * ), Y( * )
196* ..
197*
198* ===================================================================
199*
200* .. Parameters ..
201 COMPLEX*16 CZERO, CONE
202 parameter( czero = ( 0.0d+0, 0.0d+0 ),
203 $ cone = ( 1.0d+0, 0.0d+0 ) )
204* ..
205* .. Local Scalars ..
206 LOGICAL LQUERY
207 INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
208 $ NB4, NP
209* ..
210* .. External Subroutines ..
211 EXTERNAL xerbla, zcopy, zgemv, zggqrf, ztrtrs, zunmqr,
212 $ zunmrq
213* ..
214* .. External Functions ..
215 INTEGER ILAENV
216 EXTERNAL ilaenv
217* ..
218* .. Intrinsic Functions ..
219 INTRINSIC int, max, min
220* ..
221* .. Executable Statements ..
222*
223* Test the input parameters
224*
225 info = 0
226 np = min( n, p )
227 lquery = ( lwork.EQ.-1 )
228 IF( n.LT.0 ) THEN
229 info = -1
230 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
231 info = -2
232 ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
233 info = -3
234 ELSE IF( lda.LT.max( 1, n ) ) THEN
235 info = -5
236 ELSE IF( ldb.LT.max( 1, n ) ) THEN
237 info = -7
238 END IF
239*
240* Calculate workspace
241*
242 IF( info.EQ.0) THEN
243 IF( n.EQ.0 ) THEN
244 lwkmin = 1
245 lwkopt = 1
246 ELSE
247 nb1 = ilaenv( 1, 'ZGEQRF', ' ', n, m, -1, -1 )
248 nb2 = ilaenv( 1, 'ZGERQF', ' ', n, m, -1, -1 )
249 nb3 = ilaenv( 1, 'ZUNMQR', ' ', n, m, p, -1 )
250 nb4 = ilaenv( 1, 'ZUNMRQ', ' ', n, m, p, -1 )
251 nb = max( nb1, nb2, nb3, nb4 )
252 lwkmin = m + n + p
253 lwkopt = m + np + max( n, p )*nb
254 END IF
255 work( 1 ) = lwkopt
256*
257 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
258 info = -12
259 END IF
260 END IF
261*
262 IF( info.NE.0 ) THEN
263 CALL xerbla( 'ZGGGLM', -info )
264 RETURN
265 ELSE IF( lquery ) THEN
266 RETURN
267 END IF
268*
269* Quick return if possible
270*
271 IF( n.EQ.0 ) THEN
272 DO i = 1, m
273 x(i) = czero
274 END DO
275 DO i = 1, p
276 y(i) = czero
277 END DO
278 RETURN
279 END IF
280*
281* Compute the GQR factorization of matrices A and B:
282*
283* Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M
284* ( 0 ) N-M ( 0 T22 ) N-M
285* M M+P-N N-M
286*
287* where R11 and T22 are upper triangular, and Q and Z are
288* unitary.
289*
290 CALL zggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
291 $ work( m+np+1 ), lwork-m-np, info )
292 lopt = dble( work( m+np+1 ) )
293*
294* Update left-hand-side vector d = Q**H*d = ( d1 ) M
295* ( d2 ) N-M
296*
297 CALL zunmqr( 'Left', 'Conjugate transpose', n, 1, m, a, lda, work,
298 $ d, max( 1, n ), work( m+np+1 ), lwork-m-np, info )
299 lopt = max( lopt, int( work( m+np+1 ) ) )
300*
301* Solve T22*y2 = d2 for y2
302*
303 IF( n.GT.m ) THEN
304 CALL ztrtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
305 $ b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
306*
307 IF( info.GT.0 ) THEN
308 info = 1
309 RETURN
310 END IF
311*
312 CALL zcopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
313 END IF
314*
315* Set y1 = 0
316*
317 DO 10 i = 1, m + p - n
318 y( i ) = czero
319 10 CONTINUE
320*
321* Update d1 = d1 - T12*y2
322*
323 CALL zgemv( 'No transpose', m, n-m, -cone, b( 1, m+p-n+1 ), ldb,
324 $ y( m+p-n+1 ), 1, cone, d, 1 )
325*
326* Solve triangular system: R11*x = d1
327*
328 IF( m.GT.0 ) THEN
329 CALL ztrtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a, lda,
330 $ d, m, info )
331*
332 IF( info.GT.0 ) THEN
333 info = 2
334 RETURN
335 END IF
336*
337* Copy D to X
338*
339 CALL zcopy( m, d, 1, x, 1 )
340 END IF
341*
342* Backward transformation y = Z**H *y
343*
344 CALL zunmrq( 'Left', 'Conjugate transpose', p, 1, np,
345 $ b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
346 $ max( 1, p ), work( m+np+1 ), lwork-m-np, info )
347 work( 1 ) = m + np + max( lopt, int( work( m+np+1 ) ) )
348*
349 RETURN
350*
351* End of ZGGGLM
352*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
zunmqr
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:167
ztrtrs
subroutine ztrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
ZTRTRS
Definition ztrtrs.f:140
zggqrf
subroutine zggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
ZGGQRF
Definition zggqrf.f:215
zunmrq
subroutine zunmrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMRQ
Definition zunmrq.f:167
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:158
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ zhbev()

subroutine zhbev ( character jobz,
character uplo,
integer n,
integer kd,
complex*16, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZHBEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHBEV computes all the eigenvalues and, optionally, eigenvectors of
!> a complex Hermitian band 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]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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 150 of file zhbev.f.

152*
153* -- LAPACK driver routine --
154* -- LAPACK is a software package provided by Univ. of Tennessee, --
155* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
156*
157* .. Scalar Arguments ..
158 CHARACTER JOBZ, UPLO
159 INTEGER INFO, KD, LDAB, LDZ, N
160* ..
161* .. Array Arguments ..
162 DOUBLE PRECISION RWORK( * ), W( * )
163 COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
164* ..
165*
166* =====================================================================
167*
168* .. Parameters ..
169 DOUBLE PRECISION ZERO, ONE
170 parameter( zero = 0.0d0, one = 1.0d0 )
171* ..
172* .. Local Scalars ..
173 LOGICAL LOWER, WANTZ
174 INTEGER IINFO, IMAX, INDE, INDRWK, ISCALE
175 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
176 $ SMLNUM
177* ..
178* .. External Functions ..
179 LOGICAL LSAME
180 DOUBLE PRECISION DLAMCH, ZLANHB
181 EXTERNAL lsame, dlamch, zlanhb
182* ..
183* .. External Subroutines ..
184 EXTERNAL dscal, dsterf, xerbla, zhbtrd, zlascl, zsteqr
185* ..
186* .. Intrinsic Functions ..
187 INTRINSIC sqrt
188* ..
189* .. Executable Statements ..
190*
191* Test the input parameters.
192*
193 wantz = lsame( jobz, 'V' )
194 lower = lsame( uplo, 'L' )
195*
196 info = 0
197 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
198 info = -1
199 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
200 info = -2
201 ELSE IF( n.LT.0 ) THEN
202 info = -3
203 ELSE IF( kd.LT.0 ) THEN
204 info = -4
205 ELSE IF( ldab.LT.kd+1 ) THEN
206 info = -6
207 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
208 info = -9
209 END IF
210*
211 IF( info.NE.0 ) THEN
212 CALL xerbla( 'ZHBEV ', -info )
213 RETURN
214 END IF
215*
216* Quick return if possible
217*
218 IF( n.EQ.0 )
219 $ RETURN
220*
221 IF( n.EQ.1 ) THEN
222 IF( lower ) THEN
223 w( 1 ) = dble( ab( 1, 1 ) )
224 ELSE
225 w( 1 ) = dble( ab( kd+1, 1 ) )
226 END IF
227 IF( wantz )
228 $ z( 1, 1 ) = one
229 RETURN
230 END IF
231*
232* Get machine constants.
233*
234 safmin = dlamch( 'Safe minimum' )
235 eps = dlamch( 'Precision' )
236 smlnum = safmin / eps
237 bignum = one / smlnum
238 rmin = sqrt( smlnum )
239 rmax = sqrt( bignum )
240*
241* Scale matrix to allowable range, if necessary.
242*
243 anrm = zlanhb( 'M', uplo, n, kd, ab, ldab, rwork )
244 iscale = 0
245 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
246 iscale = 1
247 sigma = rmin / anrm
248 ELSE IF( anrm.GT.rmax ) THEN
249 iscale = 1
250 sigma = rmax / anrm
251 END IF
252 IF( iscale.EQ.1 ) THEN
253 IF( lower ) THEN
254 CALL zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
255 ELSE
256 CALL zlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
257 END IF
258 END IF
259*
260* Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
261*
262 inde = 1
263 CALL zhbtrd( jobz, uplo, n, kd, ab, ldab, w, rwork( inde ), z,
264 $ ldz, work, iinfo )
265*
266* For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEQR.
267*
268 IF( .NOT.wantz ) THEN
269 CALL dsterf( n, w, rwork( inde ), info )
270 ELSE
271 indrwk = inde + n
272 CALL zsteqr( jobz, n, w, rwork( inde ), z, ldz,
273 $ rwork( indrwk ), 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 dscal( imax, one / sigma, w, 1 )
285 END IF
286*
287 RETURN
288*
289* End of ZHBEV
290*
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
zlanhb
double precision function zlanhb(norm, uplo, n, k, ab, ldab, work)
ZLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlanhb.f:132
zlascl
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:143
zsteqr
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
zhbtrd
subroutine zhbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
ZHBTRD
Definition zhbtrd.f:163
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69

◆ zhbev_2stage()

subroutine zhbev_2stage ( character jobz,
character uplo,
integer n,
integer kd,
complex*16, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer info )

ZHBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHBEV_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
!> a complex Hermitian band 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]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension LWORK
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK. LWORK >= 1, when N <= 1;
!>          otherwise  
!>          If JOBZ = 'N' and N > 1, LWORK must be queried.
!>                                   LWORK = MAX(1, dimension) where
!>                                   dimension = (2KD+1)*N + KD*NTHREADS
!>                                   where KD is the size of the band.
!>                                   NTHREADS is the number of threads used when
!>                                   openMP compilation is enabled, otherwise =1.
!>          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK, 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 DOUBLE PRECISION 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 209 of file zhbev_2stage.f.

211*
212 IMPLICIT NONE
213*
214* -- LAPACK driver routine --
215* -- LAPACK is a software package provided by Univ. of Tennessee, --
216* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217*
218* .. Scalar Arguments ..
219 CHARACTER JOBZ, UPLO
220 INTEGER INFO, KD, LDAB, LDZ, N, LWORK
221* ..
222* .. Array Arguments ..
223 DOUBLE PRECISION RWORK( * ), W( * )
224 COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
225* ..
226*
227* =====================================================================
228*
229* .. Parameters ..
230 DOUBLE PRECISION ZERO, ONE
231 parameter( zero = 0.0d0, one = 1.0d0 )
232* ..
233* .. Local Scalars ..
234 LOGICAL LOWER, WANTZ, LQUERY
235 INTEGER IINFO, IMAX, INDE, INDWRK, INDRWK, ISCALE,
236 $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS
237 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
238 $ SMLNUM
239* ..
240* .. External Functions ..
241 LOGICAL LSAME
242 INTEGER ILAENV2STAGE
243 DOUBLE PRECISION DLAMCH, ZLANHB
244 EXTERNAL lsame, dlamch, zlanhb, ilaenv2stage
245* ..
246* .. External Subroutines ..
247 EXTERNAL dscal, dsterf, xerbla, zlascl, zsteqr,
248 $ zhetrd_2stage, zhetrd_hb2st
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC dble, sqrt
252* ..
253* .. Executable Statements ..
254*
255* Test the input parameters.
256*
257 wantz = lsame( jobz, 'V' )
258 lower = lsame( uplo, 'L' )
259 lquery = ( lwork.EQ.-1 )
260*
261 info = 0
262 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
263 info = -1
264 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
265 info = -2
266 ELSE IF( n.LT.0 ) THEN
267 info = -3
268 ELSE IF( kd.LT.0 ) THEN
269 info = -4
270 ELSE IF( ldab.LT.kd+1 ) THEN
271 info = -6
272 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
273 info = -9
274 END IF
275*
276 IF( info.EQ.0 ) THEN
277 IF( n.LE.1 ) THEN
278 lwmin = 1
279 work( 1 ) = lwmin
280 ELSE
281 ib = ilaenv2stage( 2, 'ZHETRD_HB2ST', jobz,
282 $ n, kd, -1, -1 )
283 lhtrd = ilaenv2stage( 3, 'ZHETRD_HB2ST', jobz,
284 $ n, kd, ib, -1 )
285 lwtrd = ilaenv2stage( 4, 'ZHETRD_HB2ST', jobz,
286 $ n, kd, ib, -1 )
287 lwmin = lhtrd + lwtrd
288 work( 1 ) = lwmin
289 ENDIF
290*
291 IF( lwork.LT.lwmin .AND. .NOT.lquery )
292 $ info = -11
293 END IF
294*
295 IF( info.NE.0 ) THEN
296 CALL xerbla( 'ZHBEV_2STAGE ', -info )
297 RETURN
298 ELSE IF( lquery ) THEN
299 RETURN
300 END IF
301*
302* Quick return if possible
303*
304 IF( n.EQ.0 )
305 $ RETURN
306*
307 IF( n.EQ.1 ) THEN
308 IF( lower ) THEN
309 w( 1 ) = dble( ab( 1, 1 ) )
310 ELSE
311 w( 1 ) = dble( ab( kd+1, 1 ) )
312 END IF
313 IF( wantz )
314 $ z( 1, 1 ) = one
315 RETURN
316 END IF
317*
318* Get machine constants.
319*
320 safmin = dlamch( 'Safe minimum' )
321 eps = dlamch( 'Precision' )
322 smlnum = safmin / eps
323 bignum = one / smlnum
324 rmin = sqrt( smlnum )
325 rmax = sqrt( bignum )
326*
327* Scale matrix to allowable range, if necessary.
328*
329 anrm = zlanhb( 'M', uplo, n, kd, ab, ldab, rwork )
330 iscale = 0
331 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
332 iscale = 1
333 sigma = rmin / anrm
334 ELSE IF( anrm.GT.rmax ) THEN
335 iscale = 1
336 sigma = rmax / anrm
337 END IF
338 IF( iscale.EQ.1 ) THEN
339 IF( lower ) THEN
340 CALL zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
341 ELSE
342 CALL zlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
343 END IF
344 END IF
345*
346* Call ZHBTRD_HB2ST to reduce Hermitian band matrix to tridiagonal form.
347*
348 inde = 1
349 indhous = 1
350 indwrk = indhous + lhtrd
351 llwork = lwork - indwrk + 1
352*
353 CALL zhetrd_hb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
354 $ rwork( inde ), work( indhous ), lhtrd,
355 $ work( indwrk ), llwork, iinfo )
356*
357* For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEQR.
358*
359 IF( .NOT.wantz ) THEN
360 CALL dsterf( n, w, rwork( inde ), info )
361 ELSE
362 indrwk = inde + n
363 CALL zsteqr( jobz, n, w, rwork( inde ), z, ldz,
364 $ rwork( indrwk ), info )
365 END IF
366*
367* If matrix was scaled, then rescale eigenvalues appropriately.
368*
369 IF( iscale.EQ.1 ) THEN
370 IF( info.EQ.0 ) THEN
371 imax = n
372 ELSE
373 imax = info - 1
374 END IF
375 CALL dscal( imax, one / sigma, w, 1 )
376 END IF
377*
378* Set WORK(1) to optimal workspace size.
379*
380 work( 1 ) = lwmin
381*
382 RETURN
383*
384* End of ZHBEV_2STAGE
385*
ilaenv2stage
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
Definition ilaenv2stage.f:149
zhetrd_2stage
subroutine zhetrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
ZHETRD_2STAGE
Definition zhetrd_2stage.f:224
zhetrd_hb2st
subroutine zhetrd_hb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
ZHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T
Definition zhetrd_hb2st.F:230

◆ zhbevd()

subroutine zhbevd ( character jobz,
character uplo,
integer n,
integer kd,
complex*16, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

ZHBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHBEVD computes all the eigenvalues and, optionally, eigenvectors of
!> a complex Hermitian band 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]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 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 at least N.
!>          If JOBZ = 'V' and N > 1, LWORK must be at least 2*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 DOUBLE PRECISION array,
!>                                         dimension (LRWORK)
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of 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 array IWORK.
!>          If JOBZ = 'N' or 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, 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 213 of file zhbevd.f.

215*
216* -- LAPACK driver routine --
217* -- LAPACK is a software package provided by Univ. of Tennessee, --
218* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
219*
220* .. Scalar Arguments ..
221 CHARACTER JOBZ, UPLO
222 INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
223* ..
224* .. Array Arguments ..
225 INTEGER IWORK( * )
226 DOUBLE PRECISION RWORK( * ), W( * )
227 COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
228* ..
229*
230* =====================================================================
231*
232* .. Parameters ..
233 DOUBLE PRECISION ZERO, ONE
234 parameter( zero = 0.0d0, one = 1.0d0 )
235 COMPLEX*16 CZERO, CONE
236 parameter( czero = ( 0.0d0, 0.0d0 ),
237 $ cone = ( 1.0d0, 0.0d0 ) )
238* ..
239* .. Local Scalars ..
240 LOGICAL LOWER, LQUERY, WANTZ
241 INTEGER IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE,
242 $ LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN
243 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
244 $ SMLNUM
245* ..
246* .. External Functions ..
247 LOGICAL LSAME
248 DOUBLE PRECISION DLAMCH, ZLANHB
249 EXTERNAL lsame, dlamch, zlanhb
250* ..
251* .. External Subroutines ..
252 EXTERNAL dscal, dsterf, xerbla, zgemm, zhbtrd, zlacpy,
253 $ zlascl, zstedc
254* ..
255* .. Intrinsic Functions ..
256 INTRINSIC sqrt
257* ..
258* .. Executable Statements ..
259*
260* Test the input parameters.
261*
262 wantz = lsame( jobz, 'V' )
263 lower = lsame( uplo, 'L' )
264 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 .OR. lrwork.EQ.-1 )
265*
266 info = 0
267 IF( n.LE.1 ) THEN
268 lwmin = 1
269 lrwmin = 1
270 liwmin = 1
271 ELSE
272 IF( wantz ) THEN
273 lwmin = 2*n**2
274 lrwmin = 1 + 5*n + 2*n**2
275 liwmin = 3 + 5*n
276 ELSE
277 lwmin = n
278 lrwmin = n
279 liwmin = 1
280 END IF
281 END IF
282 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
283 info = -1
284 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
285 info = -2
286 ELSE IF( n.LT.0 ) THEN
287 info = -3
288 ELSE IF( kd.LT.0 ) THEN
289 info = -4
290 ELSE IF( ldab.LT.kd+1 ) THEN
291 info = -6
292 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
293 info = -9
294 END IF
295*
296 IF( info.EQ.0 ) THEN
297 work( 1 ) = lwmin
298 rwork( 1 ) = lrwmin
299 iwork( 1 ) = liwmin
300*
301 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
302 info = -11
303 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
304 info = -13
305 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
306 info = -15
307 END IF
308 END IF
309*
310 IF( info.NE.0 ) THEN
311 CALL xerbla( 'ZHBEVD', -info )
312 RETURN
313 ELSE IF( lquery ) THEN
314 RETURN
315 END IF
316*
317* Quick return if possible
318*
319 IF( n.EQ.0 )
320 $ RETURN
321*
322 IF( n.EQ.1 ) THEN
323 w( 1 ) = dble( ab( 1, 1 ) )
324 IF( wantz )
325 $ z( 1, 1 ) = cone
326 RETURN
327 END IF
328*
329* Get machine constants.
330*
331 safmin = dlamch( 'Safe minimum' )
332 eps = dlamch( 'Precision' )
333 smlnum = safmin / eps
334 bignum = one / smlnum
335 rmin = sqrt( smlnum )
336 rmax = sqrt( bignum )
337*
338* Scale matrix to allowable range, if necessary.
339*
340 anrm = zlanhb( 'M', uplo, n, kd, ab, ldab, rwork )
341 iscale = 0
342 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
343 iscale = 1
344 sigma = rmin / anrm
345 ELSE IF( anrm.GT.rmax ) THEN
346 iscale = 1
347 sigma = rmax / anrm
348 END IF
349 IF( iscale.EQ.1 ) THEN
350 IF( lower ) THEN
351 CALL zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
352 ELSE
353 CALL zlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
354 END IF
355 END IF
356*
357* Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
358*
359 inde = 1
360 indwrk = inde + n
361 indwk2 = 1 + n*n
362 llwk2 = lwork - indwk2 + 1
363 llrwk = lrwork - indwrk + 1
364 CALL zhbtrd( jobz, uplo, n, kd, ab, ldab, w, rwork( inde ), z,
365 $ ldz, work, iinfo )
366*
367* For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEDC.
368*
369 IF( .NOT.wantz ) THEN
370 CALL dsterf( n, w, rwork( inde ), info )
371 ELSE
372 CALL zstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
373 $ llwk2, rwork( indwrk ), llrwk, iwork, liwork,
374 $ info )
375 CALL zgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
376 $ work( indwk2 ), n )
377 CALL zlacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
378 END IF
379*
380* If matrix was scaled, then rescale eigenvalues appropriately.
381*
382 IF( iscale.EQ.1 ) THEN
383 IF( info.EQ.0 ) THEN
384 imax = n
385 ELSE
386 imax = info - 1
387 END IF
388 CALL dscal( imax, one / sigma, w, 1 )
389 END IF
390*
391 work( 1 ) = lwmin
392 rwork( 1 ) = lrwmin
393 iwork( 1 ) = liwmin
394 RETURN
395*
396* End of ZHBEVD
397*
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zstedc
subroutine zstedc(compz, n, d, e, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
ZSTEDC
Definition zstedc.f:212
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:187

◆ zhbevd_2stage()

subroutine zhbevd_2stage ( character jobz,
character uplo,
integer n,
integer kd,
complex*16, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

ZHBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHBEVD_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
!> a complex Hermitian band matrix A using the 2stage technique for
!> the reduction to tridiagonal.  If eigenvectors are desired, it
!> uses a divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>                  Not available in this release.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 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 = (2KD+1)*N + KD*NTHREADS
!>                                   where KD is the size of the band.
!>                                   NTHREADS is the number of threads used when
!>                                   openMP compilation is enabled, otherwise =1.
!>          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK, 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 DOUBLE PRECISION array,
!>                                         dimension (LRWORK)
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of 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 array IWORK.
!>          If JOBZ = 'N' or 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, 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 257 of file zhbevd_2stage.f.

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

◆ zhbevx()

subroutine zhbevx ( character jobz,
character range,
character uplo,
integer n,
integer kd,
complex*16, dimension( ldab, * ) ab,
integer ldab,
complex*16, dimension( ldq, * ) q,
integer ldq,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

ZHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHBEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a complex Hermitian band 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]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]Q
!>          Q is COMPLEX*16 array, dimension (LDQ, N)
!>          If JOBZ = 'V', the N-by-N unitary matrix used in the
!>                          reduction to tridiagonal form.
!>          If JOBZ = 'N', the array Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.  If JOBZ = 'V', then
!>          LDQ >= max(1,N).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AB to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 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*16 array, dimension (N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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 264 of file zhbevx.f.

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

◆ zhbevx_2stage()

subroutine zhbevx_2stage ( character jobz,
character range,
character uplo,
integer n,
integer kd,
complex*16, dimension( ldab, * ) ab,
integer ldab,
complex*16, dimension( ldq, * ) q,
integer ldq,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

ZHBEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHBEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
!> of a complex Hermitian band 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]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]Q
!>          Q is COMPLEX*16 array, dimension (LDQ, N)
!>          If JOBZ = 'V', the N-by-N unitary matrix used in the
!>                          reduction to tridiagonal form.
!>          If JOBZ = 'N', the array Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.  If JOBZ = 'V', then
!>          LDQ >= max(1,N).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AB to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 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*16 array, dimension (LWORK)
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK. LWORK >= 1, when N <= 1;
!>          otherwise  
!>          If JOBZ = 'N' and N > 1, LWORK must be queried.
!>                                   LWORK = MAX(1, dimension) where
!>                                   dimension = (2KD+1)*N + KD*NTHREADS
!>                                   where KD is the size of the band.
!>                                   NTHREADS is the number of threads used when
!>                                   openMP compilation is enabled, otherwise =1.
!>          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK, 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 DOUBLE PRECISION 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 323 of file zhbevx_2stage.f.

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

◆ zhbgv()

subroutine zhbgv ( character jobz,
character uplo,
integer n,
integer ka,
integer kb,
complex*16, dimension( ldab, * ) ab,
integer ldab,
complex*16, dimension( ldbb, * ) bb,
integer ldbb,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZHBGV

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

Purpose:
!>
!> ZHBGV computes all the eigenvalues, and optionally, the eigenvectors
!> of a complex generalized Hermitian-definite banded eigenproblem, of
!> the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian
!> and banded, and B is also positive definite.
!>
Parameters
[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]KA
!>          KA is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'. KA >= 0.
!>
[in]KB
!>          KB is INTEGER
!>          The number of superdiagonals of the matrix B if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'. KB >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first ka+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
!>
!>          On exit, the contents of AB are destroyed.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KA+1.
!>
[in,out]BB
!>          BB is COMPLEX*16 array, dimension (LDBB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix B, stored in the first kb+1 rows of the array.  The
!>          j-th column of B is stored in the j-th column of the array BB
!>          as follows:
!>          if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
!>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
!>
!>          On exit, the factor S from the split Cholesky factorization
!>          B = S**H*S, as returned by ZPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
!>          eigenvectors, with the i-th column of Z holding the
!>          eigenvector associated with W(i). The eigenvectors are
!>          normalized so that Z**H*B*Z = I.
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= N.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (3*N)
!>
[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 i is:
!>             <= N:  the algorithm 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 ZPBSTF
!>                    returned INFO = i: 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 181 of file zhbgv.f.

183*
184* -- LAPACK driver routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER JOBZ, UPLO
190 INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, N
191* ..
192* .. Array Arguments ..
193 DOUBLE PRECISION RWORK( * ), W( * )
194 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
195 $ Z( LDZ, * )
196* ..
197*
198* =====================================================================
199*
200* .. Local Scalars ..
201 LOGICAL UPPER, WANTZ
202 CHARACTER VECT
203 INTEGER IINFO, INDE, INDWRK
204* ..
205* .. External Functions ..
206 LOGICAL LSAME
207 EXTERNAL lsame
208* ..
209* .. External Subroutines ..
210 EXTERNAL dsterf, xerbla, zhbgst, zhbtrd, zpbstf, zsteqr
211* ..
212* .. Executable Statements ..
213*
214* Test the input parameters.
215*
216 wantz = lsame( jobz, 'V' )
217 upper = lsame( uplo, 'U' )
218*
219 info = 0
220 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
221 info = -1
222 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
223 info = -2
224 ELSE IF( n.LT.0 ) THEN
225 info = -3
226 ELSE IF( ka.LT.0 ) THEN
227 info = -4
228 ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
229 info = -5
230 ELSE IF( ldab.LT.ka+1 ) THEN
231 info = -7
232 ELSE IF( ldbb.LT.kb+1 ) THEN
233 info = -9
234 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
235 info = -12
236 END IF
237 IF( info.NE.0 ) THEN
238 CALL xerbla( 'ZHBGV ', -info )
239 RETURN
240 END IF
241*
242* Quick return if possible
243*
244 IF( n.EQ.0 )
245 $ RETURN
246*
247* Form a split Cholesky factorization of B.
248*
249 CALL zpbstf( uplo, n, kb, bb, ldbb, info )
250 IF( info.NE.0 ) THEN
251 info = n + info
252 RETURN
253 END IF
254*
255* Transform problem to standard eigenvalue problem.
256*
257 inde = 1
258 indwrk = inde + n
259 CALL zhbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,
260 $ work, rwork( indwrk ), iinfo )
261*
262* Reduce to tridiagonal form.
263*
264 IF( wantz ) THEN
265 vect = 'U'
266 ELSE
267 vect = 'N'
268 END IF
269 CALL zhbtrd( vect, uplo, n, ka, ab, ldab, w, rwork( inde ), z,
270 $ ldz, work, iinfo )
271*
272* For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEQR.
273*
274 IF( .NOT.wantz ) THEN
275 CALL dsterf( n, w, rwork( inde ), info )
276 ELSE
277 CALL zsteqr( jobz, n, w, rwork( inde ), z, ldz,
278 $ rwork( indwrk ), info )
279 END IF
280 RETURN
281*
282* End of ZHBGV
283*
zhbgst
subroutine zhbgst(vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x, ldx, work, rwork, info)
ZHBGST
Definition zhbgst.f:165
zpbstf
subroutine zpbstf(uplo, n, kd, ab, ldab, info)
ZPBSTF
Definition zpbstf.f:153

◆ zhbgvd()

subroutine zhbgvd ( character jobz,
character uplo,
integer n,
integer ka,
integer kb,
complex*16, dimension( ldab, * ) ab,
integer ldab,
complex*16, dimension( ldbb, * ) bb,
integer ldbb,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

ZHBGVD

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

Purpose:
!>
!> ZHBGVD computes all the eigenvalues, and optionally, the eigenvectors
!> of a complex generalized Hermitian-definite banded eigenproblem, of
!> the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian
!> and banded, 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]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]KA
!>          KA is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'. KA >= 0.
!>
[in]KB
!>          KB is INTEGER
!>          The number of superdiagonals of the matrix B if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'. KB >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first ka+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
!>
!>          On exit, the contents of AB are destroyed.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KA+1.
!>
[in,out]BB
!>          BB is COMPLEX*16 array, dimension (LDBB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix B, stored in the first kb+1 rows of the array.  The
!>          j-th column of B is stored in the j-th column of the array BB
!>          as follows:
!>          if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
!>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
!>
!>          On exit, the factor S from the split Cholesky factorization
!>          B = S**H*S, as returned by ZPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
!>          eigenvectors, with the i-th column of Z holding the
!>          eigenvector associated with W(i). The eigenvectors are
!>          normalized so that Z**H*B*Z = I.
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= N.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If N <= 1,               LWORK >= 1.
!>          If JOBZ = 'N' and N > 1, LWORK >= N.
!>          If JOBZ = 'V' and N > 1, LWORK >= 2*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 DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO=0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of 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 array IWORK.
!>          If JOBZ = 'N' or 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:  if INFO = i, and i is:
!>             <= N:  the algorithm 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 ZPBSTF
!>                    returned INFO = i: 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 249 of file zhbgvd.f.

252*
253* -- LAPACK driver routine --
254* -- LAPACK is a software package provided by Univ. of Tennessee, --
255* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
256*
257* .. Scalar Arguments ..
258 CHARACTER JOBZ, UPLO
259 INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK,
260 $ LWORK, N
261* ..
262* .. Array Arguments ..
263 INTEGER IWORK( * )
264 DOUBLE PRECISION RWORK( * ), W( * )
265 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
266 $ Z( LDZ, * )
267* ..
268*
269* =====================================================================
270*
271* .. Parameters ..
272 COMPLEX*16 CONE, CZERO
273 parameter( cone = ( 1.0d+0, 0.0d+0 ),
274 $ czero = ( 0.0d+0, 0.0d+0 ) )
275* ..
276* .. Local Scalars ..
277 LOGICAL LQUERY, UPPER, WANTZ
278 CHARACTER VECT
279 INTEGER IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLRWK,
280 $ LLWK2, LRWMIN, LWMIN
281* ..
282* .. External Functions ..
283 LOGICAL LSAME
284 EXTERNAL lsame
285* ..
286* .. External Subroutines ..
287 EXTERNAL dsterf, xerbla, zgemm, zhbgst, zhbtrd, zlacpy,
288 $ zpbstf, zstedc
289* ..
290* .. Executable Statements ..
291*
292* Test the input parameters.
293*
294 wantz = lsame( jobz, 'V' )
295 upper = lsame( uplo, 'U' )
296 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
297*
298 info = 0
299 IF( n.LE.1 ) THEN
300 lwmin = 1+n
301 lrwmin = 1+n
302 liwmin = 1
303 ELSE IF( wantz ) THEN
304 lwmin = 2*n**2
305 lrwmin = 1 + 5*n + 2*n**2
306 liwmin = 3 + 5*n
307 ELSE
308 lwmin = n
309 lrwmin = n
310 liwmin = 1
311 END IF
312 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
313 info = -1
314 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
315 info = -2
316 ELSE IF( n.LT.0 ) THEN
317 info = -3
318 ELSE IF( ka.LT.0 ) THEN
319 info = -4
320 ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
321 info = -5
322 ELSE IF( ldab.LT.ka+1 ) THEN
323 info = -7
324 ELSE IF( ldbb.LT.kb+1 ) THEN
325 info = -9
326 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
327 info = -12
328 END IF
329*
330 IF( info.EQ.0 ) THEN
331 work( 1 ) = lwmin
332 rwork( 1 ) = lrwmin
333 iwork( 1 ) = liwmin
334*
335 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
336 info = -14
337 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
338 info = -16
339 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
340 info = -18
341 END IF
342 END IF
343*
344 IF( info.NE.0 ) THEN
345 CALL xerbla( 'ZHBGVD', -info )
346 RETURN
347 ELSE IF( lquery ) THEN
348 RETURN
349 END IF
350*
351* Quick return if possible
352*
353 IF( n.EQ.0 )
354 $ RETURN
355*
356* Form a split Cholesky factorization of B.
357*
358 CALL zpbstf( uplo, n, kb, bb, ldbb, info )
359 IF( info.NE.0 ) THEN
360 info = n + info
361 RETURN
362 END IF
363*
364* Transform problem to standard eigenvalue problem.
365*
366 inde = 1
367 indwrk = inde + n
368 indwk2 = 1 + n*n
369 llwk2 = lwork - indwk2 + 2
370 llrwk = lrwork - indwrk + 2
371 CALL zhbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,
372 $ work, rwork, iinfo )
373*
374* Reduce Hermitian band matrix to tridiagonal form.
375*
376 IF( wantz ) THEN
377 vect = 'U'
378 ELSE
379 vect = 'N'
380 END IF
381 CALL zhbtrd( vect, uplo, n, ka, ab, ldab, w, rwork( inde ), z,
382 $ ldz, work, iinfo )
383*
384* For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEDC.
385*
386 IF( .NOT.wantz ) THEN
387 CALL dsterf( n, w, rwork( inde ), info )
388 ELSE
389 CALL zstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
390 $ llwk2, rwork( indwrk ), llrwk, iwork, liwork,
391 $ info )
392 CALL zgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
393 $ work( indwk2 ), n )
394 CALL zlacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
395 END IF
396*
397 work( 1 ) = lwmin
398 rwork( 1 ) = lrwmin
399 iwork( 1 ) = liwmin
400 RETURN
401*
402* End of ZHBGVD
403*

◆ zhbgvx()

subroutine zhbgvx ( character jobz,
character range,
character uplo,
integer n,
integer ka,
integer kb,
complex*16, dimension( ldab, * ) ab,
integer ldab,
complex*16, dimension( ldbb, * ) bb,
integer ldbb,
complex*16, dimension( ldq, * ) q,
integer ldq,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

ZHBGVX

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

Purpose:
!>
!> ZHBGVX computes all the eigenvalues, and optionally, the eigenvectors
!> of a complex generalized Hermitian-definite banded eigenproblem, of
!> the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian
!> and banded, and B is also positive definite.  Eigenvalues and
!> eigenvectors can be selected by specifying either all eigenvalues,
!> 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 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]KA
!>          KA is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'. KA >= 0.
!>
[in]KB
!>          KB is INTEGER
!>          The number of superdiagonals of the matrix B if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'. KB >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first ka+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
!>
!>          On exit, the contents of AB are destroyed.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KA+1.
!>
[in,out]BB
!>          BB is COMPLEX*16 array, dimension (LDBB, N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix B, stored in the first kb+1 rows of the array.  The
!>          j-th column of B is stored in the j-th column of the array BB
!>          as follows:
!>          if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
!>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
!>
!>          On exit, the factor S from the split Cholesky factorization
!>          B = S**H*S, as returned by ZPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]Q
!>          Q is COMPLEX*16 array, dimension (LDQ, N)
!>          If JOBZ = 'V', the n-by-n matrix used in the reduction of
!>          A*x = (lambda)*B*x to standard form, i.e. C*x = (lambda)*x,
!>          and consequently C to tridiagonal form.
!>          If JOBZ = 'N', the array Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.  If JOBZ = 'N',
!>          LDQ >= 1. If JOBZ = 'V', LDQ >= max(1,N).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AP to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
!>          eigenvectors, with the i-th column of Z holding the
!>          eigenvector associated with W(i). The eigenvectors are
!>          normalized so that Z**H*B*Z = I.
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= N.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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, and i is:
!>             <= N:  then i eigenvectors failed to converge.  Their
!>                    indices are stored in array IFAIL.
!>             > N:   if INFO = N + i, for 1 <= i <= N, then ZPBSTF
!>                    returned INFO = i: 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 297 of file zhbgvx.f.

300*
301* -- LAPACK driver routine --
302* -- LAPACK is a software package provided by Univ. of Tennessee, --
303* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
304*
305* .. Scalar Arguments ..
306 CHARACTER JOBZ, RANGE, UPLO
307 INTEGER IL, INFO, IU, KA, KB, LDAB, LDBB, LDQ, LDZ, M,
308 $ N
309 DOUBLE PRECISION ABSTOL, VL, VU
310* ..
311* .. Array Arguments ..
312 INTEGER IFAIL( * ), IWORK( * )
313 DOUBLE PRECISION RWORK( * ), W( * )
314 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), Q( LDQ, * ),
315 $ WORK( * ), Z( LDZ, * )
316* ..
317*
318* =====================================================================
319*
320* .. Parameters ..
321 DOUBLE PRECISION ZERO
322 parameter( zero = 0.0d+0 )
323 COMPLEX*16 CZERO, CONE
324 parameter( czero = ( 0.0d+0, 0.0d+0 ),
325 $ cone = ( 1.0d+0, 0.0d+0 ) )
326* ..
327* .. Local Scalars ..
328 LOGICAL ALLEIG, INDEIG, TEST, UPPER, VALEIG, WANTZ
329 CHARACTER ORDER, VECT
330 INTEGER I, IINFO, INDD, INDE, INDEE, INDIBL, INDISP,
331 $ INDIWK, INDRWK, INDWRK, ITMP1, J, JJ, NSPLIT
332 DOUBLE PRECISION TMP1
333* ..
334* .. External Functions ..
335 LOGICAL LSAME
336 EXTERNAL lsame
337* ..
338* .. External Subroutines ..
339 EXTERNAL dcopy, dstebz, dsterf, xerbla, zcopy, zgemv,
340 $ zhbgst, zhbtrd, zlacpy, zpbstf, zstein, zsteqr,
341 $ zswap
342* ..
343* .. Intrinsic Functions ..
344 INTRINSIC 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*
356 info = 0
357 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
358 info = -1
359 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
360 info = -2
361 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
362 info = -3
363 ELSE IF( n.LT.0 ) THEN
364 info = -4
365 ELSE IF( ka.LT.0 ) THEN
366 info = -5
367 ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
368 info = -6
369 ELSE IF( ldab.LT.ka+1 ) THEN
370 info = -8
371 ELSE IF( ldbb.LT.kb+1 ) THEN
372 info = -10
373 ELSE IF( ldq.LT.1 .OR. ( wantz .AND. ldq.LT.n ) ) THEN
374 info = -12
375 ELSE
376 IF( valeig ) THEN
377 IF( n.GT.0 .AND. vu.LE.vl )
378 $ info = -14
379 ELSE IF( indeig ) THEN
380 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
381 info = -15
382 ELSE IF ( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
383 info = -16
384 END IF
385 END IF
386 END IF
387 IF( info.EQ.0) THEN
388 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
389 info = -21
390 END IF
391 END IF
392*
393 IF( info.NE.0 ) THEN
394 CALL xerbla( 'ZHBGVX', -info )
395 RETURN
396 END IF
397*
398* Quick return if possible
399*
400 m = 0
401 IF( n.EQ.0 )
402 $ RETURN
403*
404* Form a split Cholesky factorization of B.
405*
406 CALL zpbstf( uplo, n, kb, bb, ldbb, info )
407 IF( info.NE.0 ) THEN
408 info = n + info
409 RETURN
410 END IF
411*
412* Transform problem to standard eigenvalue problem.
413*
414 CALL zhbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, q, ldq,
415 $ work, rwork, iinfo )
416*
417* Solve the standard eigenvalue problem.
418* Reduce Hermitian band matrix to tridiagonal form.
419*
420 indd = 1
421 inde = indd + n
422 indrwk = inde + n
423 indwrk = 1
424 IF( wantz ) THEN
425 vect = 'U'
426 ELSE
427 vect = 'N'
428 END IF
429 CALL zhbtrd( vect, uplo, n, ka, ab, ldab, rwork( indd ),
430 $ rwork( inde ), q, ldq, work( indwrk ), iinfo )
431*
432* If all eigenvalues are desired and ABSTOL is less than or equal
433* to zero, then call DSTERF or ZSTEQR. If this fails for some
434* eigenvalue, then try DSTEBZ.
435*
436 test = .false.
437 IF( indeig ) THEN
438 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
439 test = .true.
440 END IF
441 END IF
442 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
443 CALL dcopy( n, rwork( indd ), 1, w, 1 )
444 indee = indrwk + 2*n
445 CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
446 IF( .NOT.wantz ) THEN
447 CALL dsterf( n, w, rwork( indee ), info )
448 ELSE
449 CALL zlacpy( 'A', n, n, q, ldq, z, ldz )
450 CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
451 $ rwork( indrwk ), info )
452 IF( info.EQ.0 ) THEN
453 DO 10 i = 1, n
454 ifail( i ) = 0
455 10 CONTINUE
456 END IF
457 END IF
458 IF( info.EQ.0 ) THEN
459 m = n
460 GO TO 30
461 END IF
462 info = 0
463 END IF
464*
465* Otherwise, call DSTEBZ and, if eigenvectors are desired,
466* call ZSTEIN.
467*
468 IF( wantz ) THEN
469 order = 'B'
470 ELSE
471 order = 'E'
472 END IF
473 indibl = 1
474 indisp = indibl + n
475 indiwk = indisp + n
476 CALL dstebz( range, order, n, vl, vu, il, iu, abstol,
477 $ rwork( indd ), rwork( inde ), m, nsplit, w,
478 $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
479 $ iwork( indiwk ), info )
480*
481 IF( wantz ) THEN
482 CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
483 $ iwork( indibl ), iwork( indisp ), z, ldz,
484 $ rwork( indrwk ), iwork( indiwk ), ifail, info )
485*
486* Apply unitary matrix used in reduction to tridiagonal
487* form to eigenvectors returned by ZSTEIN.
488*
489 DO 20 j = 1, m
490 CALL zcopy( n, z( 1, j ), 1, work( 1 ), 1 )
491 CALL zgemv( 'N', n, n, cone, q, ldq, work, 1, czero,
492 $ z( 1, j ), 1 )
493 20 CONTINUE
494 END IF
495*
496 30 CONTINUE
497*
498* If eigenvalues are not in order, then sort them, along with
499* eigenvectors.
500*
501 IF( wantz ) THEN
502 DO 50 j = 1, m - 1
503 i = 0
504 tmp1 = w( j )
505 DO 40 jj = j + 1, m
506 IF( w( jj ).LT.tmp1 ) THEN
507 i = jj
508 tmp1 = w( jj )
509 END IF
510 40 CONTINUE
511*
512 IF( i.NE.0 ) THEN
513 itmp1 = iwork( indibl+i-1 )
514 w( i ) = w( j )
515 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
516 w( j ) = tmp1
517 iwork( indibl+j-1 ) = itmp1
518 CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
519 IF( info.NE.0 ) THEN
520 itmp1 = ifail( i )
521 ifail( i ) = ifail( j )
522 ifail( j ) = itmp1
523 END IF
524 END IF
525 50 CONTINUE
526 END IF
527*
528 RETURN
529*
530* End of ZHBGVX
531*

◆ zhpev()

subroutine zhpev ( character jobz,
character uplo,
integer n,
complex*16, dimension( * ) ap,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZHPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHPEV computes all the eigenvalues and, optionally, eigenvectors of a
!> complex Hermitian matrix in packed storage.
!>
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]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, AP is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
!>          and first superdiagonal of the tridiagonal matrix T overwrite
!>          the corresponding elements of A, and if UPLO = 'L', the
!>          diagonal and first subdiagonal of T overwrite the
!>          corresponding elements of A.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (max(1, 2*N-1))
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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 136 of file zhpev.f.

138*
139* -- LAPACK driver routine --
140* -- LAPACK is a software package provided by Univ. of Tennessee, --
141* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142*
143* .. Scalar Arguments ..
144 CHARACTER JOBZ, UPLO
145 INTEGER INFO, LDZ, N
146* ..
147* .. Array Arguments ..
148 DOUBLE PRECISION RWORK( * ), W( * )
149 COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * )
150* ..
151*
152* =====================================================================
153*
154* .. Parameters ..
155 DOUBLE PRECISION ZERO, ONE
156 parameter( zero = 0.0d0, one = 1.0d0 )
157* ..
158* .. Local Scalars ..
159 LOGICAL WANTZ
160 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWRK,
161 $ ISCALE
162 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
163 $ SMLNUM
164* ..
165* .. External Functions ..
166 LOGICAL LSAME
167 DOUBLE PRECISION DLAMCH, ZLANHP
168 EXTERNAL lsame, dlamch, zlanhp
169* ..
170* .. External Subroutines ..
171 EXTERNAL dscal, dsterf, xerbla, zdscal, zhptrd, zsteqr,
172 $ zupgtr
173* ..
174* .. Intrinsic Functions ..
175 INTRINSIC sqrt
176* ..
177* .. Executable Statements ..
178*
179* Test the input parameters.
180*
181 wantz = lsame( jobz, 'V' )
182*
183 info = 0
184 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
185 info = -1
186 ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
187 $ THEN
188 info = -2
189 ELSE IF( n.LT.0 ) THEN
190 info = -3
191 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
192 info = -7
193 END IF
194*
195 IF( info.NE.0 ) THEN
196 CALL xerbla( 'ZHPEV ', -info )
197 RETURN
198 END IF
199*
200* Quick return if possible
201*
202 IF( n.EQ.0 )
203 $ RETURN
204*
205 IF( n.EQ.1 ) THEN
206 w( 1 ) = dble( ap( 1 ) )
207 rwork( 1 ) = 1
208 IF( wantz )
209 $ z( 1, 1 ) = one
210 RETURN
211 END IF
212*
213* Get machine constants.
214*
215 safmin = dlamch( 'Safe minimum' )
216 eps = dlamch( 'Precision' )
217 smlnum = safmin / eps
218 bignum = one / smlnum
219 rmin = sqrt( smlnum )
220 rmax = sqrt( bignum )
221*
222* Scale matrix to allowable range, if necessary.
223*
224 anrm = zlanhp( 'M', uplo, n, ap, rwork )
225 iscale = 0
226 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
227 iscale = 1
228 sigma = rmin / anrm
229 ELSE IF( anrm.GT.rmax ) THEN
230 iscale = 1
231 sigma = rmax / anrm
232 END IF
233 IF( iscale.EQ.1 ) THEN
234 CALL zdscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
235 END IF
236*
237* Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form.
238*
239 inde = 1
240 indtau = 1
241 CALL zhptrd( uplo, n, ap, w, rwork( inde ), work( indtau ),
242 $ iinfo )
243*
244* For eigenvalues only, call DSTERF. For eigenvectors, first call
245* ZUPGTR to generate the orthogonal matrix, then call ZSTEQR.
246*
247 IF( .NOT.wantz ) THEN
248 CALL dsterf( n, w, rwork( inde ), info )
249 ELSE
250 indwrk = indtau + n
251 CALL zupgtr( uplo, n, ap, work( indtau ), z, ldz,
252 $ work( indwrk ), iinfo )
253 indrwk = inde + n
254 CALL zsteqr( jobz, n, w, rwork( inde ), z, ldz,
255 $ rwork( indrwk ), info )
256 END IF
257*
258* If matrix was scaled, then rescale eigenvalues appropriately.
259*
260 IF( iscale.EQ.1 ) THEN
261 IF( info.EQ.0 ) THEN
262 imax = n
263 ELSE
264 imax = info - 1
265 END IF
266 CALL dscal( imax, one / sigma, w, 1 )
267 END IF
268*
269 RETURN
270*
271* End of ZHPEV
272*
zlanhp
double precision function zlanhp(norm, uplo, n, ap, work)
ZLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlanhp.f:117
zupgtr
subroutine zupgtr(uplo, n, ap, tau, q, ldq, work, info)
ZUPGTR
Definition zupgtr.f:114
zhptrd
subroutine zhptrd(uplo, n, ap, d, e, tau, info)
ZHPTRD
Definition zhptrd.f:151
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78

◆ zhpevd()

subroutine zhpevd ( character jobz,
character uplo,
integer n,
complex*16, dimension( * ) ap,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

ZHPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHPEVD computes all the eigenvalues and, optionally, eigenvectors of
!> a complex Hermitian matrix A in packed storage.  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]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, AP is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
!>          and first superdiagonal of the tridiagonal matrix T overwrite
!>          the corresponding elements of A, and if UPLO = 'L', the
!>          diagonal and first subdiagonal of T overwrite the
!>          corresponding elements of A.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the required LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of array WORK.
!>          If N <= 1,               LWORK must be at least 1.
!>          If JOBZ = 'N' and N > 1, LWORK must be at least N.
!>          If JOBZ = 'V' and N > 1, LWORK must be at least 2*N.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the required 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 DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the required LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of 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 required 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 required LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of array IWORK.
!>          If JOBZ  = 'N' or 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 required 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, 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 198 of file zhpevd.f.

200*
201* -- LAPACK driver routine --
202* -- LAPACK is a software package provided by Univ. of Tennessee, --
203* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
204*
205* .. Scalar Arguments ..
206 CHARACTER JOBZ, UPLO
207 INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
208* ..
209* .. Array Arguments ..
210 INTEGER IWORK( * )
211 DOUBLE PRECISION RWORK( * ), W( * )
212 COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * )
213* ..
214*
215* =====================================================================
216*
217* .. Parameters ..
218 DOUBLE PRECISION ZERO, ONE
219 parameter( zero = 0.0d+0, one = 1.0d+0 )
220 COMPLEX*16 CONE
221 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
222* ..
223* .. Local Scalars ..
224 LOGICAL LQUERY, WANTZ
225 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWRK,
226 $ ISCALE, LIWMIN, LLRWK, LLWRK, LRWMIN, LWMIN
227 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
228 $ SMLNUM
229* ..
230* .. External Functions ..
231 LOGICAL LSAME
232 DOUBLE PRECISION DLAMCH, ZLANHP
233 EXTERNAL lsame, dlamch, zlanhp
234* ..
235* .. External Subroutines ..
236 EXTERNAL dscal, dsterf, xerbla, zdscal, zhptrd, zstedc,
237 $ zupmtr
238* ..
239* .. Intrinsic Functions ..
240 INTRINSIC sqrt
241* ..
242* .. Executable Statements ..
243*
244* Test the input parameters.
245*
246 wantz = lsame( jobz, 'V' )
247 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
248*
249 info = 0
250 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
251 info = -1
252 ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
253 $ THEN
254 info = -2
255 ELSE IF( n.LT.0 ) THEN
256 info = -3
257 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
258 info = -7
259 END IF
260*
261 IF( info.EQ.0 ) THEN
262 IF( n.LE.1 ) THEN
263 lwmin = 1
264 liwmin = 1
265 lrwmin = 1
266 ELSE
267 IF( wantz ) THEN
268 lwmin = 2*n
269 lrwmin = 1 + 5*n + 2*n**2
270 liwmin = 3 + 5*n
271 ELSE
272 lwmin = n
273 lrwmin = n
274 liwmin = 1
275 END IF
276 END IF
277 work( 1 ) = lwmin
278 rwork( 1 ) = lrwmin
279 iwork( 1 ) = liwmin
280*
281 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
282 info = -9
283 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
284 info = -11
285 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
286 info = -13
287 END IF
288 END IF
289*
290 IF( info.NE.0 ) THEN
291 CALL xerbla( 'ZHPEVD', -info )
292 RETURN
293 ELSE IF( lquery ) THEN
294 RETURN
295 END IF
296*
297* Quick return if possible
298*
299 IF( n.EQ.0 )
300 $ RETURN
301*
302 IF( n.EQ.1 ) THEN
303 w( 1 ) = dble( ap( 1 ) )
304 IF( wantz )
305 $ z( 1, 1 ) = cone
306 RETURN
307 END IF
308*
309* Get machine constants.
310*
311 safmin = dlamch( 'Safe minimum' )
312 eps = dlamch( 'Precision' )
313 smlnum = safmin / eps
314 bignum = one / smlnum
315 rmin = sqrt( smlnum )
316 rmax = sqrt( bignum )
317*
318* Scale matrix to allowable range, if necessary.
319*
320 anrm = zlanhp( 'M', uplo, n, ap, rwork )
321 iscale = 0
322 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
323 iscale = 1
324 sigma = rmin / anrm
325 ELSE IF( anrm.GT.rmax ) THEN
326 iscale = 1
327 sigma = rmax / anrm
328 END IF
329 IF( iscale.EQ.1 ) THEN
330 CALL zdscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
331 END IF
332*
333* Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form.
334*
335 inde = 1
336 indtau = 1
337 indrwk = inde + n
338 indwrk = indtau + n
339 llwrk = lwork - indwrk + 1
340 llrwk = lrwork - indrwk + 1
341 CALL zhptrd( uplo, n, ap, w, rwork( inde ), work( indtau ),
342 $ iinfo )
343*
344* For eigenvalues only, call DSTERF. For eigenvectors, first call
345* ZUPGTR to generate the orthogonal matrix, then call ZSTEDC.
346*
347 IF( .NOT.wantz ) THEN
348 CALL dsterf( n, w, rwork( inde ), info )
349 ELSE
350 CALL zstedc( 'I', n, w, rwork( inde ), z, ldz, work( indwrk ),
351 $ llwrk, rwork( indrwk ), llrwk, iwork, liwork,
352 $ info )
353 CALL zupmtr( 'L', uplo, 'N', n, n, ap, work( indtau ), z, ldz,
354 $ work( indwrk ), iinfo )
355 END IF
356*
357* If matrix was scaled, then rescale eigenvalues appropriately.
358*
359 IF( iscale.EQ.1 ) THEN
360 IF( info.EQ.0 ) THEN
361 imax = n
362 ELSE
363 imax = info - 1
364 END IF
365 CALL dscal( imax, one / sigma, w, 1 )
366 END IF
367*
368 work( 1 ) = lwmin
369 rwork( 1 ) = lrwmin
370 iwork( 1 ) = liwmin
371 RETURN
372*
373* End of ZHPEVD
374*
zupmtr
subroutine zupmtr(side, uplo, trans, m, n, ap, tau, c, ldc, work, info)
ZUPMTR
Definition zupmtr.f:150

◆ zhpevx()

subroutine zhpevx ( character jobz,
character range,
character uplo,
integer n,
complex*16, dimension( * ) ap,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

ZHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
!>
!> ZHPEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a complex Hermitian matrix A in packed storage.
!> Eigenvalues/vectors 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]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, AP is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
!>          and first superdiagonal of the tridiagonal matrix T overwrite
!>          the corresponding elements of A, and if UPLO = 'L', the
!>          diagonal and first subdiagonal of T overwrite the
!>          corresponding elements of A.
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AP to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the selected eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 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*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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 237 of file zhpevx.f.

240*
241* -- LAPACK driver routine --
242* -- LAPACK is a software package provided by Univ. of Tennessee, --
243* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
244*
245* .. Scalar Arguments ..
246 CHARACTER JOBZ, RANGE, UPLO
247 INTEGER IL, INFO, IU, LDZ, M, N
248 DOUBLE PRECISION ABSTOL, VL, VU
249* ..
250* .. Array Arguments ..
251 INTEGER IFAIL( * ), IWORK( * )
252 DOUBLE PRECISION RWORK( * ), W( * )
253 COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * )
254* ..
255*
256* =====================================================================
257*
258* .. Parameters ..
259 DOUBLE PRECISION ZERO, ONE
260 parameter( zero = 0.0d0, one = 1.0d0 )
261 COMPLEX*16 CONE
262 parameter( cone = ( 1.0d0, 0.0d0 ) )
263* ..
264* .. Local Scalars ..
265 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
266 CHARACTER ORDER
267 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
268 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
269 $ ITMP1, J, JJ, NSPLIT
270 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
271 $ SIGMA, SMLNUM, TMP1, VLL, VUU
272* ..
273* .. External Functions ..
274 LOGICAL LSAME
275 DOUBLE PRECISION DLAMCH, ZLANHP
276 EXTERNAL lsame, dlamch, zlanhp
277* ..
278* .. External Subroutines ..
279 EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
280 $ zhptrd, zstein, zsteqr, zswap, zupgtr, zupmtr
281* ..
282* .. Intrinsic Functions ..
283 INTRINSIC dble, max, min, sqrt
284* ..
285* .. Executable Statements ..
286*
287* Test the input parameters.
288*
289 wantz = lsame( jobz, 'V' )
290 alleig = lsame( range, 'A' )
291 valeig = lsame( range, 'V' )
292 indeig = lsame( range, 'I' )
293*
294 info = 0
295 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
296 info = -1
297 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
298 info = -2
299 ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
300 $ THEN
301 info = -3
302 ELSE IF( n.LT.0 ) THEN
303 info = -4
304 ELSE
305 IF( valeig ) THEN
306 IF( n.GT.0 .AND. vu.LE.vl )
307 $ info = -7
308 ELSE IF( indeig ) THEN
309 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
310 info = -8
311 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
312 info = -9
313 END IF
314 END IF
315 END IF
316 IF( info.EQ.0 ) THEN
317 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
318 $ info = -14
319 END IF
320*
321 IF( info.NE.0 ) THEN
322 CALL xerbla( 'ZHPEVX', -info )
323 RETURN
324 END IF
325*
326* Quick return if possible
327*
328 m = 0
329 IF( n.EQ.0 )
330 $ RETURN
331*
332 IF( n.EQ.1 ) THEN
333 IF( alleig .OR. indeig ) THEN
334 m = 1
335 w( 1 ) = dble( ap( 1 ) )
336 ELSE
337 IF( vl.LT.dble( ap( 1 ) ) .AND. vu.GE.dble( ap( 1 ) ) ) THEN
338 m = 1
339 w( 1 ) = dble( ap( 1 ) )
340 END IF
341 END IF
342 IF( wantz )
343 $ z( 1, 1 ) = cone
344 RETURN
345 END IF
346*
347* Get machine constants.
348*
349 safmin = dlamch( 'Safe minimum' )
350 eps = dlamch( 'Precision' )
351 smlnum = safmin / eps
352 bignum = one / smlnum
353 rmin = sqrt( smlnum )
354 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
355*
356* Scale matrix to allowable range, if necessary.
357*
358 iscale = 0
359 abstll = abstol
360 IF( valeig ) THEN
361 vll = vl
362 vuu = vu
363 ELSE
364 vll = zero
365 vuu = zero
366 END IF
367 anrm = zlanhp( 'M', uplo, n, ap, rwork )
368 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
369 iscale = 1
370 sigma = rmin / anrm
371 ELSE IF( anrm.GT.rmax ) THEN
372 iscale = 1
373 sigma = rmax / anrm
374 END IF
375 IF( iscale.EQ.1 ) THEN
376 CALL zdscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
377 IF( abstol.GT.0 )
378 $ abstll = abstol*sigma
379 IF( valeig ) THEN
380 vll = vl*sigma
381 vuu = vu*sigma
382 END IF
383 END IF
384*
385* Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form.
386*
387 indd = 1
388 inde = indd + n
389 indrwk = inde + n
390 indtau = 1
391 indwrk = indtau + n
392 CALL zhptrd( uplo, n, ap, rwork( indd ), rwork( inde ),
393 $ work( indtau ), iinfo )
394*
395* If all eigenvalues are desired and ABSTOL is less than or equal
396* to zero, then call DSTERF or ZUPGTR and ZSTEQR. If this fails
397* for some eigenvalue, then try DSTEBZ.
398*
399 test = .false.
400 IF (indeig) THEN
401 IF (il.EQ.1 .AND. iu.EQ.n) THEN
402 test = .true.
403 END IF
404 END IF
405 IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
406 CALL dcopy( n, rwork( indd ), 1, w, 1 )
407 indee = indrwk + 2*n
408 IF( .NOT.wantz ) THEN
409 CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
410 CALL dsterf( n, w, rwork( indee ), info )
411 ELSE
412 CALL zupgtr( uplo, n, ap, work( indtau ), z, ldz,
413 $ work( indwrk ), iinfo )
414 CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
415 CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
416 $ rwork( indrwk ), info )
417 IF( info.EQ.0 ) THEN
418 DO 10 i = 1, n
419 ifail( i ) = 0
420 10 CONTINUE
421 END IF
422 END IF
423 IF( info.EQ.0 ) THEN
424 m = n
425 GO TO 20
426 END IF
427 info = 0
428 END IF
429*
430* Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
431*
432 IF( wantz ) THEN
433 order = 'B'
434 ELSE
435 order = 'E'
436 END IF
437 indibl = 1
438 indisp = indibl + n
439 indiwk = indisp + n
440 CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
441 $ rwork( indd ), rwork( inde ), m, nsplit, w,
442 $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
443 $ iwork( indiwk ), info )
444*
445 IF( wantz ) THEN
446 CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
447 $ iwork( indibl ), iwork( indisp ), z, ldz,
448 $ rwork( indrwk ), iwork( indiwk ), ifail, info )
449*
450* Apply unitary matrix used in reduction to tridiagonal
451* form to eigenvectors returned by ZSTEIN.
452*
453 indwrk = indtau + n
454 CALL zupmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
455 $ work( indwrk ), iinfo )
456 END IF
457*
458* If matrix was scaled, then rescale eigenvalues appropriately.
459*
460 20 CONTINUE
461 IF( iscale.EQ.1 ) THEN
462 IF( info.EQ.0 ) THEN
463 imax = m
464 ELSE
465 imax = info - 1
466 END IF
467 CALL dscal( imax, one / sigma, w, 1 )
468 END IF
469*
470* If eigenvalues are not in order, then sort them, along with
471* eigenvectors.
472*
473 IF( wantz ) THEN
474 DO 40 j = 1, m - 1
475 i = 0
476 tmp1 = w( j )
477 DO 30 jj = j + 1, m
478 IF( w( jj ).LT.tmp1 ) THEN
479 i = jj
480 tmp1 = w( jj )
481 END IF
482 30 CONTINUE
483*
484 IF( i.NE.0 ) THEN
485 itmp1 = iwork( indibl+i-1 )
486 w( i ) = w( j )
487 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
488 w( j ) = tmp1
489 iwork( indibl+j-1 ) = itmp1
490 CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
491 IF( info.NE.0 ) THEN
492 itmp1 = ifail( i )
493 ifail( i ) = ifail( j )
494 ifail( j ) = itmp1
495 END IF
496 END IF
497 40 CONTINUE
498 END IF
499*
500 RETURN
501*
502* End of ZHPEVX
503*

◆ zhpgv()

subroutine zhpgv ( integer itype,
character jobz,
character uplo,
integer n,
complex*16, dimension( * ) ap,
complex*16, dimension( * ) bp,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZHPGV

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

Purpose:
!>
!> ZHPGV 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, stored in packed format,
!> 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]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, the contents of AP are destroyed.
!>
[in,out]BP
!>          BP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          B, packed columnwise in a linear array.  The j-th column of B
!>          is stored in the array BP as follows:
!>          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
!>          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
!>
!>          On exit, the triangular factor U or L from the Cholesky
!>          factorization B = U**H*U or B = L*L**H, in the same storage
!>          format as B.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z 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 Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (max(1, 2*N-1))
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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:  ZPPTRF or ZHPEV returned an error code:
!>             <= N:  if INFO = i, ZHPEV failed to converge;
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not convergeto 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 163 of file zhpgv.f.

165*
166* -- LAPACK driver routine --
167* -- LAPACK is a software package provided by Univ. of Tennessee, --
168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169*
170* .. Scalar Arguments ..
171 CHARACTER JOBZ, UPLO
172 INTEGER INFO, ITYPE, LDZ, N
173* ..
174* .. Array Arguments ..
175 DOUBLE PRECISION RWORK( * ), W( * )
176 COMPLEX*16 AP( * ), BP( * ), WORK( * ), Z( LDZ, * )
177* ..
178*
179* =====================================================================
180*
181* .. Local Scalars ..
182 LOGICAL UPPER, WANTZ
183 CHARACTER TRANS
184 INTEGER J, NEIG
185* ..
186* .. External Functions ..
187 LOGICAL LSAME
188 EXTERNAL lsame
189* ..
190* .. External Subroutines ..
191 EXTERNAL xerbla, zhpev, zhpgst, zpptrf, ztpmv, ztpsv
192* ..
193* .. Executable Statements ..
194*
195* Test the input parameters.
196*
197 wantz = lsame( jobz, 'V' )
198 upper = lsame( uplo, 'U' )
199*
200 info = 0
201 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
202 info = -1
203 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
204 info = -2
205 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
206 info = -3
207 ELSE IF( n.LT.0 ) THEN
208 info = -4
209 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
210 info = -9
211 END IF
212 IF( info.NE.0 ) THEN
213 CALL xerbla( 'ZHPGV ', -info )
214 RETURN
215 END IF
216*
217* Quick return if possible
218*
219 IF( n.EQ.0 )
220 $ RETURN
221*
222* Form a Cholesky factorization of B.
223*
224 CALL zpptrf( uplo, n, bp, info )
225 IF( info.NE.0 ) THEN
226 info = n + info
227 RETURN
228 END IF
229*
230* Transform problem to standard eigenvalue problem and solve.
231*
232 CALL zhpgst( itype, uplo, n, ap, bp, info )
233 CALL zhpev( jobz, uplo, n, ap, w, z, ldz, work, rwork, info )
234*
235 IF( wantz ) THEN
236*
237* Backtransform eigenvectors to the original problem.
238*
239 neig = n
240 IF( info.GT.0 )
241 $ neig = info - 1
242 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
243*
244* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
245* backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
246*
247 IF( upper ) THEN
248 trans = 'N'
249 ELSE
250 trans = 'C'
251 END IF
252*
253 DO 10 j = 1, neig
254 CALL ztpsv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
255 $ 1 )
256 10 CONTINUE
257*
258 ELSE IF( itype.EQ.3 ) THEN
259*
260* For B*A*x=(lambda)*x;
261* backtransform eigenvectors: x = L*y or U**H *y
262*
263 IF( upper ) THEN
264 trans = 'C'
265 ELSE
266 trans = 'N'
267 END IF
268*
269 DO 20 j = 1, neig
270 CALL ztpmv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
271 $ 1 )
272 20 CONTINUE
273 END IF
274 END IF
275 RETURN
276*
277* End of ZHPGV
278*
zhpgst
subroutine zhpgst(itype, uplo, n, ap, bp, info)
ZHPGST
Definition zhpgst.f:113
zpptrf
subroutine zpptrf(uplo, n, ap, info)
ZPPTRF
Definition zpptrf.f:119
zhpev
subroutine zhpev(jobz, uplo, n, ap, w, z, ldz, work, rwork, info)
ZHPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
Definition zhpev.f:138
ztpsv
subroutine ztpsv(uplo, trans, diag, n, ap, x, incx)
ZTPSV
Definition ztpsv.f:144
ztpmv
subroutine ztpmv(uplo, trans, diag, n, ap, x, incx)
ZTPMV
Definition ztpmv.f:142

◆ zhpgvd()

subroutine zhpgvd ( integer itype,
character jobz,
character uplo,
integer n,
complex*16, dimension( * ) ap,
complex*16, dimension( * ) bp,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

ZHPGVD

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

Purpose:
!>
!> ZHPGVD 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, stored in packed format, 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]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, the contents of AP are destroyed.
!>
[in,out]BP
!>          BP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          B, packed columnwise in a linear array.  The j-th column of B
!>          is stored in the array BP as follows:
!>          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
!>          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
!>
!>          On exit, the triangular factor U or L from the Cholesky
!>          factorization B = U**H*U or B = L*L**H, in the same storage
!>          format as B.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z 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 Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the required LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If N <= 1,               LWORK >= 1.
!>          If JOBZ = 'N' and N > 1, LWORK >= N.
!>          If JOBZ = 'V' and N > 1, LWORK >= 2*N.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the required 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 DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the required LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of 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 required 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 required LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of array IWORK.
!>          If JOBZ  = 'N' or 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 required 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:  ZPPTRF or ZHPEVD returned an error code:
!>             <= N:  if INFO = i, ZHPEVD failed to converge;
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not convergeto 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.
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 229 of file zhpgvd.f.

231*
232* -- LAPACK driver routine --
233* -- LAPACK is a software package provided by Univ. of Tennessee, --
234* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
235*
236* .. Scalar Arguments ..
237 CHARACTER JOBZ, UPLO
238 INTEGER INFO, ITYPE, LDZ, LIWORK, LRWORK, LWORK, N
239* ..
240* .. Array Arguments ..
241 INTEGER IWORK( * )
242 DOUBLE PRECISION RWORK( * ), W( * )
243 COMPLEX*16 AP( * ), BP( * ), WORK( * ), Z( LDZ, * )
244* ..
245*
246* =====================================================================
247*
248* .. Local Scalars ..
249 LOGICAL LQUERY, UPPER, WANTZ
250 CHARACTER TRANS
251 INTEGER J, LIWMIN, LRWMIN, LWMIN, NEIG
252* ..
253* .. External Functions ..
254 LOGICAL LSAME
255 EXTERNAL lsame
256* ..
257* .. External Subroutines ..
258 EXTERNAL xerbla, zhpevd, zhpgst, zpptrf, ztpmv, ztpsv
259* ..
260* .. Intrinsic Functions ..
261 INTRINSIC dble, max
262* ..
263* .. Executable Statements ..
264*
265* Test the input parameters.
266*
267 wantz = lsame( jobz, 'V' )
268 upper = lsame( uplo, 'U' )
269 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
270*
271 info = 0
272 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
273 info = -1
274 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
275 info = -2
276 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
277 info = -3
278 ELSE IF( n.LT.0 ) THEN
279 info = -4
280 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
281 info = -9
282 END IF
283*
284 IF( info.EQ.0 ) THEN
285 IF( n.LE.1 ) THEN
286 lwmin = 1
287 liwmin = 1
288 lrwmin = 1
289 ELSE
290 IF( wantz ) THEN
291 lwmin = 2*n
292 lrwmin = 1 + 5*n + 2*n**2
293 liwmin = 3 + 5*n
294 ELSE
295 lwmin = n
296 lrwmin = n
297 liwmin = 1
298 END IF
299 END IF
300*
301 work( 1 ) = lwmin
302 rwork( 1 ) = lrwmin
303 iwork( 1 ) = liwmin
304 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
305 info = -11
306 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
307 info = -13
308 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
309 info = -15
310 END IF
311 END IF
312*
313 IF( info.NE.0 ) THEN
314 CALL xerbla( 'ZHPGVD', -info )
315 RETURN
316 ELSE IF( lquery ) THEN
317 RETURN
318 END IF
319*
320* Quick return if possible
321*
322 IF( n.EQ.0 )
323 $ RETURN
324*
325* Form a Cholesky factorization of B.
326*
327 CALL zpptrf( uplo, n, bp, info )
328 IF( info.NE.0 ) THEN
329 info = n + info
330 RETURN
331 END IF
332*
333* Transform problem to standard eigenvalue problem and solve.
334*
335 CALL zhpgst( itype, uplo, n, ap, bp, info )
336 CALL zhpevd( jobz, uplo, n, ap, w, z, ldz, work, lwork, rwork,
337 $ lrwork, iwork, liwork, info )
338 lwmin = max( dble( lwmin ), dble( work( 1 ) ) )
339 lrwmin = max( dble( lrwmin ), dble( rwork( 1 ) ) )
340 liwmin = max( dble( liwmin ), dble( iwork( 1 ) ) )
341*
342 IF( wantz ) THEN
343*
344* Backtransform eigenvectors to the original problem.
345*
346 neig = n
347 IF( info.GT.0 )
348 $ neig = info - 1
349 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
350*
351* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
352* backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
353*
354 IF( upper ) THEN
355 trans = 'N'
356 ELSE
357 trans = 'C'
358 END IF
359*
360 DO 10 j = 1, neig
361 CALL ztpsv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
362 $ 1 )
363 10 CONTINUE
364*
365 ELSE IF( itype.EQ.3 ) THEN
366*
367* For B*A*x=(lambda)*x;
368* backtransform eigenvectors: x = L*y or U**H *y
369*
370 IF( upper ) THEN
371 trans = 'C'
372 ELSE
373 trans = 'N'
374 END IF
375*
376 DO 20 j = 1, neig
377 CALL ztpmv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
378 $ 1 )
379 20 CONTINUE
380 END IF
381 END IF
382*
383 work( 1 ) = lwmin
384 rwork( 1 ) = lrwmin
385 iwork( 1 ) = liwmin
386 RETURN
387*
388* End of ZHPGVD
389*
zhpevd
subroutine zhpevd(jobz, uplo, n, ap, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
ZHPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition zhpevd.f:200

◆ zhpgvx()

subroutine zhpgvx ( integer itype,
character jobz,
character range,
character uplo,
integer n,
complex*16, dimension( * ) ap,
complex*16, dimension( * ) bp,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

ZHPGVX

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

Purpose:
!>
!> ZHPGVX 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, stored in packed format, 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]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, the contents of AP are destroyed.
!>
[in,out]BP
!>          BP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          B, packed columnwise in a linear array.  The j-th column of B
!>          is stored in the array BP as follows:
!>          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
!>          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
!>
!>          On exit, the triangular factor U or L from the Cholesky
!>          factorization B = U**H*U or B = L*L**H, in the same storage
!>          format as B.
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AP to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          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**H*B*Z = I;
!>          if ITYPE = 3, Z**H*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*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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:  ZPPTRF or ZHPEVX returned an error code:
!>             <= N:  if INFO = i, ZHPEVX 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 274 of file zhpgvx.f.

277*
278* -- LAPACK driver routine --
279* -- LAPACK is a software package provided by Univ. of Tennessee, --
280* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
281*
282* .. Scalar Arguments ..
283 CHARACTER JOBZ, RANGE, UPLO
284 INTEGER IL, INFO, ITYPE, IU, LDZ, M, N
285 DOUBLE PRECISION ABSTOL, VL, VU
286* ..
287* .. Array Arguments ..
288 INTEGER IFAIL( * ), IWORK( * )
289 DOUBLE PRECISION RWORK( * ), W( * )
290 COMPLEX*16 AP( * ), BP( * ), WORK( * ), Z( LDZ, * )
291* ..
292*
293* =====================================================================
294*
295* .. Local Scalars ..
296 LOGICAL ALLEIG, INDEIG, UPPER, VALEIG, WANTZ
297 CHARACTER TRANS
298 INTEGER J
299* ..
300* .. External Functions ..
301 LOGICAL LSAME
302 EXTERNAL lsame
303* ..
304* .. External Subroutines ..
305 EXTERNAL xerbla, zhpevx, zhpgst, zpptrf, ztpmv, ztpsv
306* ..
307* .. Intrinsic Functions ..
308 INTRINSIC min
309* ..
310* .. Executable Statements ..
311*
312* Test the input parameters.
313*
314 wantz = lsame( jobz, 'V' )
315 upper = lsame( uplo, 'U' )
316 alleig = lsame( range, 'A' )
317 valeig = lsame( range, 'V' )
318 indeig = lsame( range, 'I' )
319*
320 info = 0
321 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
322 info = -1
323 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
324 info = -2
325 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
326 info = -3
327 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
328 info = -4
329 ELSE IF( n.LT.0 ) THEN
330 info = -5
331 ELSE
332 IF( valeig ) THEN
333 IF( n.GT.0 .AND. vu.LE.vl ) THEN
334 info = -9
335 END IF
336 ELSE IF( indeig ) THEN
337 IF( il.LT.1 ) THEN
338 info = -10
339 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
340 info = -11
341 END IF
342 END IF
343 END IF
344 IF( info.EQ.0 ) THEN
345 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
346 info = -16
347 END IF
348 END IF
349*
350 IF( info.NE.0 ) THEN
351 CALL xerbla( 'ZHPGVX', -info )
352 RETURN
353 END IF
354*
355* Quick return if possible
356*
357 IF( n.EQ.0 )
358 $ RETURN
359*
360* Form a Cholesky factorization of B.
361*
362 CALL zpptrf( uplo, n, bp, info )
363 IF( info.NE.0 ) THEN
364 info = n + info
365 RETURN
366 END IF
367*
368* Transform problem to standard eigenvalue problem and solve.
369*
370 CALL zhpgst( itype, uplo, n, ap, bp, info )
371 CALL zhpevx( jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m,
372 $ w, z, ldz, work, rwork, iwork, ifail, info )
373*
374 IF( wantz ) THEN
375*
376* Backtransform eigenvectors to the original problem.
377*
378 IF( info.GT.0 )
379 $ m = info - 1
380 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
381*
382* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
383* backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
384*
385 IF( upper ) THEN
386 trans = 'N'
387 ELSE
388 trans = 'C'
389 END IF
390*
391 DO 10 j = 1, m
392 CALL ztpsv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
393 $ 1 )
394 10 CONTINUE
395*
396 ELSE IF( itype.EQ.3 ) THEN
397*
398* For B*A*x=(lambda)*x;
399* backtransform eigenvectors: x = L*y or U**H *y
400*
401 IF( upper ) THEN
402 trans = 'C'
403 ELSE
404 trans = 'N'
405 END IF
406*
407 DO 20 j = 1, m
408 CALL ztpmv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
409 $ 1 )
410 20 CONTINUE
411 END IF
412 END IF
413*
414 RETURN
415*
416* End of ZHPGVX
417*
zhpevx
subroutine zhpevx(jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
ZHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition zhpevx.f:240