1DSTEMR(1) LAPACK computational routine (version 3.1) DSTEMR(1)
2
6 DSTEMR - selected eigenvalues and, optionally, eigenvectors of a real
7 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‐
37 ing either an interval (VL,VU] or a range of indices IL:IU for the
38 desired eigenvalues.
39 Depending on the number of desired eigenvalues, these are computed
41 either by bisection or the dqds algorithm. Numerically orthogonal
42 eigenvectors are computed by the use of various suitable L D L^T fac‐
43 torizations near clusters of close eigenvalues (referred to as RRRs,
44 Relatively Robust Representations). An informal sketch of the algorithm
45 follows.
46 For each unreduced block (submatrix) of T,
48 (a) Compute T - sigma I = L D L^T, so that L and D
49 define all the wanted eigenvalues to high relative accuracy.
50 This means that small relative changes in the entries of D and L
51 cause only small relative changes in the eigenvalues and
52 eigenvectors. The standard (unfactored) representation of the
53 tridiagonal matrix T does not have this property in general.
54 (b) Compute the eigenvalues to suitable accuracy.
55 If the eigenvectors are desired, the algorithm attains full
56 accuracy of the computed eigenvalues only right before
57 the corresponding vectors have to be computed, see steps c) and
58 d).
59 (c) For each cluster of close eigenvalues, select a new
60 shift close to the cluster, find a new factorization, and refine
61 the shifted eigenvalues to suitable accuracy.
62 (d) For each eigenvalue with a large enough relative separation com‐
63 pute
64 the corresponding eigenvector by forming a rank revealing
65 twisted
66 factorization. Go back to (c) for any clusters that remain.
67 For more details, see:
69 - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representa‐
70 tions
71 to compute orthogonal eigenvectors of symmetric tridiagonal matri‐
72 ces,"
73 Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
74 - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
75 Relative Gaps," SIAM Journal on Matrix Analysis and Applications,
76 Vol. 25,
77 2004. Also LAPACK Working Note 154.
78 - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
79 tridiagonal eigenvalue/eigenvector problem",
80 Computer Science Division Technical Report No. UCB/CSD-97-971,
81 UC Berkeley, May 1997.
82 Notes:
84 1.DSTEMR works only on machines which follow IEEE-754
85 floating-point standard in their handling of infinities and NaNs. This
86 permits the use of efficient inner loops avoiding a check for zero
87 divisors.
88
91 JOBZ (input) CHARACTER*1
92 = 'N': Compute eigenvalues only;
93 = 'V': Compute eigenvalues and eigenvectors.
94 RANGE (input) CHARACTER*1
96 = 'A': all eigenvalues will be found.
97 = 'V': all eigenvalues in the half-open interval (VL,VU] will
98 be found. = 'I': the IL-th through IU-th eigenvalues will be
99 found.
100 N (input) INTEGER
102 The order of the matrix. N >= 0.
103 D (input/output) DOUBLE PRECISION array, dimension (N)
105 On entry, the N diagonal elements of the tridiagonal matrix T.
106 On exit, D is overwritten.
107 E (input/output) DOUBLE PRECISION array, dimension (N)
109 On entry, the (N-1) subdiagonal elements of the tridiagonal
110 matrix T in elements 1 to N-1 of E. E(N) need not be set on
111 input, but is used internally as workspace. On exit, E is
112 overwritten.
113 VL (input) DOUBLE PRECISION
115 VU (input) DOUBLE PRECISION If RANGE='V', the lower and
116 upper bounds of the interval to be searched for eigenvalues. VL
117 < VU. Not referenced if RANGE = 'A' or 'I'.
118 IL (input) INTEGER
120 IU (input) INTEGER If RANGE='I', the indices (in ascending
121 order) of the smallest and largest eigenvalues to be returned.
122 1 <= IL <= IU <= N, if N > 0. Not referenced if RANGE = 'A' or
123 'V'.
124 M (output) INTEGER
126 The total number of eigenvalues found. 0 <= M <= N. If RANGE
127 = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
128 W (output) DOUBLE PRECISION array, dimension (N)
130 The first M elements contain the selected eigenvalues in
131 ascending order.
132 Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
134 If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
135 contain the orthonormal eigenvectors of the matrix T corre‐
136 sponding to the selected eigenvalues, with the i-th column of Z
137 holding the eigenvector associated with W(i). If JOBZ = 'N',
138 then Z is not referenced. Note: the user must ensure that at
139 least max(1,M) columns are supplied in the array Z; if RANGE =
140 'V', the exact value of M is not known in advance and can be
141 computed with a workspace query by setting NZC = -1, see below.
142 LDZ (input) INTEGER
144 The leading dimension of the array Z. LDZ >= 1, and if JOBZ =
145 'V', then LDZ >= max(1,N).
146 NZC (input) INTEGER
148 The number of eigenvectors to be held in the array Z. If RANGE
149 = 'A', then NZC >= max(1,N). If RANGE = 'V', then NZC >= the
150 number of eigenvalues in (VL,VU]. If RANGE = 'I', then NZC >=
151 IU-IL+1. If NZC = -1, then a workspace query is assumed; the
152 routine calculates the number of columns of the array Z that
153 are needed to hold the eigenvectors. This value is returned as
154 the first entry of the Z array, and no error message related to
155 NZC is issued by XERBLA.
156 ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
158 The support of the eigenvectors in Z, i.e., the indices indi‐
159 cating the nonzero elements in Z. The i-th computed eigenvector
160 is nonzero only in elements ISUPPZ( 2*i-1 ) through ISUPPZ( 2*i
161 ). This is relevant in the case when the matrix is split.
162 ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
163 TRYRAC (input/output) LOGICAL
165 If TRYRAC.EQ..TRUE., indicates that the code should check
166 whether the tridiagonal matrix defines its eigenvalues to high
167 relative accuracy. If so, the code uses relative-accuracy pre‐
168 serving algorithms that might be (a bit) slower depending on
169 the matrix. If the matrix does not define its eigenvalues to
170 high relative accuracy, the code can uses possibly faster algo‐
171 rithms. If TRYRAC.EQ..FALSE., the code is not required to
172 guarantee relatively accurate eigenvalues and can use the
173 fastest possible techniques. On exit, a .TRUE. TRYRAC will be
174 set to .FALSE. if the matrix does not define its eigenvalues to
175 high relative accuracy.
176 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
178 On exit, if INFO = 0, WORK(1) returns the optimal (and minimal)
179 LWORK.
180 LWORK (input) INTEGER
182 The dimension of the array WORK. LWORK >= max(1,18*N) if JOBZ =
183 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. If LWORK = -1,
184 then a workspace query is assumed; the routine only calculates
185 the optimal size of the WORK array, returns this value as the
186 first entry of the WORK array, and no error message related to
187 LWORK is issued by XERBLA.
188 IWORK (workspace/output) INTEGER array, dimension (LIWORK)
190 On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
191 LIWORK (input) INTEGER
193 The dimension of the array IWORK. LIWORK >= max(1,10*N) if the
194 eigenvectors are desired, and LIWORK >= max(1,8*N) if only the
195 eigenvalues are to be computed. If LIWORK = -1, then a
196 workspace query is assumed; the routine only calculates the
197 optimal size of the IWORK array, returns this value as the
198 first entry of the IWORK array, and no error message related to
199 LIWORK is issued by XERBLA.
200 INFO (output) INTEGER
202 On exit, INFO = 0: successful exit
203 < 0: if INFO = -i, the i-th argument had an illegal value
204 > 0: if INFO = 1X, internal error in DLARRE, if INFO = 2X,
205 internal error in DLARRV. Here, the digit X = ABS( IINFO ) <
206 10, where IINFO is the nonzero error code returned by DLARRE or
207 DLARRV, respectively.
208
210 Based on contributions by
211 Beresford Parlett, University of California, Berkeley, USA
212 Jim Demmel, University of California, Berkeley, USA
213 Inderjit Dhillon, University of Texas, Austin, USA
214 Osni Marques, LBNL/NERSC, USA
215 Christof Voemel, University of California, Berkeley, USA
216 LAPACK computational routine (verNsoivoenmb3e.r1)2006 DSTEMR(1)