1CSTEMR(1) LAPACK computational routine (version 3.2) CSTEMR(1)
2
6 CSTEMR - computes selected eigenvalues and, optionally, eigenvectors of
7 a real symmetric tridiagonal matrix T
8
10 SUBROUTINE CSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ,
11 NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK,
12 INFO )
13 IMPLICIT NONE
15 CHARACTER JOBZ, RANGE
17 LOGICAL TRYRAC
19 INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
21 REAL VL, VU
23 INTEGER ISUPPZ( * ), IWORK( * )
25 REAL D( * ), E( * ), W( * ), WORK( * )
27 COMPLEX Z( LDZ, * )
29
31 CSTEMR computes selected eigenvalues and, optionally, eigenvectors of a
32 real symmetric tridiagonal matrix T. Any such unreduced matrix has a
33 well defined set of pairwise different real eigenvalues, the corre‐
34 sponding real eigenvectors are pairwise orthogonal.
35 The spectrum may be computed either completely or partially by specify‐
36 ing either an interval (VL,VU] or a range of indices IL:IU for the
37 desired eigenvalues.
38 Depending on the number of desired eigenvalues, these are computed
39 either by bisection or the dqds algorithm. Numerically orthogonal
40 eigenvectors are computed by the use of various suitable L D L^T fac‐
41 torizations near clusters of close eigenvalues (referred to as RRRs,
42 Relatively Robust Representations). An informal sketch of the algorithm
43 follows. For each unreduced block (submatrix) of T,
44 (a) Compute T - sigma I = L D L^T, so that L and D
45 define all the wanted eigenvalues to high relative accuracy.
46 This means that small relative changes in the entries of D and L
47 cause only small relative changes in the eigenvalues and
48 eigenvectors. The standard (unfactored) representation of the
49 tridiagonal matrix T does not have this property in general.
50 (b) Compute the eigenvalues to suitable accuracy.
51 If the eigenvectors are desired, the algorithm attains full
52 accuracy of the computed eigenvalues only right before
53 the corresponding vectors have to be computed, see steps c) and
54 d).
55 (c) For each cluster of close eigenvalues, select a new
56 shift close to the cluster, find a new factorization, and refine
57 the shifted eigenvalues to suitable accuracy.
58 (d) For each eigenvalue with a large enough relative separation com‐
59 pute
60 the corresponding eigenvector by forming a rank revealing
61 twisted
62 factorization. Go back to (c) for any clusters that remain. For
63 more details, see:
64 - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representa‐
65 tions
66 to compute orthogonal eigenvectors of symmetric tridiagonal matri‐
67 ces,"
68 Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
69 - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
70 Relative Gaps," SIAM Journal on Matrix Analysis and Applications,
71 Vol. 25,
72 2004. Also LAPACK Working Note 154.
73 - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
74 tridiagonal eigenvalue/eigenvector problem",
75 Computer Science Division Technical Report No. UCB/CSD-97-971,
76 UC Berkeley, May 1997.
77 Further Details
78 floating-point standard in their handling of infinities and NaNs. This
79 permits the use of efficient inner loops avoiding a check for zero
80 divisors.
81 2. LAPACK routines can be used to reduce a complex Hermitean matrix to
82 real symmetric tridiagonal form.
83 (Any complex Hermitean tridiagonal matrix has real values on its diago‐
84 nal and potentially complex numbers on its off-diagonals. By applying a
85 similarity transform with an appropriate diagonal matrix
86 diag(1,e^{i hy_1},
87 ... , e^{i hy_{n-1}}),
88 the complex Hermitean matrix can be transformed into a real symmetric
89 matrix and complex arithmetic can be entirely avoided.)
90 While the eigenvectors of the real symmetric tridiagonal matrix are
91 real, the eigenvectors of original complex Hermitean matrix have com‐
92 plex entries in general.
93 Since LAPACK drivers overwrite the matrix data with the eigenvectors,
94 CSTEMR accepts complex workspace to facilitate interoperability with
95 CUNMTR or CUPMTR.
96
98 JOBZ (input) CHARACTER*1
99 = 'N': Compute eigenvalues only;
100 = 'V': Compute eigenvalues and eigenvectors.
101 RANGE (input) CHARACTER*1
103 = 'A': all eigenvalues will be found.
104 = 'V': all eigenvalues in the half-open interval (VL,VU] will
105 be found. = 'I': the IL-th through IU-th eigenvalues will be
106 found.
107 N (input) INTEGER
109 The order of the matrix. N >= 0.
110 D (input/output) REAL array, dimension (N)
112 On entry, the N diagonal elements of the tridiagonal matrix T.
113 On exit, D is overwritten.
114 E (input/output) REAL array, dimension (N)
116 On entry, the (N-1) subdiagonal elements of the tridiagonal
117 matrix T in elements 1 to N-1 of E. E(N) need not be set on
118 input, but is used internally as workspace. On exit, E is
119 overwritten.
120 VL (input) REAL
122 VU (input) REAL If RANGE='V', the lower and upper bounds
123 of the interval to be searched for eigenvalues. VL < VU. Not
124 referenced if RANGE = 'A' or 'I'.
125 IL (input) INTEGER
127 IU (input) INTEGER If RANGE='I', the indices (in ascending
128 order) of the smallest and largest eigenvalues to be returned.
129 1 <= IL <= IU <= N, if N > 0. Not referenced if RANGE = 'A' or
130 'V'.
131 M (output) INTEGER
133 The total number of eigenvalues found. 0 <= M <= N. If RANGE
134 = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
135 W (output) REAL array, dimension (N)
137 The first M elements contain the selected eigenvalues in
138 ascending order.
139 Z (output) COMPLEX array, dimension (LDZ, max(1,M) )
141 If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
142 contain the orthonormal eigenvectors of the matrix T corre‐
143 sponding to the selected eigenvalues, with the i-th column of Z
144 holding the eigenvector associated with W(i). If JOBZ = 'N',
145 then Z is not referenced. Note: the user must ensure that at
146 least max(1,M) columns are supplied in the array Z; if RANGE =
147 'V', the exact value of M is not known in advance and can be
148 computed with a workspace query by setting NZC = -1, see below.
149 LDZ (input) INTEGER
151 The leading dimension of the array Z. LDZ >= 1, and if JOBZ =
152 'V', then LDZ >= max(1,N).
153 NZC (input) INTEGER
155 The number of eigenvectors to be held in the array Z. If RANGE
156 = 'A', then NZC >= max(1,N). If RANGE = 'V', then NZC >= the
157 number of eigenvalues in (VL,VU]. If RANGE = 'I', then NZC >=
158 IU-IL+1. If NZC = -1, then a workspace query is assumed; the
159 routine calculates the number of columns of the array Z that
160 are needed to hold the eigenvectors. This value is returned as
161 the first entry of the Z array, and no error message related to
162 NZC is issued by XERBLA.
163 ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
165 The support of the eigenvectors in Z, i.e., the indices indi‐
166 cating the nonzero elements in Z. The i-th computed eigenvector
167 is nonzero only in elements ISUPPZ( 2*i-1 ) through ISUPPZ( 2*i
168 ). This is relevant in the case when the matrix is split.
169 ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
170 TRYRAC (input/output) LOGICAL
172 If TRYRAC.EQ..TRUE., indicates that the code should check
173 whether the tridiagonal matrix defines its eigenvalues to high
174 relative accuracy. If so, the code uses relative-accuracy pre‐
175 serving algorithms that might be (a bit) slower depending on
176 the matrix. If the matrix does not define its eigenvalues to
177 high relative accuracy, the code can uses possibly faster algo‐
178 rithms. If TRYRAC.EQ..FALSE., the code is not required to
179 guarantee relatively accurate eigenvalues and can use the
180 fastest possible techniques. On exit, a .TRUE. TRYRAC will be
181 set to .FALSE. if the matrix does not define its eigenvalues to
182 high relative accuracy.
183 WORK (workspace/output) REAL array, dimension (LWORK)
185 On exit, if INFO = 0, WORK(1) returns the optimal (and minimal)
186 LWORK.
187 LWORK (input) INTEGER
189 The dimension of the array WORK. LWORK >= max(1,18*N) if JOBZ =
190 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. If LWORK = -1,
191 then a workspace query is assumed; the routine only calculates
192 the optimal size of the WORK array, returns this value as the
193 first entry of the WORK array, and no error message related to
194 LWORK is issued by XERBLA.
195 IWORK (workspace/output) INTEGER array, dimension (LIWORK)
197 On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
198 LIWORK (input) INTEGER
200 The dimension of the array IWORK. LIWORK >= max(1,10*N) if the
201 eigenvectors are desired, and LIWORK >= max(1,8*N) if only the
202 eigenvalues are to be computed. If LIWORK = -1, then a
203 workspace query is assumed; the routine only calculates the
204 optimal size of the IWORK array, returns this value as the
205 first entry of the IWORK array, and no error message related to
206 LIWORK is issued by XERBLA.
207 INFO (output) INTEGER
209 On exit, INFO = 0: successful exit
210 < 0: if INFO = -i, the i-th argument had an illegal value
211 > 0: if INFO = 1X, internal error in SLARRE, if INFO = 2X,
212 internal error in CLARRV. Here, the digit X = ABS( IINFO ) <
213 10, where IINFO is the nonzero error code returned by SLARRE or
214 CLARRV, respectively.
215
217 Based on contributions by
218 Beresford Parlett, University of California, Berkeley, USA
219 Jim Demmel, University of California, Berkeley, USA
220 Inderjit Dhillon, University of Texas, Austin, USA
221 Osni Marques, LBNL/NERSC, USA
222 Christof Voemel, University of California, Berkeley, USA
223 LAPACK computational routine (verNsoivoenmb3e.r2)2008 CSTEMR(1)