1DSTEMR(1) LAPACK computational routine (version 3.2) DSTEMR(1)
2
6 DSTEMR - computes selected eigenvalues and, optionally, eigenvectors of
7 a real symmetric tridiagonal matrix T
8
10 SUBROUTINE DSTEMR( 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 DOUBLE PRECISION VL, VU
23 INTEGER ISUPPZ( * ), IWORK( * )
25 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
27 DOUBLE PRECISION Z( LDZ, * )
29
31 DSTEMR 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
83 JOBZ (input) CHARACTER*1
84 = 'N': Compute eigenvalues only;
85 = 'V': Compute eigenvalues and eigenvectors.
86 RANGE (input) CHARACTER*1
88 = 'A': all eigenvalues will be found.
89 = 'V': all eigenvalues in the half-open interval (VL,VU] will
90 be found. = 'I': the IL-th through IU-th eigenvalues will be
91 found.
92 N (input) INTEGER
94 The order of the matrix. N >= 0.
95 D (input/output) DOUBLE PRECISION array, dimension (N)
97 On entry, the N diagonal elements of the tridiagonal matrix T.
98 On exit, D is overwritten.
99 E (input/output) DOUBLE PRECISION array, dimension (N)
101 On entry, the (N-1) subdiagonal elements of the tridiagonal
102 matrix T in elements 1 to N-1 of E. E(N) need not be set on
103 input, but is used internally as workspace. On exit, E is
104 overwritten.
105 VL (input) DOUBLE PRECISION
107 VU (input) DOUBLE PRECISION If RANGE='V', the lower and
108 upper bounds of the interval to be searched for eigenvalues. VL
109 < VU. Not referenced if RANGE = 'A' or 'I'.
110 IL (input) INTEGER
112 IU (input) INTEGER If RANGE='I', the indices (in ascending
113 order) of the smallest and largest eigenvalues to be returned.
114 1 <= IL <= IU <= N, if N > 0. Not referenced if RANGE = 'A' or
115 'V'.
116 M (output) INTEGER
118 The total number of eigenvalues found. 0 <= M <= N. If RANGE
119 = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
120 W (output) DOUBLE PRECISION array, dimension (N)
122 The first M elements contain the selected eigenvalues in
123 ascending order.
124 Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
126 If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
127 contain the orthonormal eigenvectors of the matrix T corre‐
128 sponding to the selected eigenvalues, with the i-th column of Z
129 holding the eigenvector associated with W(i). If JOBZ = 'N',
130 then Z is not referenced. Note: the user must ensure that at
131 least max(1,M) columns are supplied in the array Z; if RANGE =
132 'V', the exact value of M is not known in advance and can be
133 computed with a workspace query by setting NZC = -1, see below.
134 LDZ (input) INTEGER
136 The leading dimension of the array Z. LDZ >= 1, and if JOBZ =
137 'V', then LDZ >= max(1,N).
138 NZC (input) INTEGER
140 The number of eigenvectors to be held in the array Z. If RANGE
141 = 'A', then NZC >= max(1,N). If RANGE = 'V', then NZC >= the
142 number of eigenvalues in (VL,VU]. If RANGE = 'I', then NZC >=
143 IU-IL+1. If NZC = -1, then a workspace query is assumed; the
144 routine calculates the number of columns of the array Z that
145 are needed to hold the eigenvectors. This value is returned as
146 the first entry of the Z array, and no error message related to
147 NZC is issued by XERBLA.
148 ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
150 The support of the eigenvectors in Z, i.e., the indices indi‐
151 cating the nonzero elements in Z. The i-th computed eigenvector
152 is nonzero only in elements ISUPPZ( 2*i-1 ) through ISUPPZ( 2*i
153 ). This is relevant in the case when the matrix is split.
154 ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
155 TRYRAC (input/output) LOGICAL
157 If TRYRAC.EQ..TRUE., indicates that the code should check
158 whether the tridiagonal matrix defines its eigenvalues to high
159 relative accuracy. If so, the code uses relative-accuracy pre‐
160 serving algorithms that might be (a bit) slower depending on
161 the matrix. If the matrix does not define its eigenvalues to
162 high relative accuracy, the code can uses possibly faster algo‐
163 rithms. If TRYRAC.EQ..FALSE., the code is not required to
164 guarantee relatively accurate eigenvalues and can use the
165 fastest possible techniques. On exit, a .TRUE. TRYRAC will be
166 set to .FALSE. if the matrix does not define its eigenvalues to
167 high relative accuracy.
168 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
170 On exit, if INFO = 0, WORK(1) returns the optimal (and minimal)
171 LWORK.
172 LWORK (input) INTEGER
174 The dimension of the array WORK. LWORK >= max(1,18*N) if JOBZ =
175 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. If LWORK = -1,
176 then a workspace query is assumed; the routine only calculates
177 the optimal size of the WORK array, returns this value as the
178 first entry of the WORK array, and no error message related to
179 LWORK is issued by XERBLA.
180 IWORK (workspace/output) INTEGER array, dimension (LIWORK)
182 On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
183 LIWORK (input) INTEGER
185 The dimension of the array IWORK. LIWORK >= max(1,10*N) if the
186 eigenvectors are desired, and LIWORK >= max(1,8*N) if only the
187 eigenvalues are to be computed. If LIWORK = -1, then a
188 workspace query is assumed; the routine only calculates the
189 optimal size of the IWORK array, returns this value as the
190 first entry of the IWORK array, and no error message related to
191 LIWORK is issued by XERBLA.
192 INFO (output) INTEGER
194 On exit, INFO = 0: successful exit
195 < 0: if INFO = -i, the i-th argument had an illegal value
196 > 0: if INFO = 1X, internal error in DLARRE, if INFO = 2X,
197 internal error in DLARRV. Here, the digit X = ABS( IINFO ) <
198 10, where IINFO is the nonzero error code returned by DLARRE or
199 DLARRV, respectively.
200
202 Based on contributions by
203 Beresford Parlett, University of California, Berkeley, USA
204 Jim Demmel, University of California, Berkeley, USA
205 Inderjit Dhillon, University of Texas, Austin, USA
206 Osni Marques, LBNL/NERSC, USA
207 Christof Voemel, University of California, Berkeley, USA
208 LAPACK computational routine (verNsoivoenmb3e.r2)2008 DSTEMR(1)