LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
ssterf
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
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
cpbstf
subroutine cpbstf(UPLO, N, KD, AB, LDAB, INFO)
CPBSTF
Definition: cpbstf.f:153
chbtrd
subroutine chbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
CHBTRD
Definition: chbtrd.f:163
chbgst
subroutine chbgst(VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X, LDX, WORK, RWORK, INFO)
CHBGST
Definition: chbgst.f:165
cstedc
subroutine cstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CSTEDC
Definition: cstedc.f:212
chbgvd
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