LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
dsbevd.f
Go to the documentation of this file.
1 *> \brief <b> DSBEVD 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 DSBEVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbevd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbevd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbevd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
22 * LWORK, IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DSBEVD computes all the eigenvalues and, optionally, eigenvectors of
40 *> a real symmetric band matrix A. If eigenvectors are desired, it uses
41 *> a divide and conquer algorithm.
42 *>
43 *> The divide and conquer algorithm makes very mild assumptions about
44 *> floating point arithmetic. It will work on machines with a guard
45 *> digit in add/subtract, or on those binary machines without guard
46 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
47 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
48 *> without guard digits, but we know of none.
49 *> \endverbatim
50 *
51 * Arguments:
52 * ==========
53 *
54 *> \param[in] JOBZ
55 *> \verbatim
56 *> JOBZ is CHARACTER*1
57 *> = 'N': Compute eigenvalues only;
58 *> = 'V': Compute eigenvalues and eigenvectors.
59 *> \endverbatim
60 *>
61 *> \param[in] UPLO
62 *> \verbatim
63 *> UPLO is CHARACTER*1
64 *> = 'U': Upper triangle of A is stored;
65 *> = 'L': Lower triangle of A is stored.
66 *> \endverbatim
67 *>
68 *> \param[in] N
69 *> \verbatim
70 *> N is INTEGER
71 *> The order of the matrix A. N >= 0.
72 *> \endverbatim
73 *>
74 *> \param[in] KD
75 *> \verbatim
76 *> KD is INTEGER
77 *> The number of superdiagonals of the matrix A if UPLO = 'U',
78 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
79 *> \endverbatim
80 *>
81 *> \param[in,out] AB
82 *> \verbatim
83 *> AB is DOUBLE PRECISION array, dimension (LDAB, N)
84 *> On entry, the upper or lower triangle of the symmetric band
85 *> matrix A, stored in the first KD+1 rows of the array. The
86 *> j-th column of A is stored in the j-th column of the array AB
87 *> as follows:
88 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
89 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
90 *>
91 *> On exit, AB is overwritten by values generated during the
92 *> reduction to tridiagonal form. If UPLO = 'U', the first
93 *> superdiagonal and the diagonal of the tridiagonal matrix T
94 *> are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
95 *> the diagonal and first subdiagonal of T are returned in the
96 *> first two rows of AB.
97 *> \endverbatim
98 *>
99 *> \param[in] LDAB
100 *> \verbatim
101 *> LDAB is INTEGER
102 *> The leading dimension of the array AB. LDAB >= KD + 1.
103 *> \endverbatim
104 *>
105 *> \param[out] W
106 *> \verbatim
107 *> W is DOUBLE PRECISION array, dimension (N)
108 *> If INFO = 0, the eigenvalues in ascending order.
109 *> \endverbatim
110 *>
111 *> \param[out] Z
112 *> \verbatim
113 *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
114 *> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
115 *> eigenvectors of the matrix A, with the i-th column of Z
116 *> holding the eigenvector associated with W(i).
117 *> If JOBZ = 'N', then Z is not referenced.
118 *> \endverbatim
119 *>
120 *> \param[in] LDZ
121 *> \verbatim
122 *> LDZ is INTEGER
123 *> The leading dimension of the array Z. LDZ >= 1, and if
124 *> JOBZ = 'V', LDZ >= max(1,N).
125 *> \endverbatim
126 *>
127 *> \param[out] WORK
128 *> \verbatim
129 *> WORK is DOUBLE PRECISION array,
130 *> dimension (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 > 2, LWORK must be at least 2*N.
140 *> If JOBZ = 'V' and N > 2, LWORK must be at least
141 *> ( 1 + 5*N + 2*N**2 ).
142 *>
143 *> If LWORK = -1, then a workspace query is assumed; the routine
144 *> only calculates the optimal sizes of the WORK and IWORK
145 *> arrays, returns these values as the first entries of the WORK
146 *> and IWORK arrays, and no error message related to LWORK or
147 *> LIWORK is issued by XERBLA.
148 *> \endverbatim
149 *>
150 *> \param[out] IWORK
151 *> \verbatim
152 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
153 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
154 *> \endverbatim
155 *>
156 *> \param[in] LIWORK
157 *> \verbatim
158 *> LIWORK is INTEGER
159 *> The dimension of the array IWORK.
160 *> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
161 *> If JOBZ = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
162 *>
163 *> If LIWORK = -1, then a workspace query is assumed; the
164 *> routine only calculates the optimal sizes of the WORK and
165 *> IWORK arrays, returns these values as the first entries of
166 *> the WORK and IWORK arrays, and no error message related to
167 *> LWORK or LIWORK is issued by XERBLA.
168 *> \endverbatim
169 *>
170 *> \param[out] INFO
171 *> \verbatim
172 *> INFO is INTEGER
173 *> = 0: successful exit
174 *> < 0: if INFO = -i, the i-th argument had an illegal value
175 *> > 0: if INFO = i, the algorithm failed to converge; i
176 *> off-diagonal elements of an intermediate tridiagonal
177 *> form did not converge to zero.
178 *> \endverbatim
179 *
180 * Authors:
181 * ========
182 *
183 *> \author Univ. of Tennessee
184 *> \author Univ. of California Berkeley
185 *> \author Univ. of Colorado Denver
186 *> \author NAG Ltd.
187 *
188 *> \date November 2011
189 *
190 *> \ingroup doubleOTHEReigen
191 *
192 * =====================================================================
193  SUBROUTINE dsbevd( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
194  $ lwork, iwork, liwork, info )
195 *
196 * -- LAPACK driver routine (version 3.4.0) --
197 * -- LAPACK is a software package provided by Univ. of Tennessee, --
198 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
199 * November 2011
200 *
201 * .. Scalar Arguments ..
202  CHARACTER JOBZ, UPLO
203  INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
204 * ..
205 * .. Array Arguments ..
206  INTEGER IWORK( * )
207  DOUBLE PRECISION AB( ldab, * ), W( * ), WORK( * ), Z( ldz, * )
208 * ..
209 *
210 * =====================================================================
211 *
212 * .. Parameters ..
213  DOUBLE PRECISION ZERO, ONE
214  parameter ( zero = 0.0d+0, one = 1.0d+0 )
215 * ..
216 * .. Local Scalars ..
217  LOGICAL LOWER, LQUERY, WANTZ
218  INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
219  $ llwrk2, lwmin
220  DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
221  $ smlnum
222 * ..
223 * .. External Functions ..
224  LOGICAL LSAME
225  DOUBLE PRECISION DLAMCH, DLANSB
226  EXTERNAL lsame, dlamch, dlansb
227 * ..
228 * .. External Subroutines ..
229  EXTERNAL dgemm, dlacpy, dlascl, dsbtrd, dscal, dstedc,
230  $ dsterf, xerbla
231 * ..
232 * .. Intrinsic Functions ..
233  INTRINSIC sqrt
234 * ..
235 * .. Executable Statements ..
236 *
237 * Test the input parameters.
238 *
239  wantz = lsame( jobz, 'V' )
240  lower = lsame( uplo, 'L' )
241  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
242 *
243  info = 0
244  IF( n.LE.1 ) THEN
245  liwmin = 1
246  lwmin = 1
247  ELSE
248  IF( wantz ) THEN
249  liwmin = 3 + 5*n
250  lwmin = 1 + 5*n + 2*n**2
251  ELSE
252  liwmin = 1
253  lwmin = 2*n
254  END IF
255  END IF
256  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
257  info = -1
258  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
259  info = -2
260  ELSE IF( n.LT.0 ) THEN
261  info = -3
262  ELSE IF( kd.LT.0 ) THEN
263  info = -4
264  ELSE IF( ldab.LT.kd+1 ) THEN
265  info = -6
266  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
267  info = -9
268  END IF
269 *
270  IF( info.EQ.0 ) THEN
271  work( 1 ) = lwmin
272  iwork( 1 ) = liwmin
273 *
274  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
275  info = -11
276  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
277  info = -13
278  END IF
279  END IF
280 *
281  IF( info.NE.0 ) THEN
282  CALL xerbla( 'DSBEVD', -info )
283  RETURN
284  ELSE IF( lquery ) THEN
285  RETURN
286  END IF
287 *
288 * Quick return if possible
289 *
290  IF( n.EQ.0 )
291  $ RETURN
292 *
293  IF( n.EQ.1 ) THEN
294  w( 1 ) = ab( 1, 1 )
295  IF( wantz )
296  $ z( 1, 1 ) = one
297  RETURN
298  END IF
299 *
300 * Get machine constants.
301 *
302  safmin = dlamch( 'Safe minimum' )
303  eps = dlamch( 'Precision' )
304  smlnum = safmin / eps
305  bignum = one / smlnum
306  rmin = sqrt( smlnum )
307  rmax = sqrt( bignum )
308 *
309 * Scale matrix to allowable range, if necessary.
310 *
311  anrm = dlansb( 'M', uplo, n, kd, ab, ldab, work )
312  iscale = 0
313  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
314  iscale = 1
315  sigma = rmin / anrm
316  ELSE IF( anrm.GT.rmax ) THEN
317  iscale = 1
318  sigma = rmax / anrm
319  END IF
320  IF( iscale.EQ.1 ) THEN
321  IF( lower ) THEN
322  CALL dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
323  ELSE
324  CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
325  END IF
326  END IF
327 *
328 * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
329 *
330  inde = 1
331  indwrk = inde + n
332  indwk2 = indwrk + n*n
333  llwrk2 = lwork - indwk2 + 1
334  CALL dsbtrd( jobz, uplo, n, kd, ab, ldab, w, work( inde ), z, ldz,
335  $ work( indwrk ), iinfo )
336 *
337 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
338 *
339  IF( .NOT.wantz ) THEN
340  CALL dsterf( n, w, work( inde ), info )
341  ELSE
342  CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
343  $ work( indwk2 ), llwrk2, iwork, liwork, info )
344  CALL dgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
345  $ zero, work( indwk2 ), n )
346  CALL dlacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
347  END IF
348 *
349 * If matrix was scaled, then rescale eigenvalues appropriately.
350 *
351  IF( iscale.EQ.1 )
352  $ CALL dscal( n, one / sigma, w, 1 )
353 *
354  work( 1 ) = lwmin
355  iwork( 1 ) = liwmin
356  RETURN
357 *
358 * End of DSBEVD
359 *
360  END
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine dsbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
DSBTRD
Definition: dsbtrd.f:165
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dsbevd(JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matric...
Definition: dsbevd.f:195
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:191