LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zhpgvx.f
Go to the documentation of this file.
1 *> \brief \b ZHPGVX
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHPGVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpgvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpgvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpgvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHPGVX( ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, VL, VU,
22 * IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK,
23 * IWORK, IFAIL, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, ITYPE, IU, LDZ, M, N
28 * DOUBLE PRECISION ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * DOUBLE PRECISION RWORK( * ), W( * )
33 * COMPLEX*16 AP( * ), BP( * ), WORK( * ), Z( LDZ, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZHPGVX computes selected eigenvalues and, optionally, eigenvectors
43 *> of a complex generalized Hermitian-definite eigenproblem, of the form
44 *> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and
45 *> B are assumed to be Hermitian, stored in packed format, and B is also
46 *> positive definite. Eigenvalues and eigenvectors can be selected by
47 *> specifying either a range of values or a range of indices for the
48 *> desired eigenvalues.
49 *> \endverbatim
50 *
51 * Arguments:
52 * ==========
53 *
54 *> \param[in] ITYPE
55 *> \verbatim
56 *> ITYPE is INTEGER
57 *> Specifies the problem type to be solved:
58 *> = 1: A*x = (lambda)*B*x
59 *> = 2: A*B*x = (lambda)*x
60 *> = 3: B*A*x = (lambda)*x
61 *> \endverbatim
62 *>
63 *> \param[in] JOBZ
64 *> \verbatim
65 *> JOBZ is CHARACTER*1
66 *> = 'N': Compute eigenvalues only;
67 *> = 'V': Compute eigenvalues and eigenvectors.
68 *> \endverbatim
69 *>
70 *> \param[in] RANGE
71 *> \verbatim
72 *> RANGE is CHARACTER*1
73 *> = 'A': all eigenvalues will be found;
74 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
75 *> will be found;
76 *> = 'I': the IL-th through IU-th eigenvalues will be found.
77 *> \endverbatim
78 *>
79 *> \param[in] UPLO
80 *> \verbatim
81 *> UPLO is CHARACTER*1
82 *> = 'U': Upper triangles of A and B are stored;
83 *> = 'L': Lower triangles of A and B are stored.
84 *> \endverbatim
85 *>
86 *> \param[in] N
87 *> \verbatim
88 *> N is INTEGER
89 *> The order of the matrices A and B. N >= 0.
90 *> \endverbatim
91 *>
92 *> \param[in,out] AP
93 *> \verbatim
94 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
95 *> On entry, the upper or lower triangle of the Hermitian matrix
96 *> A, packed columnwise in a linear array. The j-th column of A
97 *> is stored in the array AP as follows:
98 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
99 *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
100 *>
101 *> On exit, the contents of AP are destroyed.
102 *> \endverbatim
103 *>
104 *> \param[in,out] BP
105 *> \verbatim
106 *> BP is COMPLEX*16 array, dimension (N*(N+1)/2)
107 *> On entry, the upper or lower triangle of the Hermitian matrix
108 *> B, packed columnwise in a linear array. The j-th column of B
109 *> is stored in the array BP as follows:
110 *> if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
111 *> if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
112 *>
113 *> On exit, the triangular factor U or L from the Cholesky
114 *> factorization B = U**H*U or B = L*L**H, in the same storage
115 *> format as B.
116 *> \endverbatim
117 *>
118 *> \param[in] VL
119 *> \verbatim
120 *> VL is DOUBLE PRECISION
121 *>
122 *> If RANGE='V', the lower bound of the interval to
123 *> be searched for eigenvalues. VL < VU.
124 *> Not referenced if RANGE = 'A' or 'I'.
125 *> \endverbatim
126 *>
127 *> \param[in] VU
128 *> \verbatim
129 *> VU is DOUBLE PRECISION
130 *>
131 *> If RANGE='V', the upper bound of the interval to
132 *> be searched for eigenvalues. VL < VU.
133 *> Not referenced if RANGE = 'A' or 'I'.
134 *> \endverbatim
135 *>
136 *> \param[in] IL
137 *> \verbatim
138 *> IL is INTEGER
139 *>
140 *> If RANGE='I', the index of the
141 *> smallest eigenvalue to be returned.
142 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
143 *> Not referenced if RANGE = 'A' or 'V'.
144 *> \endverbatim
145 *>
146 *> \param[in] IU
147 *> \verbatim
148 *> IU is INTEGER
149 *>
150 *> If RANGE='I', the index of the
151 *> largest eigenvalue to be returned.
152 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
153 *> Not referenced if RANGE = 'A' or 'V'.
154 *> \endverbatim
155 *>
156 *> \param[in] ABSTOL
157 *> \verbatim
158 *> ABSTOL is DOUBLE PRECISION
159 *> The absolute error tolerance for the eigenvalues.
160 *> An approximate eigenvalue is accepted as converged
161 *> when it is determined to lie in an interval [a,b]
162 *> of width less than or equal to
163 *>
164 *> ABSTOL + EPS * max( |a|,|b| ) ,
165 *>
166 *> where EPS is the machine precision. If ABSTOL is less than
167 *> or equal to zero, then EPS*|T| will be used in its place,
168 *> where |T| is the 1-norm of the tridiagonal matrix obtained
169 *> by reducing AP to tridiagonal form.
170 *>
171 *> Eigenvalues will be computed most accurately when ABSTOL is
172 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
173 *> If this routine returns with INFO>0, indicating that some
174 *> eigenvectors did not converge, try setting ABSTOL to
175 *> 2*DLAMCH('S').
176 *> \endverbatim
177 *>
178 *> \param[out] M
179 *> \verbatim
180 *> M is INTEGER
181 *> The total number of eigenvalues found. 0 <= M <= N.
182 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
183 *> \endverbatim
184 *>
185 *> \param[out] W
186 *> \verbatim
187 *> W is DOUBLE PRECISION array, dimension (N)
188 *> On normal exit, the first M elements contain the selected
189 *> eigenvalues in ascending order.
190 *> \endverbatim
191 *>
192 *> \param[out] Z
193 *> \verbatim
194 *> Z is COMPLEX*16 array, dimension (LDZ, N)
195 *> If JOBZ = 'N', then Z is not referenced.
196 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
197 *> contain the orthonormal eigenvectors of the matrix A
198 *> corresponding to the selected eigenvalues, with the i-th
199 *> column of Z holding the eigenvector associated with W(i).
200 *> The eigenvectors are normalized as follows:
201 *> if ITYPE = 1 or 2, Z**H*B*Z = I;
202 *> if ITYPE = 3, Z**H*inv(B)*Z = I.
203 *>
204 *> If an eigenvector fails to converge, then that column of Z
205 *> contains the latest approximation to the eigenvector, and the
206 *> index of the eigenvector is returned in IFAIL.
207 *> Note: the user must ensure that at least max(1,M) columns are
208 *> supplied in the array Z; if RANGE = 'V', the exact value of M
209 *> is not known in advance and an upper bound must be used.
210 *> \endverbatim
211 *>
212 *> \param[in] LDZ
213 *> \verbatim
214 *> LDZ is INTEGER
215 *> The leading dimension of the array Z. LDZ >= 1, and if
216 *> JOBZ = 'V', LDZ >= max(1,N).
217 *> \endverbatim
218 *>
219 *> \param[out] WORK
220 *> \verbatim
221 *> WORK is COMPLEX*16 array, dimension (2*N)
222 *> \endverbatim
223 *>
224 *> \param[out] RWORK
225 *> \verbatim
226 *> RWORK is DOUBLE PRECISION array, dimension (7*N)
227 *> \endverbatim
228 *>
229 *> \param[out] IWORK
230 *> \verbatim
231 *> IWORK is INTEGER array, dimension (5*N)
232 *> \endverbatim
233 *>
234 *> \param[out] IFAIL
235 *> \verbatim
236 *> IFAIL is INTEGER array, dimension (N)
237 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
238 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
239 *> indices of the eigenvectors that failed to converge.
240 *> If JOBZ = 'N', then IFAIL is not referenced.
241 *> \endverbatim
242 *>
243 *> \param[out] INFO
244 *> \verbatim
245 *> INFO is INTEGER
246 *> = 0: successful exit
247 *> < 0: if INFO = -i, the i-th argument had an illegal value
248 *> > 0: ZPPTRF or ZHPEVX returned an error code:
249 *> <= N: if INFO = i, ZHPEVX failed to converge;
250 *> i eigenvectors failed to converge. Their indices
251 *> are stored in array IFAIL.
252 *> > N: if INFO = N + i, for 1 <= i <= n, then the leading
253 *> minor of order i of B is not positive definite.
254 *> The factorization of B could not be completed and
255 *> no eigenvalues or eigenvectors were computed.
256 *> \endverbatim
257 *
258 * Authors:
259 * ========
260 *
261 *> \author Univ. of Tennessee
262 *> \author Univ. of California Berkeley
263 *> \author Univ. of Colorado Denver
264 *> \author NAG Ltd.
265 *
266 *> \date June 2016
267 *
268 *> \ingroup complex16OTHEReigen
269 *
270 *> \par Contributors:
271 * ==================
272 *>
273 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
274 *
275 * =====================================================================
276  SUBROUTINE zhpgvx( ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, VL, VU,
277  $ il, iu, abstol, m, w, z, ldz, work, rwork,
278  $ iwork, ifail, info )
279 *
280 * -- LAPACK driver routine (version 3.6.1) --
281 * -- LAPACK is a software package provided by Univ. of Tennessee, --
282 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
283 * June 2016
284 *
285 * .. Scalar Arguments ..
286  CHARACTER JOBZ, RANGE, UPLO
287  INTEGER IL, INFO, ITYPE, IU, LDZ, M, N
288  DOUBLE PRECISION ABSTOL, VL, VU
289 * ..
290 * .. Array Arguments ..
291  INTEGER IFAIL( * ), IWORK( * )
292  DOUBLE PRECISION RWORK( * ), W( * )
293  COMPLEX*16 AP( * ), BP( * ), WORK( * ), Z( ldz, * )
294 * ..
295 *
296 * =====================================================================
297 *
298 * .. Local Scalars ..
299  LOGICAL ALLEIG, INDEIG, UPPER, VALEIG, WANTZ
300  CHARACTER TRANS
301  INTEGER J
302 * ..
303 * .. External Functions ..
304  LOGICAL LSAME
305  EXTERNAL lsame
306 * ..
307 * .. External Subroutines ..
308  EXTERNAL xerbla, zhpevx, zhpgst, zpptrf, ztpmv, ztpsv
309 * ..
310 * .. Intrinsic Functions ..
311  INTRINSIC min
312 * ..
313 * .. Executable Statements ..
314 *
315 * Test the input parameters.
316 *
317  wantz = lsame( jobz, 'V' )
318  upper = lsame( uplo, 'U' )
319  alleig = lsame( range, 'A' )
320  valeig = lsame( range, 'V' )
321  indeig = lsame( range, 'I' )
322 *
323  info = 0
324  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
325  info = -1
326  ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
327  info = -2
328  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
329  info = -3
330  ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
331  info = -4
332  ELSE IF( n.LT.0 ) THEN
333  info = -5
334  ELSE
335  IF( valeig ) THEN
336  IF( n.GT.0 .AND. vu.LE.vl ) THEN
337  info = -9
338  END IF
339  ELSE IF( indeig ) THEN
340  IF( il.LT.1 ) THEN
341  info = -10
342  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
343  info = -11
344  END IF
345  END IF
346  END IF
347  IF( info.EQ.0 ) THEN
348  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
349  info = -16
350  END IF
351  END IF
352 *
353  IF( info.NE.0 ) THEN
354  CALL xerbla( 'ZHPGVX', -info )
355  RETURN
356  END IF
357 *
358 * Quick return if possible
359 *
360  IF( n.EQ.0 )
361  $ RETURN
362 *
363 * Form a Cholesky factorization of B.
364 *
365  CALL zpptrf( uplo, n, bp, info )
366  IF( info.NE.0 ) THEN
367  info = n + info
368  RETURN
369  END IF
370 *
371 * Transform problem to standard eigenvalue problem and solve.
372 *
373  CALL zhpgst( itype, uplo, n, ap, bp, info )
374  CALL zhpevx( jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m,
375  $ w, z, ldz, work, rwork, iwork, ifail, info )
376 *
377  IF( wantz ) THEN
378 *
379 * Backtransform eigenvectors to the original problem.
380 *
381  IF( info.GT.0 )
382  $ m = info - 1
383  IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
384 *
385 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
386 * backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
387 *
388  IF( upper ) THEN
389  trans = 'N'
390  ELSE
391  trans = 'C'
392  END IF
393 *
394  DO 10 j = 1, m
395  CALL ztpsv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
396  $ 1 )
397  10 CONTINUE
398 *
399  ELSE IF( itype.EQ.3 ) THEN
400 *
401 * For B*A*x=(lambda)*x;
402 * backtransform eigenvectors: x = L*y or U**H *y
403 *
404  IF( upper ) THEN
405  trans = 'C'
406  ELSE
407  trans = 'N'
408  END IF
409 *
410  DO 20 j = 1, m
411  CALL ztpmv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
412  $ 1 )
413  20 CONTINUE
414  END IF
415  END IF
416 *
417  RETURN
418 *
419 * End of ZHPGVX
420 *
421  END
zhpgvx
subroutine zhpgvx(ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO)
ZHPGVX
Definition: zhpgvx.f:279
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
ztpsv
subroutine ztpsv(UPLO, TRANS, DIAG, N, AP, X, INCX)
ZTPSV
Definition: ztpsv.f:146
ztpmv
subroutine ztpmv(UPLO, TRANS, DIAG, N, AP, X, INCX)
ZTPMV
Definition: ztpmv.f:144
zhpevx
subroutine zhpevx(JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO)
ZHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matric...
Definition: zhpevx.f:242
zpptrf
subroutine zpptrf(UPLO, N, AP, INFO)
ZPPTRF
Definition: zpptrf.f:121
zhpgst
subroutine zhpgst(ITYPE, UPLO, N, AP, BP, INFO)
ZHPGST
Definition: zhpgst.f:115