1ZSTEMR ‐ selected eigenvalues and, optionally, eigenvectors of a
2real symmetric tridiagonal matrix T SUBROUTINE ZSTEMR( JOBZ,
3RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ,
4TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO )
5 IMPLICIT NONE
6 CHARACTER JOBZ, RANGE
7 LOGICAL TRYRAC
8 INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
9 DOUBLE PRECISION VL, VU
10 INTEGER ISUPPZ( * ), IWORK( * )
11 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
12 COMPLEX*16 Z( LDZ, * ) ZSTEMR computes selected eigenvalues
13and, optionally, eigenvectors of a real symmetric tridiagonal ma‐
14trix T. Any such unreduced matrix has a well defined set of pair‐
15wise different real eigenvalues, the corresponding real eigenvec‐
16tors are pairwise orthogonal.
17The spectrum may be computed either completely or partially by
19specifying either an interval (VL,VU] or a range of indices IL:IU
20for the desired eigenvalues.
21Depending on the number of desired eigenvalues, these are comput‐
23ed either by bisection or the dqds algorithm. Numerically orthog‐
24onal eigenvectors are computed by the use of various suitable L D
25L^T factorizations near clusters of close eigenvalues (referred
26to as RRRs, Relatively Robust Representations). An informal
27sketch of the algorithm follows.
28For each unreduced block (submatrix) of T,
30 (a) Compute T ‐ sigma I = L D L^T, so that L and D
31 define all the wanted eigenvalues to high relative accura‐
32cy.
33 This means that small relative changes in the entries of D
34and L
35 cause only small relative changes in the eigenvalues and
36 eigenvectors. The standard (unfactored) representation of
37the
38 tridiagonal matrix T does not have this property in gener‐
39al.
40 (b) Compute the eigenvalues to suitable accuracy.
41 If the eigenvectors are desired, the algorithm attains
42full
43 accuracy of the computed eigenvalues only right before
44 the corresponding vectors have to be computed, see steps
45c) and d).
46 (c) For each cluster of close eigenvalues, select a new
47 shift close to the cluster, find a new factorization, and
48refine
49 the shifted eigenvalues to suitable accuracy.
50 (d) For each eigenvalue with a large enough relative separa‐
51tion compute
52 the corresponding eigenvector by forming a rank revealing
53twisted
54 factorization. Go back to (c) for any clusters that re‐
55main.
56For more details, see:
58‐ Inderjit S. Dhillon and Beresford N. Parlett: "Multiple repre‐
59sentations
60 to compute orthogonal eigenvectors of symmetric tridiagonal ma‐
61trices,"
62 Linear Algebra and its Applications, 387(1), pp. 1‐28, August
632004. ‐ Inderjit Dhillon and Beresford Parlett: "Orthogonal
64Eigenvectors and
65 Relative Gaps," SIAM Journal on Matrix Analysis and Applica‐
66tions, Vol. 25,
67 2004. Also LAPACK Working Note 154.
68‐ Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
69 tridiagonal eigenvalue/eigenvector problem",
70 Computer Science Division Technical Report No. UCB/CSD‐97‐971,
71 UC Berkeley, May 1997.
72Notes:
741.ZSTEMR works only on machines which follow IEEE‐754
75floating‐point standard in their handling of infinities and NaNs.
76This permits the use of efficient inner loops avoiding a check
77for zero divisors.
782. LAPACK routines can be used to reduce a complex Hermitean ma‐
80trix to real symmetric tridiagonal form.
81(Any complex Hermitean tridiagonal matrix has real values on its
83diagonal and potentially complex numbers on its off‐diagonals. By
84applying a similarity transform with an appropriate diagonal ma‐
85trix
86diag(1,e^{i hy_1},
87... , e^{i hy_{n‐1}}),
88the complex Hermitean matrix can be transformed into a real sym‐
89metric matrix and complex arithmetic can be entirely avoided.)
90While the eigenvectors of the real symmetric tridiagonal matrix
92are real, the eigenvectors of original complex Hermitean matrix
93have complex entries in general.
94Since LAPACK drivers overwrite the matrix data with the eigenvec‐
95tors, ZSTEMR accepts complex workspace to facilitate interoper‐
96ability with ZUNMTR or ZUPMTR.
97JOBZ (input) CHARACTER*1 = 'N': Compute eigenvalues only;
99= 'V': Compute eigenvalues and eigenvectors. RANGE (input)
100CHARACTER*1
101= 'A': all eigenvalues will be found.
102= 'V': all eigenvalues in the half‐open interval (VL,VU] will be
103found. = 'I': the IL‐th through IU‐th eigenvalues will be found.
104N (input) INTEGER The order of the matrix. N >= 0. D
105(input/output) DOUBLE PRECISION array, dimension (N) On entry,
106the N diagonal elements of the tridiagonal matrix T. On exit, D
107is overwritten. E (input/output) DOUBLE PRECISION array,
108dimension (N) On entry, the (N‐1) subdiagonal elements of the
109tridiagonal matrix T in elements 1 to N‐1 of E. E(N) need not be
110set on input, but is used internally as workspace. On exit, E is
111overwritten. VL (input) DOUBLE PRECISION VU (input)
112DOUBLE PRECISION If RANGE='V', the lower and upper bounds of the
113interval to be searched for eigenvalues. VL < VU. Not referenced
114if RANGE = 'A' or 'I'. IL (input) INTEGER IU (input)
115INTEGER If RANGE='I', the indices (in ascending order) of the
116smallest and largest eigenvalues to be returned. 1 <= IL <= IU
117<= N, if N > 0. Not referenced if RANGE = 'A' or 'V'. M
118(output) INTEGER The total number of eigenvalues found. 0 <= M
119<= N. If RANGE = 'A', M = N, and if RANGE = 'I', M = IU‐IL+1. W
120(output) DOUBLE PRECISION array, dimension (N) The first M ele‐
121ments contain the selected eigenvalues in ascending order. Z
122(output) COMPLEX*16 array, dimension (LDZ, max(1,M) ) If JOBZ =
123'V', and if INFO = 0, then the first M columns of Z contain the
124orthonormal eigenvectors of the matrix T corresponding to the se‐
125lected eigenvalues, with the i‐th column of Z holding the eigen‐
126vector associated with W(i). If JOBZ = 'N', then Z is not refer‐
127enced. Note: the user must ensure that at least max(1,M) columns
128are supplied in the array Z; if RANGE = 'V', the exact value of M
129is not known in advance and can be computed with a workspace
130query by setting NZC = ‐1, see below. LDZ (input) INTEGER
131The leading dimension of the array Z. LDZ >= 1, and if JOBZ =
132'V', then LDZ >= max(1,N). NZC (input) INTEGER The number of
133eigenvectors to be held in the array Z. If RANGE = 'A', then NZC
134>= max(1,N). If RANGE = 'V', then NZC >= the number of eigenval‐
135ues in (VL,VU]. If RANGE = 'I', then NZC >= IU‐IL+1. If NZC =
136‐1, then a workspace query is assumed; the routine calculates the
137number of columns of the array Z that are needed to hold the
138eigenvectors. This value is returned as the first entry of the Z
139array, and no error message related to NZC is issued by XERBLA.
140ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) The sup‐
141port of the eigenvectors in Z, i.e., the indices indicating the
142nonzero elements in Z. The i‐th computed eigenvector is nonzero
143only in elements ISUPPZ( 2*i‐1 ) through ISUPPZ( 2*i ). This is
144relevant in the case when the matrix is split. ISUPPZ is only ac‐
145cessed when JOBZ is 'V' and N > 0. TRYRAC (input/output) LOGI‐
146CAL If TRYRAC.EQ..TRUE., indicates that the code should check
147whether the tridiagonal matrix defines its eigenvalues to high
148relative accuracy. If so, the code uses relative‐accuracy pre‐
149serving algorithms that might be (a bit) slower depending on the
150matrix. If the matrix does not define its eigenvalues to high
151relative accuracy, the code can uses possibly faster algorithms.
152If TRYRAC.EQ..FALSE., the code is not required to guarantee rela‐
153tively accurate eigenvalues and can use the fastest possible
154techniques. On exit, a .TRUE. TRYRAC will be set to .FALSE. if
155the matrix does not define its eigenvalues to high relative accu‐
156racy. WORK (workspace/output) DOUBLE PRECISION array, dimen‐
157sion (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal
158(and minimal) LWORK. LWORK (input) INTEGER The dimension of
159the array WORK. LWORK >= max(1,18*N) if JOBZ = 'V', and LWORK >=
160max(1,12*N) if JOBZ = 'N'. If LWORK = ‐1, then a workspace query
161is assumed; the routine only calculates the optimal size of the
162WORK array, returns this value as the first entry of the WORK ar‐
163ray, and no error message related to LWORK is issued by XERBLA.
164IWORK (workspace/output) INTEGER array, dimension (LIWORK) On
165exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. LIWORK
166(input) INTEGER The dimension of the array IWORK. LIWORK >=
167max(1,10*N) if the eigenvectors are desired, and LIWORK >=
168max(1,8*N) if only the eigenvalues are to be computed. If LIWORK
169= ‐1, then a workspace query is assumed; the routine only calcu‐
170lates the optimal size of the IWORK array, returns this value as
171the first entry of the IWORK array, and no error message related
172to LIWORK is issued by XERBLA. INFO (output) INTEGER On exit,
173INFO = 0: successful exit
174< 0: if INFO = ‐i, the i‐th argument had an illegal value
175> 0: if INFO = 1X, internal error in DLARRE, if INFO = 2X, in‐
176ternal error in ZLARRV. Here, the digit X = ABS( IINFO ) < 10,
177where IINFO is the nonzero error code returned by DLARRE or ZLAR‐
178RV, respectively. Based on contributions by
179 Beresford Parlett, University of California, Berkeley, USA
180 Jim Demmel, University of California, Berkeley, USA
181 Inderjit Dhillon, University of Texas, Austin, USA
182 Osni Marques, LBNL/NERSC, USA
183 Christof Voemel, University of California, Berkeley, USA
184