LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
chbgvd.f
Go to the documentation of this file.
1 *> \brief \b CHBGVD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHBGVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbgvd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbgvd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbgvd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
22 * Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK,
23 * LIWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, UPLO
27 * INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK,
28 * $ LWORK, N
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IWORK( * )
32 * REAL RWORK( * ), W( * )
33 * COMPLEX AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
34 * $ Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> CHBGVD computes all the eigenvalues, and optionally, the eigenvectors
44 *> of a complex generalized Hermitian-definite banded eigenproblem, of
45 *> the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian
46 *> and banded, and B is also positive definite. If eigenvectors are
47 *> desired, it uses a divide and conquer algorithm.
48 *>
49 *> The divide and conquer algorithm makes very mild assumptions about
50 *> floating point arithmetic. It will work on machines with a guard
51 *> digit in add/subtract, or on those binary machines without guard
52 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
53 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
54 *> without guard digits, but we know of none.
55 *> \endverbatim
56 *
57 * Arguments:
58 * ==========
59 *
60 *> \param[in] JOBZ
61 *> \verbatim
62 *> JOBZ is CHARACTER*1
63 *> = 'N': Compute eigenvalues only;
64 *> = 'V': Compute eigenvalues and eigenvectors.
65 *> \endverbatim
66 *>
67 *> \param[in] UPLO
68 *> \verbatim
69 *> UPLO is CHARACTER*1
70 *> = 'U': Upper triangles of A and B are stored;
71 *> = 'L': Lower triangles of A and B are stored.
72 *> \endverbatim
73 *>
74 *> \param[in] N
75 *> \verbatim
76 *> N is INTEGER
77 *> The order of the matrices A and B. N >= 0.
78 *> \endverbatim
79 *>
80 *> \param[in] KA
81 *> \verbatim
82 *> KA is INTEGER
83 *> The number of superdiagonals of the matrix A if UPLO = 'U',
84 *> or the number of subdiagonals if UPLO = 'L'. KA >= 0.
85 *> \endverbatim
86 *>
87 *> \param[in] KB
88 *> \verbatim
89 *> KB is INTEGER
90 *> The number of superdiagonals of the matrix B if UPLO = 'U',
91 *> or the number of subdiagonals if UPLO = 'L'. KB >= 0.
92 *> \endverbatim
93 *>
94 *> \param[in,out] AB
95 *> \verbatim
96 *> AB is COMPLEX array, dimension (LDAB, N)
97 *> On entry, the upper or lower triangle of the Hermitian band
98 *> matrix A, stored in the first ka+1 rows of the array. The
99 *> j-th column of A is stored in the j-th column of the array AB
100 *> as follows:
101 *> if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
102 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
103 *>
104 *> On exit, the contents of AB are destroyed.
105 *> \endverbatim
106 *>
107 *> \param[in] LDAB
108 *> \verbatim
109 *> LDAB is INTEGER
110 *> The leading dimension of the array AB. LDAB >= KA+1.
111 *> \endverbatim
112 *>
113 *> \param[in,out] BB
114 *> \verbatim
115 *> BB is COMPLEX array, dimension (LDBB, N)
116 *> On entry, the upper or lower triangle of the Hermitian band
117 *> matrix B, stored in the first kb+1 rows of the array. The
118 *> j-th column of B is stored in the j-th column of the array BB
119 *> as follows:
120 *> if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
121 *> if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb).
122 *>
123 *> On exit, the factor S from the split Cholesky factorization
124 *> B = S**H*S, as returned by CPBSTF.
125 *> \endverbatim
126 *>
127 *> \param[in] LDBB
128 *> \verbatim
129 *> LDBB is INTEGER
130 *> The leading dimension of the array BB. LDBB >= KB+1.
131 *> \endverbatim
132 *>
133 *> \param[out] W
134 *> \verbatim
135 *> W is REAL array, dimension (N)
136 *> If INFO = 0, the eigenvalues in ascending order.
137 *> \endverbatim
138 *>
139 *> \param[out] Z
140 *> \verbatim
141 *> Z is COMPLEX array, dimension (LDZ, N)
142 *> If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
143 *> eigenvectors, with the i-th column of Z holding the
144 *> eigenvector associated with W(i). The eigenvectors are
145 *> normalized so that Z**H*B*Z = I.
146 *> If JOBZ = 'N', then Z is not referenced.
147 *> \endverbatim
148 *>
149 *> \param[in] LDZ
150 *> \verbatim
151 *> LDZ is INTEGER
152 *> The leading dimension of the array Z. LDZ >= 1, and if
153 *> JOBZ = 'V', LDZ >= N.
154 *> \endverbatim
155 *>
156 *> \param[out] WORK
157 *> \verbatim
158 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
159 *> On exit, if INFO=0, WORK(1) returns the optimal LWORK.
160 *> \endverbatim
161 *>
162 *> \param[in] LWORK
163 *> \verbatim
164 *> LWORK is INTEGER
165 *> The dimension of the array WORK.
166 *> If N <= 1, LWORK >= 1.
167 *> If JOBZ = 'N' and N > 1, LWORK >= N.
168 *> If JOBZ = 'V' and N > 1, LWORK >= 2*N**2.
169 *>
170 *> If LWORK = -1, then a workspace query is assumed; the routine
171 *> only calculates the optimal sizes of the WORK, RWORK and
172 *> IWORK arrays, returns these values as the first entries of
173 *> the WORK, RWORK and IWORK arrays, and no error message
174 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
175 *> \endverbatim
176 *>
177 *> \param[out] RWORK
178 *> \verbatim
179 *> RWORK is REAL array, dimension (MAX(1,LRWORK))
180 *> On exit, if INFO=0, RWORK(1) returns the optimal LRWORK.
181 *> \endverbatim
182 *>
183 *> \param[in] LRWORK
184 *> \verbatim
185 *> LRWORK is INTEGER
186 *> The dimension of array RWORK.
187 *> If N <= 1, LRWORK >= 1.
188 *> If JOBZ = 'N' and N > 1, LRWORK >= N.
189 *> If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
190 *>
191 *> If LRWORK = -1, then a workspace query is assumed; the
192 *> routine only calculates the optimal sizes of the WORK, RWORK
193 *> and IWORK arrays, returns these values as the first entries
194 *> of the WORK, RWORK and IWORK arrays, and no error message
195 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
196 *> \endverbatim
197 *>
198 *> \param[out] IWORK
199 *> \verbatim
200 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
201 *> On exit, if INFO=0, IWORK(1) returns the optimal LIWORK.
202 *> \endverbatim
203 *>
204 *> \param[in] LIWORK
205 *> \verbatim
206 *> LIWORK is INTEGER
207 *> The dimension of array IWORK.
208 *> If JOBZ = 'N' or N <= 1, LIWORK >= 1.
209 *> If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
210 *>
211 *> If LIWORK = -1, then a workspace query is assumed; the
212 *> routine only calculates the optimal sizes of the WORK, RWORK
213 *> and IWORK arrays, returns these values as the first entries
214 *> of the WORK, RWORK and IWORK arrays, and no error message
215 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
216 *> \endverbatim
217 *>
218 *> \param[out] INFO
219 *> \verbatim
220 *> INFO is INTEGER
221 *> = 0: successful exit
222 *> < 0: if INFO = -i, the i-th argument had an illegal value
223 *> > 0: if INFO = i, and i is:
224 *> <= N: the algorithm failed to converge:
225 *> i off-diagonal elements of an intermediate
226 *> tridiagonal form did not converge to zero;
227 *> > N: if INFO = N + i, for 1 <= i <= N, then CPBSTF
228 *> returned INFO = i: B is not positive definite.
229 *> The factorization of B could not be completed and
230 *> no eigenvalues or eigenvectors were computed.
231 *> \endverbatim
232 *
233 * Authors:
234 * ========
235 *
236 *> \author Univ. of Tennessee
237 *> \author Univ. of California Berkeley
238 *> \author Univ. of Colorado Denver
239 *> \author NAG Ltd.
240 *
241 *> \ingroup complexOTHEReigen
242 *
243 *> \par Contributors:
244 * ==================
245 *>
246 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
247 *
248 * =====================================================================
249  SUBROUTINE chbgvd( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
250  $ Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK,
251  $ LIWORK, INFO )
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 *
404  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
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
subroutine cpbstf(UPLO, N, KD, AB, LDAB, INFO)
CPBSTF
Definition: cpbstf.f:153
subroutine chbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
CHBTRD
Definition: chbtrd.f:163
subroutine chbgst(VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X, LDX, WORK, RWORK, INFO)
CHBGST
Definition: chbgst.f:165
subroutine cstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CSTEDC
Definition: cstedc.f:212
subroutine chbgvd(JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CHBGVD
Definition: chbgvd.f:252