LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
chbevd.f
Go to the documentation of this file.
1 *> \brief <b> CHBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHBEVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbevd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbevd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbevd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
22 * LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * REAL RWORK( * ), W( * )
31 * COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CHBEVD computes all the eigenvalues and, optionally, eigenvectors of
41 *> a complex Hermitian band matrix A. If eigenvectors are desired, it
42 *> uses a divide and conquer algorithm.
43 *>
44 *> The divide and conquer algorithm makes very mild assumptions about
45 *> floating point arithmetic. It will work on machines with a guard
46 *> digit in add/subtract, or on those binary machines without guard
47 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
48 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
49 *> without guard digits, but we know of none.
50 *> \endverbatim
51 *
52 * Arguments:
53 * ==========
54 *
55 *> \param[in] JOBZ
56 *> \verbatim
57 *> JOBZ is CHARACTER*1
58 *> = 'N': Compute eigenvalues only;
59 *> = 'V': Compute eigenvalues and eigenvectors.
60 *> \endverbatim
61 *>
62 *> \param[in] UPLO
63 *> \verbatim
64 *> UPLO is CHARACTER*1
65 *> = 'U': Upper triangle of A is stored;
66 *> = 'L': Lower triangle of A is stored.
67 *> \endverbatim
68 *>
69 *> \param[in] N
70 *> \verbatim
71 *> N is INTEGER
72 *> The order of the matrix A. N >= 0.
73 *> \endverbatim
74 *>
75 *> \param[in] KD
76 *> \verbatim
77 *> KD is INTEGER
78 *> The number of superdiagonals of the matrix A if UPLO = 'U',
79 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
80 *> \endverbatim
81 *>
82 *> \param[in,out] AB
83 *> \verbatim
84 *> AB is COMPLEX array, dimension (LDAB, N)
85 *> On entry, the upper or lower triangle of the Hermitian band
86 *> matrix A, stored in the first KD+1 rows of the array. The
87 *> j-th column of A is stored in the j-th column of the array AB
88 *> as follows:
89 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
90 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
91 *>
92 *> On exit, AB is overwritten by values generated during the
93 *> reduction to tridiagonal form. If UPLO = 'U', the first
94 *> superdiagonal and the diagonal of the tridiagonal matrix T
95 *> are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
96 *> the diagonal and first subdiagonal of T are returned in the
97 *> first two rows of AB.
98 *> \endverbatim
99 *>
100 *> \param[in] LDAB
101 *> \verbatim
102 *> LDAB is INTEGER
103 *> The leading dimension of the array AB. LDAB >= KD + 1.
104 *> \endverbatim
105 *>
106 *> \param[out] W
107 *> \verbatim
108 *> W is REAL array, dimension (N)
109 *> If INFO = 0, the eigenvalues in ascending order.
110 *> \endverbatim
111 *>
112 *> \param[out] Z
113 *> \verbatim
114 *> Z is COMPLEX array, dimension (LDZ, N)
115 *> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
116 *> eigenvectors of the matrix A, with the i-th column of Z
117 *> holding the eigenvector associated with W(i).
118 *> If JOBZ = 'N', then Z is not referenced.
119 *> \endverbatim
120 *>
121 *> \param[in] LDZ
122 *> \verbatim
123 *> LDZ is INTEGER
124 *> The leading dimension of the array Z. LDZ >= 1, and if
125 *> JOBZ = 'V', LDZ >= max(1,N).
126 *> \endverbatim
127 *>
128 *> \param[out] WORK
129 *> \verbatim
130 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
131 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
132 *> \endverbatim
133 *>
134 *> \param[in] LWORK
135 *> \verbatim
136 *> LWORK is INTEGER
137 *> The dimension of the array WORK.
138 *> If N <= 1, LWORK must be at least 1.
139 *> If JOBZ = 'N' and N > 1, LWORK must be at least N.
140 *> If JOBZ = 'V' and N > 1, LWORK must be at least 2*N**2.
141 *>
142 *> If LWORK = -1, then a workspace query is assumed; the routine
143 *> only calculates the optimal sizes of the WORK, RWORK and
144 *> IWORK arrays, returns these values as the first entries of
145 *> the WORK, RWORK and IWORK arrays, and no error message
146 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
147 *> \endverbatim
148 *>
149 *> \param[out] RWORK
150 *> \verbatim
151 *> RWORK is REAL array,
152 *> dimension (LRWORK)
153 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
154 *> \endverbatim
155 *>
156 *> \param[in] LRWORK
157 *> \verbatim
158 *> LRWORK is INTEGER
159 *> The dimension of array RWORK.
160 *> If N <= 1, LRWORK must be at least 1.
161 *> If JOBZ = 'N' and N > 1, LRWORK must be at least N.
162 *> If JOBZ = 'V' and N > 1, LRWORK must be at least
163 *> 1 + 5*N + 2*N**2.
164 *>
165 *> If LRWORK = -1, then a workspace query is assumed; the
166 *> routine only calculates the optimal sizes of the WORK, RWORK
167 *> and IWORK arrays, returns these values as the first entries
168 *> of the WORK, RWORK and IWORK arrays, and no error message
169 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
170 *> \endverbatim
171 *>
172 *> \param[out] IWORK
173 *> \verbatim
174 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
175 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
176 *> \endverbatim
177 *>
178 *> \param[in] LIWORK
179 *> \verbatim
180 *> LIWORK is INTEGER
181 *> The dimension of array IWORK.
182 *> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
183 *> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N .
184 *>
185 *> If LIWORK = -1, then a workspace query is assumed; the
186 *> routine only calculates the optimal sizes of the WORK, RWORK
187 *> and IWORK arrays, returns these values as the first entries
188 *> of the WORK, RWORK and IWORK arrays, and no error message
189 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
190 *> \endverbatim
191 *>
192 *> \param[out] INFO
193 *> \verbatim
194 *> INFO is INTEGER
195 *> = 0: successful exit.
196 *> < 0: if INFO = -i, the i-th argument had an illegal value.
197 *> > 0: if INFO = i, the algorithm failed to converge; i
198 *> off-diagonal elements of an intermediate tridiagonal
199 *> form did not converge to zero.
200 *> \endverbatim
201 *
202 * Authors:
203 * ========
204 *
205 *> \author Univ. of Tennessee
206 *> \author Univ. of California Berkeley
207 *> \author Univ. of Colorado Denver
208 *> \author NAG Ltd.
209 *
210 *> \date November 2011
211 *
212 *> \ingroup complexOTHEReigen
213 *
214 * =====================================================================
215  SUBROUTINE chbevd( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
216  $ lwork, rwork, lrwork, iwork, liwork, info )
217 *
218 * -- LAPACK driver routine (version 3.4.0) --
219 * -- LAPACK is a software package provided by Univ. of Tennessee, --
220 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
221 * November 2011
222 *
223 * .. Scalar Arguments ..
224  CHARACTER jobz, uplo
225  INTEGER info, kd, ldab, ldz, liwork, lrwork, lwork, n
226 * ..
227 * .. Array Arguments ..
228  INTEGER iwork( * )
229  REAL rwork( * ), w( * )
230  COMPLEX ab( ldab, * ), work( * ), z( ldz, * )
231 * ..
232 *
233 * =====================================================================
234 *
235 * .. Parameters ..
236  REAL zero, one
237  parameter( zero = 0.0e0, one = 1.0e0 )
238  COMPLEX czero, cone
239  parameter( czero = ( 0.0e0, 0.0e0 ),
240  $ cone = ( 1.0e0, 0.0e0 ) )
241 * ..
242 * .. Local Scalars ..
243  LOGICAL lower, lquery, wantz
244  INTEGER iinfo, imax, inde, indwk2, indwrk, iscale,
245  $ liwmin, llrwk, llwk2, lrwmin, lwmin
246  REAL anrm, bignum, eps, rmax, rmin, safmin, sigma,
247  $ smlnum
248 * ..
249 * .. External Functions ..
250  LOGICAL lsame
251  REAL clanhb, slamch
252  EXTERNAL lsame, clanhb, slamch
253 * ..
254 * .. External Subroutines ..
255  EXTERNAL cgemm, chbtrd, clacpy, clascl, cstedc, sscal,
256  $ ssterf, xerbla
257 * ..
258 * .. Intrinsic Functions ..
259  INTRINSIC sqrt
260 * ..
261 * .. Executable Statements ..
262 *
263 * Test the input parameters.
264 *
265  wantz = lsame( jobz, 'V' )
266  lower = lsame( uplo, 'L' )
267  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 .OR. lrwork.EQ.-1 )
268 *
269  info = 0
270  IF( n.LE.1 ) THEN
271  lwmin = 1
272  lrwmin = 1
273  liwmin = 1
274  ELSE
275  IF( wantz ) THEN
276  lwmin = 2*n**2
277  lrwmin = 1 + 5*n + 2*n**2
278  liwmin = 3 + 5*n
279  ELSE
280  lwmin = n
281  lrwmin = n
282  liwmin = 1
283  END IF
284  END IF
285  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
286  info = -1
287  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
288  info = -2
289  ELSE IF( n.LT.0 ) THEN
290  info = -3
291  ELSE IF( kd.LT.0 ) THEN
292  info = -4
293  ELSE IF( ldab.LT.kd+1 ) THEN
294  info = -6
295  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
296  info = -9
297  END IF
298 *
299  IF( info.EQ.0 ) THEN
300  work( 1 ) = lwmin
301  rwork( 1 ) = lrwmin
302  iwork( 1 ) = liwmin
303 *
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( 'CHBEVD', -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  IF( n.EQ.1 ) THEN
326  w( 1 ) = ab( 1, 1 )
327  IF( wantz )
328  $ z( 1, 1 ) = cone
329  return
330  END IF
331 *
332 * Get machine constants.
333 *
334  safmin = slamch( 'Safe minimum' )
335  eps = slamch( 'Precision' )
336  smlnum = safmin / eps
337  bignum = one / smlnum
338  rmin = sqrt( smlnum )
339  rmax = sqrt( bignum )
340 *
341 * Scale matrix to allowable range, if necessary.
342 *
343  anrm = clanhb( 'M', uplo, n, kd, ab, ldab, rwork )
344  iscale = 0
345  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
346  iscale = 1
347  sigma = rmin / anrm
348  ELSE IF( anrm.GT.rmax ) THEN
349  iscale = 1
350  sigma = rmax / anrm
351  END IF
352  IF( iscale.EQ.1 ) THEN
353  IF( lower ) THEN
354  CALL clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
355  ELSE
356  CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
357  END IF
358  END IF
359 *
360 * Call CHBTRD to reduce Hermitian band matrix to tridiagonal form.
361 *
362  inde = 1
363  indwrk = inde + n
364  indwk2 = 1 + n*n
365  llwk2 = lwork - indwk2 + 1
366  llrwk = lrwork - indwrk + 1
367  CALL chbtrd( jobz, uplo, n, kd, ab, ldab, w, rwork( inde ), z,
368  $ ldz, work, iinfo )
369 *
370 * For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC.
371 *
372  IF( .NOT.wantz ) THEN
373  CALL ssterf( n, w, rwork( inde ), info )
374  ELSE
375  CALL cstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
376  $ llwk2, rwork( indwrk ), llrwk, iwork, liwork,
377  $ info )
378  CALL cgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
379  $ work( indwk2 ), n )
380  CALL clacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
381  END IF
382 *
383 * If matrix was scaled, then rescale eigenvalues appropriately.
384 *
385  IF( iscale.EQ.1 ) THEN
386  IF( info.EQ.0 ) THEN
387  imax = n
388  ELSE
389  imax = info - 1
390  END IF
391  CALL sscal( imax, one / sigma, w, 1 )
392  END IF
393 *
394  work( 1 ) = lwmin
395  rwork( 1 ) = lrwmin
396  iwork( 1 ) = liwmin
397  return
398 *
399 * End of CHBEVD
400 *
401  END