LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dspevx.f
Go to the documentation of this file.
1 *> \brief <b> DSPEVX 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 DSPEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dspevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dspevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dspevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
22 * ABSTOL, M, W, Z, LDZ, WORK, IWORK, IFAIL,
23 * INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, IU, LDZ, M, N
28 * DOUBLE PRECISION ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DSPEVX computes selected eigenvalues and, optionally, eigenvectors
42 *> of a real symmetric matrix A in packed storage. Eigenvalues/vectors
43 *> can be selected by specifying either a range of values or a range of
44 *> indices for the desired eigenvalues.
45 *> \endverbatim
46 *
47 * Arguments:
48 * ==========
49 *
50 *> \param[in] JOBZ
51 *> \verbatim
52 *> JOBZ is CHARACTER*1
53 *> = 'N': Compute eigenvalues only;
54 *> = 'V': Compute eigenvalues and eigenvectors.
55 *> \endverbatim
56 *>
57 *> \param[in] RANGE
58 *> \verbatim
59 *> RANGE is CHARACTER*1
60 *> = 'A': all eigenvalues will be found;
61 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
62 *> will be found;
63 *> = 'I': the IL-th through IU-th eigenvalues will be found.
64 *> \endverbatim
65 *>
66 *> \param[in] UPLO
67 *> \verbatim
68 *> UPLO is CHARACTER*1
69 *> = 'U': Upper triangle of A is stored;
70 *> = 'L': Lower triangle of A is stored.
71 *> \endverbatim
72 *>
73 *> \param[in] N
74 *> \verbatim
75 *> N is INTEGER
76 *> The order of the matrix A. N >= 0.
77 *> \endverbatim
78 *>
79 *> \param[in,out] AP
80 *> \verbatim
81 *> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
82 *> On entry, the upper or lower triangle of the symmetric matrix
83 *> A, packed columnwise in a linear array. The j-th column of A
84 *> is stored in the array AP as follows:
85 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
86 *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
87 *>
88 *> On exit, AP is overwritten by values generated during the
89 *> reduction to tridiagonal form. If UPLO = 'U', the diagonal
90 *> and first superdiagonal of the tridiagonal matrix T overwrite
91 *> the corresponding elements of A, and if UPLO = 'L', the
92 *> diagonal and first subdiagonal of T overwrite the
93 *> corresponding elements of A.
94 *> \endverbatim
95 *>
96 *> \param[in] VL
97 *> \verbatim
98 *> VL is DOUBLE PRECISION
99 *> \endverbatim
100 *>
101 *> \param[in] VU
102 *> \verbatim
103 *> VU is DOUBLE PRECISION
104 *> If RANGE='V', the lower and upper bounds of the interval to
105 *> be searched for eigenvalues. VL < VU.
106 *> Not referenced if RANGE = 'A' or 'I'.
107 *> \endverbatim
108 *>
109 *> \param[in] IL
110 *> \verbatim
111 *> IL is INTEGER
112 *> \endverbatim
113 *>
114 *> \param[in] IU
115 *> \verbatim
116 *> IU is INTEGER
117 *> If RANGE='I', the indices (in ascending order) of the
118 *> smallest and largest eigenvalues to be returned.
119 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
120 *> Not referenced if RANGE = 'A' or 'V'.
121 *> \endverbatim
122 *>
123 *> \param[in] ABSTOL
124 *> \verbatim
125 *> ABSTOL is DOUBLE PRECISION
126 *> The absolute error tolerance for the eigenvalues.
127 *> An approximate eigenvalue is accepted as converged
128 *> when it is determined to lie in an interval [a,b]
129 *> of width less than or equal to
130 *>
131 *> ABSTOL + EPS * max( |a|,|b| ) ,
132 *>
133 *> where EPS is the machine precision. If ABSTOL is less than
134 *> or equal to zero, then EPS*|T| will be used in its place,
135 *> where |T| is the 1-norm of the tridiagonal matrix obtained
136 *> by reducing AP to tridiagonal form.
137 *>
138 *> Eigenvalues will be computed most accurately when ABSTOL is
139 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
140 *> If this routine returns with INFO>0, indicating that some
141 *> eigenvectors did not converge, try setting ABSTOL to
142 *> 2*DLAMCH('S').
143 *>
144 *> See "Computing Small Singular Values of Bidiagonal Matrices
145 *> with Guaranteed High Relative Accuracy," by Demmel and
146 *> Kahan, LAPACK Working Note #3.
147 *> \endverbatim
148 *>
149 *> \param[out] M
150 *> \verbatim
151 *> M is INTEGER
152 *> The total number of eigenvalues found. 0 <= M <= N.
153 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
154 *> \endverbatim
155 *>
156 *> \param[out] W
157 *> \verbatim
158 *> W is DOUBLE PRECISION array, dimension (N)
159 *> If INFO = 0, the selected eigenvalues in ascending order.
160 *> \endverbatim
161 *>
162 *> \param[out] Z
163 *> \verbatim
164 *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
165 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
166 *> contain the orthonormal eigenvectors of the matrix A
167 *> corresponding to the selected eigenvalues, with the i-th
168 *> column of Z holding the eigenvector associated with W(i).
169 *> If an eigenvector fails to converge, then that column of Z
170 *> contains the latest approximation to the eigenvector, and the
171 *> index of the eigenvector is returned in IFAIL.
172 *> If JOBZ = 'N', then Z is not referenced.
173 *> Note: the user must ensure that at least max(1,M) columns are
174 *> supplied in the array Z; if RANGE = 'V', the exact value of M
175 *> is not known in advance and an upper bound must be used.
176 *> \endverbatim
177 *>
178 *> \param[in] LDZ
179 *> \verbatim
180 *> LDZ is INTEGER
181 *> The leading dimension of the array Z. LDZ >= 1, and if
182 *> JOBZ = 'V', LDZ >= max(1,N).
183 *> \endverbatim
184 *>
185 *> \param[out] WORK
186 *> \verbatim
187 *> WORK is DOUBLE PRECISION array, dimension (8*N)
188 *> \endverbatim
189 *>
190 *> \param[out] IWORK
191 *> \verbatim
192 *> IWORK is INTEGER array, dimension (5*N)
193 *> \endverbatim
194 *>
195 *> \param[out] IFAIL
196 *> \verbatim
197 *> IFAIL is INTEGER array, dimension (N)
198 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
199 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
200 *> indices of the eigenvectors that failed to converge.
201 *> If JOBZ = 'N', then IFAIL is not referenced.
202 *> \endverbatim
203 *>
204 *> \param[out] INFO
205 *> \verbatim
206 *> INFO is INTEGER
207 *> = 0: successful exit
208 *> < 0: if INFO = -i, the i-th argument had an illegal value
209 *> > 0: if INFO = i, then i eigenvectors failed to converge.
210 *> Their indices are stored in array IFAIL.
211 *> \endverbatim
212 *
213 * Authors:
214 * ========
215 *
216 *> \author Univ. of Tennessee
217 *> \author Univ. of California Berkeley
218 *> \author Univ. of Colorado Denver
219 *> \author NAG Ltd.
220 *
221 *> \date November 2011
222 *
223 *> \ingroup doubleOTHEReigen
224 *
225 * =====================================================================
226  SUBROUTINE dspevx( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
227  $ abstol, m, w, z, ldz, work, iwork, ifail,
228  $ info )
229 *
230 * -- LAPACK driver routine (version 3.4.0) --
231 * -- LAPACK is a software package provided by Univ. of Tennessee, --
232 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233 * November 2011
234 *
235 * .. Scalar Arguments ..
236  CHARACTER jobz, range, uplo
237  INTEGER il, info, iu, ldz, m, n
238  DOUBLE PRECISION abstol, vl, vu
239 * ..
240 * .. Array Arguments ..
241  INTEGER ifail( * ), iwork( * )
242  DOUBLE PRECISION ap( * ), w( * ), work( * ), z( ldz, * )
243 * ..
244 *
245 * =====================================================================
246 *
247 * .. Parameters ..
248  DOUBLE PRECISION zero, one
249  parameter( zero = 0.0d0, one = 1.0d0 )
250 * ..
251 * .. Local Scalars ..
252  LOGICAL alleig, indeig, test, valeig, wantz
253  CHARACTER order
254  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
255  $ indisp, indiwo, indtau, indwrk, iscale, itmp1,
256  $ j, jj, nsplit
257  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
258  $ sigma, smlnum, tmp1, vll, vuu
259 * ..
260 * .. External Functions ..
261  LOGICAL lsame
262  DOUBLE PRECISION dlamch, dlansp
263  EXTERNAL lsame, dlamch, dlansp
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL dcopy, dopgtr, dopmtr, dscal, dsptrd, dstebz,
267  $ dstein, dsteqr, dsterf, dswap, xerbla
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC max, min, sqrt
271 * ..
272 * .. Executable Statements ..
273 *
274 * Test the input parameters.
275 *
276  wantz = lsame( jobz, 'V' )
277  alleig = lsame( range, 'A' )
278  valeig = lsame( range, 'V' )
279  indeig = lsame( range, 'I' )
280 *
281  info = 0
282  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
283  info = -1
284  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
285  info = -2
286  ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
287  $ THEN
288  info = -3
289  ELSE IF( n.LT.0 ) THEN
290  info = -4
291  ELSE
292  IF( valeig ) THEN
293  IF( n.GT.0 .AND. vu.LE.vl )
294  $ info = -7
295  ELSE IF( indeig ) THEN
296  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
297  info = -8
298  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
299  info = -9
300  END IF
301  END IF
302  END IF
303  IF( info.EQ.0 ) THEN
304  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
305  $ info = -14
306  END IF
307 *
308  IF( info.NE.0 ) THEN
309  CALL xerbla( 'DSPEVX', -info )
310  return
311  END IF
312 *
313 * Quick return if possible
314 *
315  m = 0
316  IF( n.EQ.0 )
317  $ return
318 *
319  IF( n.EQ.1 ) THEN
320  IF( alleig .OR. indeig ) THEN
321  m = 1
322  w( 1 ) = ap( 1 )
323  ELSE
324  IF( vl.LT.ap( 1 ) .AND. vu.GE.ap( 1 ) ) THEN
325  m = 1
326  w( 1 ) = ap( 1 )
327  END IF
328  END IF
329  IF( wantz )
330  $ z( 1, 1 ) = one
331  return
332  END IF
333 *
334 * Get machine constants.
335 *
336  safmin = dlamch( 'Safe minimum' )
337  eps = dlamch( 'Precision' )
338  smlnum = safmin / eps
339  bignum = one / smlnum
340  rmin = sqrt( smlnum )
341  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
342 *
343 * Scale matrix to allowable range, if necessary.
344 *
345  iscale = 0
346  abstll = abstol
347  IF( valeig ) THEN
348  vll = vl
349  vuu = vu
350  ELSE
351  vll = zero
352  vuu = zero
353  END IF
354  anrm = dlansp( 'M', uplo, n, ap, work )
355  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
356  iscale = 1
357  sigma = rmin / anrm
358  ELSE IF( anrm.GT.rmax ) THEN
359  iscale = 1
360  sigma = rmax / anrm
361  END IF
362  IF( iscale.EQ.1 ) THEN
363  CALL dscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
364  IF( abstol.GT.0 )
365  $ abstll = abstol*sigma
366  IF( valeig ) THEN
367  vll = vl*sigma
368  vuu = vu*sigma
369  END IF
370  END IF
371 *
372 * Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
373 *
374  indtau = 1
375  inde = indtau + n
376  indd = inde + n
377  indwrk = indd + n
378  CALL dsptrd( uplo, n, ap, work( indd ), work( inde ),
379  $ work( indtau ), iinfo )
380 *
381 * If all eigenvalues are desired and ABSTOL is less than or equal
382 * to zero, then call DSTERF or DOPGTR and SSTEQR. If this fails
383 * for some eigenvalue, then try DSTEBZ.
384 *
385  test = .false.
386  IF (indeig) THEN
387  IF (il.EQ.1 .AND. iu.EQ.n) THEN
388  test = .true.
389  END IF
390  END IF
391  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
392  CALL dcopy( n, work( indd ), 1, w, 1 )
393  indee = indwrk + 2*n
394  IF( .NOT.wantz ) THEN
395  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
396  CALL dsterf( n, w, work( indee ), info )
397  ELSE
398  CALL dopgtr( uplo, n, ap, work( indtau ), z, ldz,
399  $ work( indwrk ), iinfo )
400  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
401  CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
402  $ work( indwrk ), info )
403  IF( info.EQ.0 ) THEN
404  DO 10 i = 1, n
405  ifail( i ) = 0
406  10 continue
407  END IF
408  END IF
409  IF( info.EQ.0 ) THEN
410  m = n
411  go to 20
412  END IF
413  info = 0
414  END IF
415 *
416 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
417 *
418  IF( wantz ) THEN
419  order = 'B'
420  ELSE
421  order = 'E'
422  END IF
423  indibl = 1
424  indisp = indibl + n
425  indiwo = indisp + n
426  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
427  $ work( indd ), work( inde ), m, nsplit, w,
428  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
429  $ iwork( indiwo ), info )
430 *
431  IF( wantz ) THEN
432  CALL dstein( n, work( indd ), work( inde ), m, w,
433  $ iwork( indibl ), iwork( indisp ), z, ldz,
434  $ work( indwrk ), iwork( indiwo ), ifail, info )
435 *
436 * Apply orthogonal matrix used in reduction to tridiagonal
437 * form to eigenvectors returned by DSTEIN.
438 *
439  CALL dopmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
440  $ work( indwrk ), iinfo )
441  END IF
442 *
443 * If matrix was scaled, then rescale eigenvalues appropriately.
444 *
445  20 continue
446  IF( iscale.EQ.1 ) THEN
447  IF( info.EQ.0 ) THEN
448  imax = m
449  ELSE
450  imax = info - 1
451  END IF
452  CALL dscal( imax, one / sigma, w, 1 )
453  END IF
454 *
455 * If eigenvalues are not in order, then sort them, along with
456 * eigenvectors.
457 *
458  IF( wantz ) THEN
459  DO 40 j = 1, m - 1
460  i = 0
461  tmp1 = w( j )
462  DO 30 jj = j + 1, m
463  IF( w( jj ).LT.tmp1 ) THEN
464  i = jj
465  tmp1 = w( jj )
466  END IF
467  30 continue
468 *
469  IF( i.NE.0 ) THEN
470  itmp1 = iwork( indibl+i-1 )
471  w( i ) = w( j )
472  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
473  w( j ) = tmp1
474  iwork( indibl+j-1 ) = itmp1
475  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
476  IF( info.NE.0 ) THEN
477  itmp1 = ifail( i )
478  ifail( i ) = ifail( j )
479  ifail( j ) = itmp1
480  END IF
481  END IF
482  40 continue
483  END IF
484 *
485  return
486 *
487 * End of DSPEVX
488 *
489  END