LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 *> \ingroup complexOTHEReigen
211 *
212 * =====================================================================
213  SUBROUTINE chbevd( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
214  $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
215 *
216 * -- LAPACK driver routine --
217 * -- LAPACK is a software package provided by Univ. of Tennessee, --
218 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
219 *
220 * .. Scalar Arguments ..
221  CHARACTER JOBZ, UPLO
222  INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
223 * ..
224 * .. Array Arguments ..
225  INTEGER IWORK( * )
226  REAL RWORK( * ), W( * )
227  COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * )
228 * ..
229 *
230 * =====================================================================
231 *
232 * .. Parameters ..
233  REAL ZERO, ONE
234  parameter( zero = 0.0e0, one = 1.0e0 )
235  COMPLEX CZERO, CONE
236  parameter( czero = ( 0.0e0, 0.0e0 ),
237  $ cone = ( 1.0e0, 0.0e0 ) )
238 * ..
239 * .. Local Scalars ..
240  LOGICAL LOWER, LQUERY, WANTZ
241  INTEGER IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE,
242  $ liwmin, llrwk, llwk2, lrwmin, lwmin
243  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
244  $ smlnum
245 * ..
246 * .. External Functions ..
247  LOGICAL LSAME
248  REAL CLANHB, SLAMCH
249  EXTERNAL lsame, clanhb, slamch
250 * ..
251 * .. External Subroutines ..
252  EXTERNAL cgemm, chbtrd, clacpy, clascl, cstedc, sscal,
253  $ ssterf, xerbla
254 * ..
255 * .. Intrinsic Functions ..
256  INTRINSIC sqrt
257 * ..
258 * .. Executable Statements ..
259 *
260 * Test the input parameters.
261 *
262  wantz = lsame( jobz, 'V' )
263  lower = lsame( uplo, 'L' )
264  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 .OR. lrwork.EQ.-1 )
265 *
266  info = 0
267  IF( n.LE.1 ) THEN
268  lwmin = 1
269  lrwmin = 1
270  liwmin = 1
271  ELSE
272  IF( wantz ) THEN
273  lwmin = 2*n**2
274  lrwmin = 1 + 5*n + 2*n**2
275  liwmin = 3 + 5*n
276  ELSE
277  lwmin = n
278  lrwmin = n
279  liwmin = 1
280  END IF
281  END IF
282  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
283  info = -1
284  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
285  info = -2
286  ELSE IF( n.LT.0 ) THEN
287  info = -3
288  ELSE IF( kd.LT.0 ) THEN
289  info = -4
290  ELSE IF( ldab.LT.kd+1 ) THEN
291  info = -6
292  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
293  info = -9
294  END IF
295 *
296  IF( info.EQ.0 ) THEN
297  work( 1 ) = lwmin
298  rwork( 1 ) = lrwmin
299  iwork( 1 ) = liwmin
300 *
301  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
302  info = -11
303  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
304  info = -13
305  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
306  info = -15
307  END IF
308  END IF
309 *
310  IF( info.NE.0 ) THEN
311  CALL xerbla( 'CHBEVD', -info )
312  RETURN
313  ELSE IF( lquery ) THEN
314  RETURN
315  END IF
316 *
317 * Quick return if possible
318 *
319  IF( n.EQ.0 )
320  $ RETURN
321 *
322  IF( n.EQ.1 ) THEN
323  w( 1 ) = real( ab( 1, 1 ) )
324  IF( wantz )
325  $ z( 1, 1 ) = cone
326  RETURN
327  END IF
328 *
329 * Get machine constants.
330 *
331  safmin = slamch( 'Safe minimum' )
332  eps = slamch( 'Precision' )
333  smlnum = safmin / eps
334  bignum = one / smlnum
335  rmin = sqrt( smlnum )
336  rmax = sqrt( bignum )
337 *
338 * Scale matrix to allowable range, if necessary.
339 *
340  anrm = clanhb( 'M', uplo, n, kd, ab, ldab, rwork )
341  iscale = 0
342  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
343  iscale = 1
344  sigma = rmin / anrm
345  ELSE IF( anrm.GT.rmax ) THEN
346  iscale = 1
347  sigma = rmax / anrm
348  END IF
349  IF( iscale.EQ.1 ) THEN
350  IF( lower ) THEN
351  CALL clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
352  ELSE
353  CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
354  END IF
355  END IF
356 *
357 * Call CHBTRD to reduce Hermitian band matrix to tridiagonal form.
358 *
359  inde = 1
360  indwrk = inde + n
361  indwk2 = 1 + n*n
362  llwk2 = lwork - indwk2 + 1
363  llrwk = lrwork - indwrk + 1
364  CALL chbtrd( jobz, uplo, n, kd, ab, ldab, w, rwork( inde ), z,
365  $ ldz, work, iinfo )
366 *
367 * For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC.
368 *
369  IF( .NOT.wantz ) THEN
370  CALL ssterf( n, w, rwork( inde ), info )
371  ELSE
372  CALL cstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
373  $ llwk2, rwork( indwrk ), llrwk, iwork, liwork,
374  $ info )
375  CALL cgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
376  $ work( indwk2 ), n )
377  CALL clacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
378  END IF
379 *
380 * If matrix was scaled, then rescale eigenvalues appropriately.
381 *
382  IF( iscale.EQ.1 ) THEN
383  IF( info.EQ.0 ) THEN
384  imax = n
385  ELSE
386  imax = info - 1
387  END IF
388  CALL sscal( imax, one / sigma, w, 1 )
389  END IF
390 *
391  work( 1 ) = lwmin
392  rwork( 1 ) = lrwmin
393  iwork( 1 ) = liwmin
394  RETURN
395 *
396 * End of CHBEVD
397 *
398  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 clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:143
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 chbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
CHBTRD
Definition: chbtrd.f:163
subroutine cstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CSTEDC
Definition: cstedc.f:212
subroutine chbevd(JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CHBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition: chbevd.f:215
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79