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