LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
chbgvd.f
Go to the documentation of this file.
1 *> \brief \b CHBGST
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 *> \date November 2011
242 *
243 *> \ingroup complexOTHEReigen
244 *
245 *> \par Contributors:
246 * ==================
247 *>
248 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
249 *
250 * =====================================================================
251  SUBROUTINE chbgvd( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
252  $ z, ldz, work, lwork, rwork, lrwork, iwork,
253  $ liwork, info )
254 *
255 * -- LAPACK driver routine (version 3.4.0) --
256 * -- LAPACK is a software package provided by Univ. of Tennessee, --
257 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
258 * November 2011
259 *
260 * .. Scalar Arguments ..
261  CHARACTER jobz, uplo
262  INTEGER info, ka, kb, ldab, ldbb, ldz, liwork, lrwork,
263  $ lwork, n
264 * ..
265 * .. Array Arguments ..
266  INTEGER iwork( * )
267  REAL rwork( * ), w( * )
268  COMPLEX ab( ldab, * ), bb( ldbb, * ), work( * ),
269  $ z( ldz, * )
270 * ..
271 *
272 * =====================================================================
273 *
274 * .. Parameters ..
275  COMPLEX cone, czero
276  parameter( cone = ( 1.0e+0, 0.0e+0 ),
277  $ czero = ( 0.0e+0, 0.0e+0 ) )
278 * ..
279 * .. Local Scalars ..
280  LOGICAL lquery, upper, wantz
281  CHARACTER vect
282  INTEGER iinfo, inde, indwk2, indwrk, liwmin, llrwk,
283  $ llwk2, lrwmin, lwmin
284 * ..
285 * .. External Functions ..
286  LOGICAL lsame
287  EXTERNAL lsame
288 * ..
289 * .. External Subroutines ..
290  EXTERNAL ssterf, xerbla, cgemm, chbgst, chbtrd, clacpy,
291  $ cpbstf, cstedc
292 * ..
293 * .. Executable Statements ..
294 *
295 * Test the input parameters.
296 *
297  wantz = lsame( jobz, 'V' )
298  upper = lsame( uplo, 'U' )
299  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
300 *
301  info = 0
302  IF( n.LE.1 ) THEN
303  lwmin = 1+n
304  lrwmin = 1+n
305  liwmin = 1
306  ELSE IF( wantz ) THEN
307  lwmin = 2*n**2
308  lrwmin = 1 + 5*n + 2*n**2
309  liwmin = 3 + 5*n
310  ELSE
311  lwmin = n
312  lrwmin = n
313  liwmin = 1
314  END IF
315  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
316  info = -1
317  ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
318  info = -2
319  ELSE IF( n.LT.0 ) THEN
320  info = -3
321  ELSE IF( ka.LT.0 ) THEN
322  info = -4
323  ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
324  info = -5
325  ELSE IF( ldab.LT.ka+1 ) THEN
326  info = -7
327  ELSE IF( ldbb.LT.kb+1 ) THEN
328  info = -9
329  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
330  info = -12
331  END IF
332 *
333  IF( info.EQ.0 ) THEN
334  work( 1 ) = lwmin
335  rwork( 1 ) = lrwmin
336  iwork( 1 ) = liwmin
337 *
338  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
339  info = -14
340  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
341  info = -16
342  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
343  info = -18
344  END IF
345  END IF
346 *
347  IF( info.NE.0 ) THEN
348  CALL xerbla( 'CHBGVD', -info )
349  return
350  ELSE IF( lquery ) THEN
351  return
352  END IF
353 *
354 * Quick return if possible
355 *
356  IF( n.EQ.0 )
357  $ return
358 *
359 * Form a split Cholesky factorization of B.
360 *
361  CALL cpbstf( uplo, n, kb, bb, ldbb, info )
362  IF( info.NE.0 ) THEN
363  info = n + info
364  return
365  END IF
366 *
367 * Transform problem to standard eigenvalue problem.
368 *
369  inde = 1
370  indwrk = inde + n
371  indwk2 = 1 + n*n
372  llwk2 = lwork - indwk2 + 2
373  llrwk = lrwork - indwrk + 2
374  CALL chbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,
375  $ work, rwork( indwrk ), iinfo )
376 *
377 * Reduce Hermitian band matrix to tridiagonal form.
378 *
379  IF( wantz ) THEN
380  vect = 'U'
381  ELSE
382  vect = 'N'
383  END IF
384  CALL chbtrd( vect, uplo, n, ka, ab, ldab, w, rwork( inde ), z,
385  $ ldz, work, iinfo )
386 *
387 * For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC.
388 *
389  IF( .NOT.wantz ) THEN
390  CALL ssterf( n, w, rwork( inde ), info )
391  ELSE
392  CALL cstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
393  $ llwk2, rwork( indwrk ), llrwk, iwork, liwork,
394  $ info )
395  CALL cgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
396  $ work( indwk2 ), n )
397  CALL clacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
398  END IF
399 *
400  work( 1 ) = lwmin
401  rwork( 1 ) = lrwmin
402  iwork( 1 ) = liwmin
403  return
404 *
405 * End of CHBGVD
406 *
407  END