1CSTEMR ‐ selected eigenvalues and, optionally, eigenvectors of a
2real symmetric tridiagonal matrix T SUBROUTINE CSTEMR( 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 REAL VL, VU
10 INTEGER ISUPPZ( * ), IWORK( * )
11 REAL D( * ), E( * ), W( * ), WORK( * )
12 COMPLEX Z( LDZ, * ) CSTEMR computes selected eigenvalues and,
13optionally, eigenvectors of a real symmetric tridiagonal matrix
14T. Any such unreduced matrix has a well defined set of pairwise
15different real eigenvalues, the corresponding real eigenvectors
16are 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.CSTEMR 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, CSTEMR accepts complex workspace to facilitate interoper‐
96ability with CUNMTR or CUPMTR.
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) REAL array, dimension (N) On entry, the N diagonal
106elements of the tridiagonal matrix T. On exit, D is overwritten.
107E (input/output) REAL array, dimension (N) On entry, the
108(N‐1) subdiagonal elements of the tridiagonal matrix T in ele‐
109ments 1 to N‐1 of E. E(N) need not be set on input, but is used
110internally as workspace. On exit, E is overwritten. VL
111(input) REAL VU (input) REAL If RANGE='V', the lower and up‐
112per bounds of the interval to be searched for eigenvalues. VL <
113VU. Not referenced if RANGE = 'A' or 'I'. IL (input) INTE‐
114GER IU (input) INTEGER If RANGE='I', the indices (in ascend‐
115ing order) of the smallest and largest eigenvalues to be re‐
116turned. 1 <= IL <= IU <= N, if N > 0. Not referenced if RANGE =
117'A' or 'V'. M (output) INTEGER The total number of eigen‐
118values found. 0 <= M <= N. If RANGE = 'A', M = N, and if RANGE
119= 'I', M = IU‐IL+1. W (output) REAL array, dimension (N)
120The first M elements contain the selected eigenvalues in ascend‐
121ing order. Z (output) COMPLEX array, dimension (LDZ,
122max(1,M) ) If JOBZ = 'V', and if INFO = 0, then the first M col‐
123umns of Z contain the orthonormal eigenvectors of the matrix T
124corresponding to the selected eigenvalues, with the i‐th column
125of Z holding the eigenvector associated with W(i). If JOBZ =
126'N', then Z is not referenced. Note: the user must ensure that
127at least max(1,M) columns are supplied in the array Z; if RANGE =
128'V', the exact value of M is not known in advance and can be com‐
129puted with a workspace query by setting NZC = ‐1, see below. LDZ
130(input) INTEGER The leading dimension of the array Z. LDZ >= 1,
131and if JOBZ = 'V', then LDZ >= max(1,N). NZC (input) INTEGER
132The number of eigenvectors to be held in the array Z. If RANGE =
133'A', then NZC >= max(1,N). If RANGE = 'V', then NZC >= the num‐
134ber of eigenvalues in (VL,VU]. If RANGE = 'I', then NZC >= IU‐
135IL+1. If NZC = ‐1, then a workspace query is assumed; the rou‐
136tine calculates the number of columns of the array Z that are
137needed to hold the eigenvectors. This value is returned as the
138first entry of the Z array, and no error message related to NZC
139is issued by XERBLA. ISUPPZ (output) INTEGER ARRAY, dimension (
1402*max(1,M) ) The support of the eigenvectors in Z, i.e., the in‐
141dices indicating the nonzero elements in Z. The i‐th computed
142eigenvector is nonzero only in elements ISUPPZ( 2*i‐1 ) through
143ISUPPZ( 2*i ). This is relevant in the case when the matrix is
144split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
145TRYRAC (input/output) LOGICAL If TRYRAC.EQ..TRUE., indicates
146that the code should check whether the tridiagonal matrix defines
147its eigenvalues to high relative accuracy. If so, the code uses
148relative‐accuracy preserving algorithms that might be (a bit)
149slower depending on the matrix. If the matrix does not define
150its eigenvalues to high relative accuracy, the code can uses pos‐
151sibly faster algorithms. If TRYRAC.EQ..FALSE., the code is not
152required to guarantee relatively accurate eigenvalues and can use
153the fastest possible techniques. On exit, a .TRUE. TRYRAC will
154be set to .FALSE. if the matrix does not define its eigenvalues
155to high relative accuracy. WORK (workspace/output) REAL ar‐
156ray, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the
157optimal (and minimal) LWORK. LWORK (input) INTEGER The dimen‐
158sion of the array WORK. LWORK >= max(1,18*N) if JOBZ = 'V', and
159LWORK >= max(1,12*N) if JOBZ = 'N'. If LWORK = ‐1, then a
160workspace query is assumed; the routine only calculates the opti‐
161mal size of the WORK array, returns this value as the first entry
162of the WORK array, and no error message related to LWORK is is‐
163sued by XERBLA. IWORK (workspace/output) INTEGER array, dimen‐
164sion (LIWORK) On exit, if INFO = 0, IWORK(1) returns the optimal
165LIWORK. LIWORK (input) INTEGER The dimension of the array
166IWORK. LIWORK >= max(1,10*N) if the eigenvectors are desired,
167and LIWORK >= max(1,8*N) if only the eigenvalues are to be com‐
168puted. If LIWORK = ‐1, then a workspace query is assumed; the
169routine only calculates the optimal size of the IWORK array, re‐
170turns this value as the first entry of the IWORK array, and no
171error message related to LIWORK is issued by XERBLA. INFO
172(output) INTEGER On exit, INFO = 0: successful exit
173< 0: if INFO = ‐i, the i‐th argument had an illegal value
174> 0: if INFO = 1X, internal error in SLARRE, if INFO = 2X, in‐
175ternal error in CLARRV. Here, the digit X = ABS( IINFO ) < 10,
176where IINFO is the nonzero error code returned by SLARRE or CLAR‐
177RV, respectively. Based on contributions by
178 Beresford Parlett, University of California, Berkeley, USA
179 Jim Demmel, University of California, Berkeley, USA
180 Inderjit Dhillon, University of Texas, Austin, USA
181 Osni Marques, LBNL/NERSC, USA
182 Christof Voemel, University of California, Berkeley, USA
183