LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dspevd.f
Go to the documentation of this file.
1 *> \brief <b> DSPEVD 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 DSPEVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dspevd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dspevd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dspevd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
22 * IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, LDZ, LIWORK, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DSPEVD computes all the eigenvalues and, optionally, eigenvectors
40 *> of a real symmetric matrix A in packed storage. If eigenvectors are
41 *> desired, it uses 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,out] AP
75 *> \verbatim
76 *> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
77 *> On entry, the upper or lower triangle of the symmetric matrix
78 *> A, packed columnwise in a linear array. The j-th column of A
79 *> is stored in the array AP as follows:
80 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
81 *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
82 *>
83 *> On exit, AP is overwritten by values generated during the
84 *> reduction to tridiagonal form. If UPLO = 'U', the diagonal
85 *> and first superdiagonal of the tridiagonal matrix T overwrite
86 *> the corresponding elements of A, and if UPLO = 'L', the
87 *> diagonal and first subdiagonal of T overwrite the
88 *> corresponding elements of A.
89 *> \endverbatim
90 *>
91 *> \param[out] W
92 *> \verbatim
93 *> W is DOUBLE PRECISION array, dimension (N)
94 *> If INFO = 0, the eigenvalues in ascending order.
95 *> \endverbatim
96 *>
97 *> \param[out] Z
98 *> \verbatim
99 *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
100 *> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
101 *> eigenvectors of the matrix A, with the i-th column of Z
102 *> holding the eigenvector associated with W(i).
103 *> If JOBZ = 'N', then Z is not referenced.
104 *> \endverbatim
105 *>
106 *> \param[in] LDZ
107 *> \verbatim
108 *> LDZ is INTEGER
109 *> The leading dimension of the array Z. LDZ >= 1, and if
110 *> JOBZ = 'V', LDZ >= max(1,N).
111 *> \endverbatim
112 *>
113 *> \param[out] WORK
114 *> \verbatim
115 *> WORK is DOUBLE PRECISION array,
116 *> dimension (LWORK)
117 *> On exit, if INFO = 0, WORK(1) returns the required LWORK.
118 *> \endverbatim
119 *>
120 *> \param[in] LWORK
121 *> \verbatim
122 *> LWORK is INTEGER
123 *> The dimension of the array WORK.
124 *> If N <= 1, LWORK must be at least 1.
125 *> If JOBZ = 'N' and N > 1, LWORK must be at least 2*N.
126 *> If JOBZ = 'V' and N > 1, LWORK must be at least
127 *> 1 + 6*N + N**2.
128 *>
129 *> If LWORK = -1, then a workspace query is assumed; the routine
130 *> only calculates the required sizes of the WORK and IWORK
131 *> arrays, returns these values as the first entries of the WORK
132 *> and IWORK arrays, and no error message related to LWORK or
133 *> LIWORK is issued by XERBLA.
134 *> \endverbatim
135 *>
136 *> \param[out] IWORK
137 *> \verbatim
138 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
139 *> On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
140 *> \endverbatim
141 *>
142 *> \param[in] LIWORK
143 *> \verbatim
144 *> LIWORK is INTEGER
145 *> The dimension of the array IWORK.
146 *> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
147 *> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
148 *>
149 *> If LIWORK = -1, then a workspace query is assumed; the
150 *> routine only calculates the required sizes of the WORK and
151 *> IWORK arrays, returns these values as the first entries of
152 *> the WORK and IWORK arrays, and no error message related to
153 *> LWORK or LIWORK is issued by XERBLA.
154 *> \endverbatim
155 *>
156 *> \param[out] INFO
157 *> \verbatim
158 *> INFO is INTEGER
159 *> = 0: successful exit
160 *> < 0: if INFO = -i, the i-th argument had an illegal value.
161 *> > 0: if INFO = i, the algorithm failed to converge; i
162 *> off-diagonal elements of an intermediate tridiagonal
163 *> form did not converge to zero.
164 *> \endverbatim
165 *
166 * Authors:
167 * ========
168 *
169 *> \author Univ. of Tennessee
170 *> \author Univ. of California Berkeley
171 *> \author Univ. of Colorado Denver
172 *> \author NAG Ltd.
173 *
174 *> \date November 2011
175 *
176 *> \ingroup doubleOTHEReigen
177 *
178 * =====================================================================
179  SUBROUTINE dspevd( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
180  $ iwork, liwork, info )
181 *
182 * -- LAPACK driver routine (version 3.4.0) --
183 * -- LAPACK is a software package provided by Univ. of Tennessee, --
184 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185 * November 2011
186 *
187 * .. Scalar Arguments ..
188  CHARACTER jobz, uplo
189  INTEGER info, ldz, liwork, lwork, n
190 * ..
191 * .. Array Arguments ..
192  INTEGER iwork( * )
193  DOUBLE PRECISION ap( * ), w( * ), work( * ), z( ldz, * )
194 * ..
195 *
196 * =====================================================================
197 *
198 * .. Parameters ..
199  DOUBLE PRECISION zero, one
200  parameter( zero = 0.0d+0, one = 1.0d+0 )
201 * ..
202 * .. Local Scalars ..
203  LOGICAL lquery, wantz
204  INTEGER iinfo, inde, indtau, indwrk, iscale, liwmin,
205  $ llwork, lwmin
206  DOUBLE PRECISION anrm, bignum, eps, rmax, rmin, safmin, sigma,
207  $ smlnum
208 * ..
209 * .. External Functions ..
210  LOGICAL lsame
211  DOUBLE PRECISION dlamch, dlansp
212  EXTERNAL lsame, dlamch, dlansp
213 * ..
214 * .. External Subroutines ..
215  EXTERNAL dopmtr, dscal, dsptrd, dstedc, dsterf, xerbla
216 * ..
217 * .. Intrinsic Functions ..
218  INTRINSIC sqrt
219 * ..
220 * .. Executable Statements ..
221 *
222 * Test the input parameters.
223 *
224  wantz = lsame( jobz, 'V' )
225  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
226 *
227  info = 0
228  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
229  info = -1
230  ELSE IF( .NOT.( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) )
231  $ THEN
232  info = -2
233  ELSE IF( n.LT.0 ) THEN
234  info = -3
235  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
236  info = -7
237  END IF
238 *
239  IF( info.EQ.0 ) THEN
240  IF( n.LE.1 ) THEN
241  liwmin = 1
242  lwmin = 1
243  ELSE
244  IF( wantz ) THEN
245  liwmin = 3 + 5*n
246  lwmin = 1 + 6*n + n**2
247  ELSE
248  liwmin = 1
249  lwmin = 2*n
250  END IF
251  END IF
252  iwork( 1 ) = liwmin
253  work( 1 ) = lwmin
254 *
255  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
256  info = -9
257  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
258  info = -11
259  END IF
260  END IF
261 *
262  IF( info.NE.0 ) THEN
263  CALL xerbla( 'DSPEVD', -info )
264  return
265  ELSE IF( lquery ) THEN
266  return
267  END IF
268 *
269 * Quick return if possible
270 *
271  IF( n.EQ.0 )
272  $ return
273 *
274  IF( n.EQ.1 ) THEN
275  w( 1 ) = ap( 1 )
276  IF( wantz )
277  $ z( 1, 1 ) = one
278  return
279  END IF
280 *
281 * Get machine constants.
282 *
283  safmin = dlamch( 'Safe minimum' )
284  eps = dlamch( 'Precision' )
285  smlnum = safmin / eps
286  bignum = one / smlnum
287  rmin = sqrt( smlnum )
288  rmax = sqrt( bignum )
289 *
290 * Scale matrix to allowable range, if necessary.
291 *
292  anrm = dlansp( 'M', uplo, n, ap, work )
293  iscale = 0
294  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
295  iscale = 1
296  sigma = rmin / anrm
297  ELSE IF( anrm.GT.rmax ) THEN
298  iscale = 1
299  sigma = rmax / anrm
300  END IF
301  IF( iscale.EQ.1 ) THEN
302  CALL dscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
303  END IF
304 *
305 * Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
306 *
307  inde = 1
308  indtau = inde + n
309  CALL dsptrd( uplo, n, ap, w, work( inde ), work( indtau ), iinfo )
310 *
311 * For eigenvalues only, call DSTERF. For eigenvectors, first call
312 * DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
313 * tridiagonal matrix, then call DOPMTR to multiply it by the
314 * Householder transformations represented in AP.
315 *
316  IF( .NOT.wantz ) THEN
317  CALL dsterf( n, w, work( inde ), info )
318  ELSE
319  indwrk = indtau + n
320  llwork = lwork - indwrk + 1
321  CALL dstedc( 'I', n, w, work( inde ), z, ldz, work( indwrk ),
322  $ llwork, iwork, liwork, info )
323  CALL dopmtr( 'L', uplo, 'N', n, n, ap, work( indtau ), z, ldz,
324  $ work( indwrk ), iinfo )
325  END IF
326 *
327 * If matrix was scaled, then rescale eigenvalues appropriately.
328 *
329  IF( iscale.EQ.1 )
330  $ CALL dscal( n, one / sigma, w, 1 )
331 *
332  work( 1 ) = lwmin
333  iwork( 1 ) = liwmin
334  return
335 *
336 * End of DSPEVD
337 *
338  END