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

Functions

subroutine cggglm (n, m, p, a, lda, b, ldb, d, x, y, work, lwork, info)
 CGGGLM
subroutine chbev (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, rwork, info)
  CHBEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chbev_2stage (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, rwork, info)
  CHBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chbevd (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
  CHBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chbevd_2stage (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
  CHBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chbevx (jobz, range, uplo, n, kd, ab, ldab, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
  CHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chbevx_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)
  CHBEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chbgv (jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z, ldz, work, rwork, info)
 CHBGV
subroutine chbgvd (jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
 CHBGVD
subroutine chbgvx (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)
 CHBGVX
subroutine chpev (jobz, uplo, n, ap, w, z, ldz, work, rwork, info)
  CHPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chpevd (jobz, uplo, n, ap, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
  CHPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chpevx (jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
  CHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine chpgv (itype, jobz, uplo, n, ap, bp, w, z, ldz, work, rwork, info)
 CHPGV
subroutine chpgvd (itype, jobz, uplo, n, ap, bp, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
 CHPGVD
subroutine chpgvx (itype, jobz, range, uplo, n, ap, bp, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
 CHPGVX

Detailed Description

This is the group of complex Other Eigenvalue routines

Function Documentation

◆ cggglm()

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

CGGGLM

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

Purpose:
!>
!> CGGGLM 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 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 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 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 array, dimension (M)
!>
[out]Y
!>          Y is COMPLEX array, dimension (P)
!>
!>          On exit, X and Y are the solutions of the GLM problem.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. 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
!>          CGEQRF, CGERQF, CUNMQR and CUNMRQ.
!>
!>          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 cggglm.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 A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195 $ X( * ), Y( * )
196* ..
197*
198* ===================================================================
199*
200* .. Parameters ..
201 COMPLEX CZERO, CONE
202 parameter( czero = ( 0.0e+0, 0.0e+0 ),
203 $ cone = ( 1.0e+0, 0.0e+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 ccopy, cgemv, cggqrf, ctrtrs, cunmqr, cunmrq,
212 $ xerbla
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, 'CGEQRF', ' ', n, m, -1, -1 )
248 nb2 = ilaenv( 1, 'CGERQF', ' ', n, m, -1, -1 )
249 nb3 = ilaenv( 1, 'CUNMQR', ' ', n, m, p, -1 )
250 nb4 = ilaenv( 1, 'CUNMRQ', ' ', 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( 'CGGGLM', -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 cggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
291 $ work( m+np+1 ), lwork-m-np, info )
292 lopt = real( work( m+np+1 ) )
293*
294* Update left-hand-side vector d = Q**H*d = ( d1 ) M
295* ( d2 ) N-M
296*
297 CALL cunmqr( '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 ctrtrs( '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 ccopy( 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 cgemv( '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 ctrtrs( '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 ccopy( m, d, 1, x, 1 )
340 END IF
341*
342* Backward transformation y = Z**H *y
343*
344 CALL cunmrq( '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 CGGGLM
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
cunmrq
subroutine cunmrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMRQ
Definition cunmrq.f:168
cunmqr
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
cggqrf
subroutine cggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
CGGQRF
Definition cggqrf.f:215
ctrtrs
subroutine ctrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
CTRTRS
Definition ctrtrs.f:140
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ chbev()

subroutine chbev ( character jobz,
character uplo,
integer n,
integer kd,
complex, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) w,
complex, dimension( ldz, * ) z,
integer ldz,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> CHBEV 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 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 array, dimension (N)
!>
[out]RWORK
!>          RWORK is REAL array, dimension (max(1,3*N-2))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 150 of file chbev.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 REAL RWORK( * ), W( * )
163 COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * )
164* ..
165*
166* =====================================================================
167*
168* .. Parameters ..
169 REAL ZERO, ONE
170 parameter( zero = 0.0e0, one = 1.0e0 )
171* ..
172* .. Local Scalars ..
173 LOGICAL LOWER, WANTZ
174 INTEGER IINFO, IMAX, INDE, INDRWK, ISCALE
175 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
176 $ SMLNUM
177* ..
178* .. External Functions ..
179 LOGICAL LSAME
180 REAL CLANHB, SLAMCH
181 EXTERNAL lsame, clanhb, slamch
182* ..
183* .. External Subroutines ..
184 EXTERNAL chbtrd, clascl, csteqr, sscal, ssterf, xerbla
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( 'CHBEV ', -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 ) = real( ab( 1, 1 ) )
224 ELSE
225 w( 1 ) = real( 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 = slamch( 'Safe minimum' )
235 eps = slamch( '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 = clanhb( '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 clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
255 ELSE
256 CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
257 END IF
258 END IF
259*
260* Call CHBTRD to reduce Hermitian band matrix to tridiagonal form.
261*
262 inde = 1
263 CALL chbtrd( jobz, uplo, n, kd, ab, ldab, w, rwork( inde ), z,
264 $ ldz, work, iinfo )
265*
266* For eigenvalues only, call SSTERF. For eigenvectors, call CSTEQR.
267*
268 IF( .NOT.wantz ) THEN
269 CALL ssterf( n, w, rwork( inde ), info )
270 ELSE
271 indrwk = inde + n
272 CALL csteqr( 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 sscal( imax, one / sigma, w, 1 )
285 END IF
286*
287 RETURN
288*
289* End of CHBEV
290*
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
clascl
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
clanhb
real function clanhb(norm, uplo, n, k, ab, ldab, work)
CLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhb.f:132
chbtrd
subroutine chbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
CHBTRD
Definition chbtrd.f:163
csteqr
subroutine csteqr(compz, n, d, e, z, ldz, work, info)
CSTEQR
Definition csteqr.f:132
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68

◆ chbev_2stage()

subroutine chbev_2stage ( character jobz,
character uplo,
integer n,
integer kd,
complex, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) w,
complex, dimension( ldz, * ) z,
integer ldz,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> CHBEV_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 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 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 REAL array, dimension (max(1,3*N-2))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  All details about the 2stage techniques are available in:
!>
!>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
!>  Parallel reduction to condensed forms for symmetric eigenvalue problems
!>  using aggregated fine-grained and memory-aware kernels. In Proceedings
!>  of 2011 International Conference for High Performance Computing,
!>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
!>  Article 8 , 11 pages.
!>  http://doi.acm.org/10.1145/2063384.2063394
!>
!>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
!>  An improved parallel singular value algorithm and its implementation 
!>  for multicore hardware, In Proceedings of 2013 International Conference
!>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
!>  Denver, Colorado, USA, 2013.
!>  Article 90, 12 pages.
!>  http://doi.acm.org/10.1145/2503210.2503292
!>
!>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
!>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
!>  calculations based on fine-grained memory aware tasks.
!>  International Journal of High Performance Computing Applications.
!>  Volume 28 Issue 2, Pages 196-209, May 2014.
!>  http://hpc.sagepub.com/content/28/2/196 
!>
!>

Definition at line 209 of file chbev_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 REAL RWORK( * ), W( * )
224 COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * )
225* ..
226*
227* =====================================================================
228*
229* .. Parameters ..
230 REAL ZERO, ONE
231 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
238 $ SMLNUM
239* ..
240* .. External Functions ..
241 LOGICAL LSAME
242 INTEGER ILAENV2STAGE
243 REAL SLAMCH, CLANHB
244 EXTERNAL lsame, slamch, clanhb, ilaenv2stage
245* ..
246* .. External Subroutines ..
247 EXTERNAL sscal, ssterf, xerbla, clascl, csteqr,
248 $ chetrd_2stage, chetrd_hb2st
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC real, 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, 'CHETRD_HB2ST', jobz,
282 $ n, kd, -1, -1 )
283 lhtrd = ilaenv2stage( 3, 'CHETRD_HB2ST', jobz,
284 $ n, kd, ib, -1 )
285 lwtrd = ilaenv2stage( 4, 'CHETRD_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( 'CHBEV_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 ) = real( ab( 1, 1 ) )
310 ELSE
311 w( 1 ) = real( 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 = slamch( 'Safe minimum' )
321 eps = slamch( '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 = clanhb( '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 clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
341 ELSE
342 CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
343 END IF
344 END IF
345*
346* Call CHBTRD_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 chetrd_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 SSTERF. For eigenvectors, call CSTEQR.
358*
359 IF( .NOT.wantz ) THEN
360 CALL ssterf( n, w, rwork( inde ), info )
361 ELSE
362 indrwk = inde + n
363 CALL csteqr( 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 sscal( 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 CHBEV_2STAGE
385*
ilaenv2stage
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
Definition ilaenv2stage.f:149
chetrd_2stage
subroutine chetrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
CHETRD_2STAGE
Definition chetrd_2stage.f:224
chetrd_hb2st
subroutine chetrd_hb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
CHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T
Definition chetrd_hb2st.F:230

◆ chbevd()

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

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

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

Purpose:
!>
!> CHBEVD 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 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 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 REAL 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 chbevd.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 REAL RWORK( * ), W( * )
227 COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * )
228* ..
229*
230* =====================================================================
231*
232* .. Parameters ..
233 REAL ZERO, ONE
234 parameter( zero = 0.0e0, one = 1.0e0 )
235 COMPLEX CZERO, CONE
236 parameter( czero = ( 0.0e0, 0.0e0 ),
237 $ cone = ( 1.0e0, 0.0e0 ) )
238* ..
239* .. Local Scalars ..
240 LOGICAL LOWER, LQUERY, WANTZ
241 INTEGER IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE,
242 $ LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN
243 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
244 $ SMLNUM
245* ..
246* .. External Functions ..
247 LOGICAL LSAME
248 REAL CLANHB, SLAMCH
249 EXTERNAL lsame, clanhb, slamch
250* ..
251* .. External Subroutines ..
252 EXTERNAL cgemm, chbtrd, clacpy, clascl, cstedc, sscal,
253 $ ssterf, xerbla
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( 'CHBEVD', -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 ) = real( ab( 1, 1 ) )
324 IF( wantz )
325 $ z( 1, 1 ) = cone
326 RETURN
327 END IF
328*
329* Get machine constants.
330*
331 safmin = slamch( 'Safe minimum' )
332 eps = slamch( '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 = clanhb( '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 clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
352 ELSE
353 CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
354 END IF
355 END IF
356*
357* Call CHBTRD 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 chbtrd( jobz, uplo, n, kd, ab, ldab, w, rwork( inde ), z,
365 $ ldz, work, iinfo )
366*
367* For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC.
368*
369 IF( .NOT.wantz ) THEN
370 CALL ssterf( n, w, rwork( inde ), info )
371 ELSE
372 CALL cstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
373 $ llwk2, rwork( indwrk ), llrwk, iwork, liwork,
374 $ info )
375 CALL cgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
376 $ work( indwk2 ), n )
377 CALL clacpy( '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 sscal( 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 CHBEVD
397*
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
cstedc
subroutine cstedc(compz, n, d, e, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
CSTEDC
Definition cstedc.f:212
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187

◆ chbevd_2stage()

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

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

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

Purpose:
!>
!> CHBEVD_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 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 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 REAL 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 chbevd_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 REAL RWORK( * ), W( * )
274 COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * )
275* ..
276*
277* =====================================================================
278*
279* .. Parameters ..
280 REAL ZERO, ONE
281 parameter( zero = 0.0e0, one = 1.0e0 )
282 COMPLEX CZERO, CONE
283 parameter( czero = ( 0.0e0, 0.0e0 ),
284 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
292 $ SMLNUM
293* ..
294* .. External Functions ..
295 LOGICAL LSAME
296 INTEGER ILAENV2STAGE
297 REAL SLAMCH, CLANHB
298 EXTERNAL lsame, slamch, clanhb, ilaenv2stage
299* ..
300* .. External Subroutines ..
301 EXTERNAL sscal, ssterf, xerbla, cgemm, clacpy,
302 $ clascl, cstedc, chetrd_hb2st
303* ..
304* .. Intrinsic Functions ..
305 INTRINSIC real, 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, 'CHETRD_HB2ST', jobz, n, kd, -1, -1 )
322 lhtrd = ilaenv2stage( 3, 'CHETRD_HB2ST', jobz, n, kd, ib, -1 )
323 lwtrd = ilaenv2stage( 4, 'CHETRD_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( 'CHBEVD_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 ) = real( ab( 1, 1 ) )
376 IF( wantz )
377 $ z( 1, 1 ) = cone
378 RETURN
379 END IF
380*
381* Get machine constants.
382*
383 safmin = slamch( 'Safe minimum' )
384 eps = slamch( '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 = clanhb( '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 clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
404 ELSE
405 CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
406 END IF
407 END IF
408*
409* Call CHBTRD_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 chetrd_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 SSTERF. For eigenvectors, call CSTEDC.
425*
426 IF( .NOT.wantz ) THEN
427 CALL ssterf( n, w, rwork( inde ), info )
428 ELSE
429 CALL cstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
430 $ llwk2, rwork( indrwk ), llrwk, iwork, liwork,
431 $ info )
432 CALL cgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
433 $ work( indwk2 ), n )
434 CALL clacpy( '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 sscal( 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 CHBEVD_2STAGE
454*

◆ chbevx()

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

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

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

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

Definition at line 264 of file chbevx.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 REAL ABSTOL, VL, VU
276* ..
277* .. Array Arguments ..
278 INTEGER IFAIL( * ), IWORK( * )
279 REAL RWORK( * ), W( * )
280 COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ),
281 $ Z( LDZ, * )
282* ..
283*
284* =====================================================================
285*
286* .. Parameters ..
287 REAL ZERO, ONE
288 parameter( zero = 0.0e0, one = 1.0e0 )
289 COMPLEX CZERO, CONE
290 parameter( czero = ( 0.0e0, 0.0e0 ),
291 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
300 $ SIGMA, SMLNUM, TMP1, VLL, VUU
301 COMPLEX CTMP1
302* ..
303* .. External Functions ..
304 LOGICAL LSAME
305 REAL CLANHB, SLAMCH
306 EXTERNAL lsame, clanhb, slamch
307* ..
308* .. External Subroutines ..
309 EXTERNAL ccopy, cgemv, chbtrd, clacpy, clascl, cstein,
310 $ csteqr, cswap, scopy, sscal, sstebz, ssterf,
311 $ xerbla
312* ..
313* .. Intrinsic Functions ..
314 INTRINSIC max, min, real, 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( 'CHBEVX', -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 = real( 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 ) = real( 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 = slamch( 'Safe minimum' )
392 eps = slamch( '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 ENDIF
409 anrm = clanhb( '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 clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
420 ELSE
421 CALL clascl( '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 CHBTRD 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 chbtrd( 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 SSTERF or CSTEQR. If this fails for some
442* eigenvalue, then try SSTEBZ.
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 scopy( n, rwork( indd ), 1, w, 1 )
452 indee = indrwk + 2*n
453 IF( .NOT.wantz ) THEN
454 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
455 CALL ssterf( n, w, rwork( indee ), info )
456 ELSE
457 CALL clacpy( 'A', n, n, q, ldq, z, ldz )
458 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
459 CALL csteqr( 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 SSTEBZ and, if eigenvectors are desired, CSTEIN.
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 sstebz( 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 cstein( 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 CSTEIN.
496*
497 DO 20 j = 1, m
498 CALL ccopy( n, z( 1, j ), 1, work( 1 ), 1 )
499 CALL cgemv( '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 sscal( 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 cswap( 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 CHBEVX
549*
sstebz
subroutine sstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
SSTEBZ
Definition sstebz.f:273
cstein
subroutine cstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
CSTEIN
Definition cstein.f:182
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82

◆ chbevx_2stage()

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

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

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

Purpose:
!>
!> CHBEVX_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 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 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 REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AB to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is COMPLEX array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (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 REAL array, dimension (7*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, then i eigenvectors failed to converge.
!>                Their indices are stored in array IFAIL.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  All details about the 2stage techniques are available in:
!>
!>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
!>  Parallel reduction to condensed forms for symmetric eigenvalue problems
!>  using aggregated fine-grained and memory-aware kernels. In Proceedings
!>  of 2011 International Conference for High Performance Computing,
!>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
!>  Article 8 , 11 pages.
!>  http://doi.acm.org/10.1145/2063384.2063394
!>
!>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
!>  An improved parallel singular value algorithm and its implementation 
!>  for multicore hardware, In Proceedings of 2013 International Conference
!>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
!>  Denver, Colorado, USA, 2013.
!>  Article 90, 12 pages.
!>  http://doi.acm.org/10.1145/2503210.2503292
!>
!>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
!>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
!>  calculations based on fine-grained memory aware tasks.
!>  International Journal of High Performance Computing Applications.
!>  Volume 28 Issue 2, Pages 196-209, May 2014.
!>  http://hpc.sagepub.com/content/28/2/196 
!>
!>

Definition at line 323 of file chbevx_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 REAL ABSTOL, VL, VU
338* ..
339* .. Array Arguments ..
340 INTEGER IFAIL( * ), IWORK( * )
341 REAL RWORK( * ), W( * )
342 COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ),
343 $ Z( LDZ, * )
344* ..
345*
346* =====================================================================
347*
348* .. Parameters ..
349 REAL ZERO, ONE
350 parameter( zero = 0.0e0, one = 1.0e0 )
351 COMPLEX CZERO, CONE
352 parameter( czero = ( 0.0e0, 0.0e0 ),
353 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
364 $ SIGMA, SMLNUM, TMP1, VLL, VUU
365 COMPLEX CTMP1
366* ..
367* .. External Functions ..
368 LOGICAL LSAME
369 INTEGER ILAENV2STAGE
370 REAL SLAMCH, CLANHB
371 EXTERNAL lsame, slamch, clanhb, ilaenv2stage
372* ..
373* .. External Subroutines ..
374 EXTERNAL scopy, sscal, sstebz, ssterf, xerbla, ccopy,
375 $ cgemv, clacpy, clascl, cstein, csteqr,
376 $ cswap, chetrd_hb2st
377* ..
378* .. Intrinsic Functions ..
379 INTRINSIC real, 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, 'CHETRD_HB2ST', jobz,
430 $ n, kd, -1, -1 )
431 lhtrd = ilaenv2stage( 3, 'CHETRD_HB2ST', jobz,
432 $ n, kd, ib, -1 )
433 lwtrd = ilaenv2stage( 4, 'CHETRD_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( 'CHBEVX_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 = real( 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 ) = real( 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 = slamch( 'Safe minimum' )
479 eps = slamch( '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 = clanhb( '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 clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
507 ELSE
508 CALL clascl( '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 CHBTRD_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 chetrd_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 SSTERF or CSTEQR. If this fails for some
534* eigenvalue, then try SSTEBZ.
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 scopy( n, rwork( indd ), 1, w, 1 )
544 indee = indrwk + 2*n
545 IF( .NOT.wantz ) THEN
546 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
547 CALL ssterf( n, w, rwork( indee ), info )
548 ELSE
549 CALL clacpy( 'A', n, n, q, ldq, z, ldz )
550 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
551 CALL csteqr( 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 SSTEBZ and, if eigenvectors are desired, CSTEIN.
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 sstebz( 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 cstein( 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 CSTEIN.
588*
589 DO 20 j = 1, m
590 CALL ccopy( n, z( 1, j ), 1, work( 1 ), 1 )
591 CALL cgemv( '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 sscal( 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 cswap( 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 CHBEVX_2STAGE
645*

◆ chbgv()

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

CHBGV

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

Purpose:
!>
!> CHBGV 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 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 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 CPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 array, dimension (N)
!>
[out]RWORK
!>          RWORK is REAL 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 CPBSTF
!>                    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 chbgv.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 REAL RWORK( * ), W( * )
194 COMPLEX 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 chbgst, chbtrd, cpbstf, csteqr, ssterf, xerbla
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( 'CHBGV ', -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 cpbstf( 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 chbgst( 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 chbtrd( vect, uplo, n, ka, ab, ldab, w, rwork( inde ), z,
270 $ ldz, work, iinfo )
271*
272* For eigenvalues only, call SSTERF. For eigenvectors, call CSTEQR.
273*
274 IF( .NOT.wantz ) THEN
275 CALL ssterf( n, w, rwork( inde ), info )
276 ELSE
277 CALL csteqr( jobz, n, w, rwork( inde ), z, ldz,
278 $ rwork( indwrk ), info )
279 END IF
280 RETURN
281*
282* End of CHBGV
283*
chbgst
subroutine chbgst(vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x, ldx, work, rwork, info)
CHBGST
Definition chbgst.f:165
cpbstf
subroutine cpbstf(uplo, n, kd, ab, ldab, info)
CPBSTF
Definition cpbstf.f:153

◆ chbgvd()

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

CHBGVD

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

Purpose:
!>
!> CHBGVD 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 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 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 CPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 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 REAL array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO=0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of 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 CPBSTF
!>                    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 chbgvd.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 REAL RWORK( * ), W( * )
265 COMPLEX AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
266 $ Z( LDZ, * )
267* ..
268*
269* =====================================================================
270*
271* .. Parameters ..
272 COMPLEX CONE, CZERO
273 parameter( cone = ( 1.0e+0, 0.0e+0 ),
274 $ czero = ( 0.0e+0, 0.0e+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 ssterf, xerbla, cgemm, chbgst, chbtrd, clacpy,
288 $ cpbstf, cstedc
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( 'CHBGVD', -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 cpbstf( 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 chbgst( 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 chbtrd( vect, uplo, n, ka, ab, ldab, w, rwork( inde ), z,
382 $ ldz, work, iinfo )
383*
384* For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC.
385*
386 IF( .NOT.wantz ) THEN
387 CALL ssterf( n, w, rwork( inde ), info )
388 ELSE
389 CALL cstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
390 $ llwk2, rwork( indwrk ), llrwk, iwork, liwork,
391 $ info )
392 CALL cgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
393 $ work( indwk2 ), n )
394 CALL clacpy( '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 CHBGVD
403*

◆ chbgvx()

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

CHBGVX

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

Purpose:
!>
!> CHBGVX 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 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 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 CPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]Q
!>          Q is COMPLEX 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 REAL
!>
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AP to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('S').
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 array, dimension (N)
!>
[out]RWORK
!>          RWORK is REAL array, dimension (7*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, 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 CPBSTF
!>                    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 chbgvx.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 REAL ABSTOL, VL, VU
310* ..
311* .. Array Arguments ..
312 INTEGER IFAIL( * ), IWORK( * )
313 REAL RWORK( * ), W( * )
314 COMPLEX AB( LDAB, * ), BB( LDBB, * ), Q( LDQ, * ),
315 $ WORK( * ), Z( LDZ, * )
316* ..
317*
318* =====================================================================
319*
320* .. Parameters ..
321 REAL ZERO
322 parameter( zero = 0.0e+0 )
323 COMPLEX CZERO, CONE
324 parameter( czero = ( 0.0e+0, 0.0e+0 ),
325 $ cone = ( 1.0e+0, 0.0e+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 REAL TMP1
333* ..
334* .. External Functions ..
335 LOGICAL LSAME
336 EXTERNAL lsame
337* ..
338* .. External Subroutines ..
339 EXTERNAL ccopy, cgemv, chbgst, chbtrd, clacpy, cpbstf,
340 $ cstein, csteqr, cswap, scopy, sstebz, ssterf,
341 $ xerbla
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( 'CHBGVX', -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 cpbstf( 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 chbgst( 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 chbtrd( 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 SSTERF or CSTEQR. If this fails for some
434* eigenvalue, then try SSTEBZ.
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 scopy( n, rwork( indd ), 1, w, 1 )
444 indee = indrwk + 2*n
445 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
446 IF( .NOT.wantz ) THEN
447 CALL ssterf( n, w, rwork( indee ), info )
448 ELSE
449 CALL clacpy( 'A', n, n, q, ldq, z, ldz )
450 CALL csteqr( 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 SSTEBZ and, if eigenvectors are desired,
466* call CSTEIN.
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 sstebz( 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 cstein( 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 CSTEIN.
488*
489 DO 20 j = 1, m
490 CALL ccopy( n, z( 1, j ), 1, work( 1 ), 1 )
491 CALL cgemv( '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 cswap( 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 CHBGVX
531*

◆ chpev()

subroutine chpev ( character jobz,
character uplo,
integer n,
complex, dimension( * ) ap,
real, dimension( * ) w,
complex, dimension( ldz, * ) z,
integer ldz,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> CHPEV 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 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 array, dimension (max(1, 2*N-1))
!>
[out]RWORK
!>          RWORK is REAL array, dimension (max(1, 3*N-2))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 136 of file chpev.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 REAL RWORK( * ), W( * )
149 COMPLEX AP( * ), WORK( * ), Z( LDZ, * )
150* ..
151*
152* =====================================================================
153*
154* .. Parameters ..
155 REAL ZERO, ONE
156 parameter( zero = 0.0e0, one = 1.0e0 )
157* ..
158* .. Local Scalars ..
159 LOGICAL WANTZ
160 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWRK,
161 $ ISCALE
162 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
163 $ SMLNUM
164* ..
165* .. External Functions ..
166 LOGICAL LSAME
167 REAL CLANHP, SLAMCH
168 EXTERNAL lsame, clanhp, slamch
169* ..
170* .. External Subroutines ..
171 EXTERNAL chptrd, csscal, csteqr, cupgtr, sscal, ssterf,
172 $ xerbla
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( 'CHPEV ', -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 ) = real( 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 = slamch( 'Safe minimum' )
216 eps = slamch( '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 = clanhp( '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 csscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
235 END IF
236*
237* Call CHPTRD to reduce Hermitian packed matrix to tridiagonal form.
238*
239 inde = 1
240 indtau = 1
241 CALL chptrd( uplo, n, ap, w, rwork( inde ), work( indtau ),
242 $ iinfo )
243*
244* For eigenvalues only, call SSTERF. For eigenvectors, first call
245* CUPGTR to generate the orthogonal matrix, then call CSTEQR.
246*
247 IF( .NOT.wantz ) THEN
248 CALL ssterf( n, w, rwork( inde ), info )
249 ELSE
250 indwrk = indtau + n
251 CALL cupgtr( uplo, n, ap, work( indtau ), z, ldz,
252 $ work( indwrk ), iinfo )
253 indrwk = inde + n
254 CALL csteqr( 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 sscal( imax, one / sigma, w, 1 )
267 END IF
268*
269 RETURN
270*
271* End of CHPEV
272*
clanhp
real function clanhp(norm, uplo, n, ap, work)
CLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhp.f:117
cupgtr
subroutine cupgtr(uplo, n, ap, tau, q, ldq, work, info)
CUPGTR
Definition cupgtr.f:114
chptrd
subroutine chptrd(uplo, n, ap, d, e, tau, info)
CHPTRD
Definition chptrd.f:151
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78

◆ chpevd()

subroutine chpevd ( character jobz,
character uplo,
integer n,
complex, dimension( * ) ap,
real, dimension( * ) w,
complex, dimension( ldz, * ) z,
integer ldz,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> CHPEVD 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 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 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 REAL 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 chpevd.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 REAL RWORK( * ), W( * )
212 COMPLEX AP( * ), WORK( * ), Z( LDZ, * )
213* ..
214*
215* =====================================================================
216*
217* .. Parameters ..
218 REAL ZERO, ONE
219 parameter( zero = 0.0e+0, one = 1.0e+0 )
220 COMPLEX CONE
221 parameter( cone = ( 1.0e+0, 0.0e+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 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
228 $ SMLNUM
229* ..
230* .. External Functions ..
231 LOGICAL LSAME
232 REAL CLANHP, SLAMCH
233 EXTERNAL lsame, clanhp, slamch
234* ..
235* .. External Subroutines ..
236 EXTERNAL chptrd, csscal, cstedc, cupmtr, sscal, ssterf,
237 $ xerbla
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( 'CHPEVD', -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 ) = real( ap( 1 ) )
304 IF( wantz )
305 $ z( 1, 1 ) = cone
306 RETURN
307 END IF
308*
309* Get machine constants.
310*
311 safmin = slamch( 'Safe minimum' )
312 eps = slamch( '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 = clanhp( '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 csscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
331 END IF
332*
333* Call CHPTRD 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 chptrd( uplo, n, ap, w, rwork( inde ), work( indtau ),
342 $ iinfo )
343*
344* For eigenvalues only, call SSTERF. For eigenvectors, first call
345* CUPGTR to generate the orthogonal matrix, then call CSTEDC.
346*
347 IF( .NOT.wantz ) THEN
348 CALL ssterf( n, w, rwork( inde ), info )
349 ELSE
350 CALL cstedc( 'I', n, w, rwork( inde ), z, ldz, work( indwrk ),
351 $ llwrk, rwork( indrwk ), llrwk, iwork, liwork,
352 $ info )
353 CALL cupmtr( '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 sscal( 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 CHPEVD
374*
cupmtr
subroutine cupmtr(side, uplo, trans, m, n, ap, tau, c, ldc, work, info)
CUPMTR
Definition cupmtr.f:150

◆ chpevx()

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

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

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

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

Definition at line 237 of file chpevx.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 REAL ABSTOL, VL, VU
249* ..
250* .. Array Arguments ..
251 INTEGER IFAIL( * ), IWORK( * )
252 REAL RWORK( * ), W( * )
253 COMPLEX AP( * ), WORK( * ), Z( LDZ, * )
254* ..
255*
256* =====================================================================
257*
258* .. Parameters ..
259 REAL ZERO, ONE
260 parameter( zero = 0.0e0, one = 1.0e0 )
261 COMPLEX CONE
262 parameter( cone = ( 1.0e0, 0.0e0 ) )
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 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
271 $ SIGMA, SMLNUM, TMP1, VLL, VUU
272* ..
273* .. External Functions ..
274 LOGICAL LSAME
275 REAL CLANHP, SLAMCH
276 EXTERNAL lsame, clanhp, slamch
277* ..
278* .. External Subroutines ..
279 EXTERNAL chptrd, csscal, cstein, csteqr, cswap, cupgtr,
280 $ cupmtr, scopy, sscal, sstebz, ssterf, xerbla
281* ..
282* .. Intrinsic Functions ..
283 INTRINSIC max, min, real, 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( 'CHPEVX', -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 ) = real( ap( 1 ) )
336 ELSE
337 IF( vl.LT.real( ap( 1 ) ) .AND. vu.GE.real( ap( 1 ) ) ) THEN
338 m = 1
339 w( 1 ) = real( 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 = slamch( 'Safe minimum' )
350 eps = slamch( '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 ENDIF
367 anrm = clanhp( '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 csscal( ( 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 CHPTRD 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 chptrd( 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 SSTERF or CUPGTR and CSTEQR. If this fails
397* for some eigenvalue, then try SSTEBZ.
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 scopy( n, rwork( indd ), 1, w, 1 )
407 indee = indrwk + 2*n
408 IF( .NOT.wantz ) THEN
409 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
410 CALL ssterf( n, w, rwork( indee ), info )
411 ELSE
412 CALL cupgtr( uplo, n, ap, work( indtau ), z, ldz,
413 $ work( indwrk ), iinfo )
414 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
415 CALL csteqr( 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 SSTEBZ and, if eigenvectors are desired, CSTEIN.
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 sstebz( 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 cstein( 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 CSTEIN.
452*
453 indwrk = indtau + n
454 CALL cupmtr( '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 sscal( 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 cswap( 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 CHPEVX
503*

◆ chpgv()

subroutine chpgv ( integer itype,
character jobz,
character uplo,
integer n,
complex, dimension( * ) ap,
complex, dimension( * ) bp,
real, dimension( * ) w,
complex, dimension( ldz, * ) z,
integer ldz,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

CHPGV

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

Purpose:
!>
!> CHPGV 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 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 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 array, dimension (max(1, 2*N-1))
!>
[out]RWORK
!>          RWORK is REAL array, dimension (max(1, 3*N-2))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  CPPTRF or CHPEV returned an error code:
!>             <= N:  if INFO = i, CHPEV 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 chpgv.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 REAL RWORK( * ), W( * )
176 COMPLEX 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 chpev, chpgst, cpptrf, ctpmv, ctpsv, xerbla
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( 'CHPGV ', -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 cpptrf( 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 chpgst( itype, uplo, n, ap, bp, info )
233 CALL chpev( 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 ctpsv( 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 ctpmv( 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 CHPGV
278*
cpptrf
subroutine cpptrf(uplo, n, ap, info)
CPPTRF
Definition cpptrf.f:119
chpgst
subroutine chpgst(itype, uplo, n, ap, bp, info)
CHPGST
Definition chpgst.f:113
chpev
subroutine chpev(jobz, uplo, n, ap, w, z, ldz, work, rwork, info)
CHPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
Definition chpev.f:138
ctpmv
subroutine ctpmv(uplo, trans, diag, n, ap, x, incx)
CTPMV
Definition ctpmv.f:142
ctpsv
subroutine ctpsv(uplo, trans, diag, n, ap, x, incx)
CTPSV
Definition ctpsv.f:144

◆ chpgvd()

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

CHPGVD

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

Purpose:
!>
!> CHPGVD 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 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 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 REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX 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 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 >= 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 REAL 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:  CPPTRF or CHPEVD returned an error code:
!>             <= N:  if INFO = i, CHPEVD 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 chpgvd.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 REAL RWORK( * ), W( * )
243 COMPLEX 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 chpevd, chpgst, cpptrf, ctpmv, ctpsv, xerbla
259* ..
260* .. Intrinsic Functions ..
261 INTRINSIC max, real
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( 'CHPGVD', -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 cpptrf( 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 chpgst( itype, uplo, n, ap, bp, info )
336 CALL chpevd( jobz, uplo, n, ap, w, z, ldz, work, lwork, rwork,
337 $ lrwork, iwork, liwork, info )
338 lwmin = max( real( lwmin ), real( work( 1 ) ) )
339 lrwmin = max( real( lrwmin ), real( rwork( 1 ) ) )
340 liwmin = max( real( liwmin ), real( 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 ctpsv( 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 ctpmv( 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 CHPGVD
389*
chpevd
subroutine chpevd(jobz, uplo, n, ap, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
CHPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition chpevd.f:200

◆ chpgvx()

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

CHPGVX

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

Purpose:
!>
!> CHPGVX 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 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 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 REAL
!>
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AP to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('S').
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX array, dimension (LDZ, 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 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL array, dimension (7*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  CPPTRF or CHPEVX returned an error code:
!>             <= N:  if INFO = i, CHPEVX 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 chpgvx.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 REAL ABSTOL, VL, VU
286* ..
287* .. Array Arguments ..
288 INTEGER IFAIL( * ), IWORK( * )
289 REAL RWORK( * ), W( * )
290 COMPLEX 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 chpevx, chpgst, cpptrf, ctpmv, ctpsv, xerbla
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( 'CHPGVX', -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 cpptrf( 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 chpgst( itype, uplo, n, ap, bp, info )
371 CALL chpevx( 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 ctpsv( 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 ctpmv( 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 CHPGVX
417*
chpevx
subroutine chpevx(jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m, w, z, ldz, work, rwork, iwork, ifail, info)
CHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition chpevx.f:240