ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
dstegr2b.f
Go to the documentation of this file.
1  SUBROUTINE dstegr2b( JOBZ, N, D, E,
2  $ M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK,
3  $ LIWORK, DOL, DOU, NEEDIL, NEEDIU,
4  $ INDWLC, PIVMIN, SCALE, WL, WU,
5  $ VSTART, FINISH, MAXCLS,
6  $ NDEPTH, PARITY, ZOFFSET, INFO )
7 *
8 * -- ScaLAPACK auxiliary routine (version 2.0) --
9 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
10 * July 4, 2010
11 *
12  IMPLICIT NONE
13 *
14 * .. Scalar Arguments ..
15  CHARACTER JOBZ
16  INTEGER DOL, DOU, INDWLC, INFO, LDZ, LIWORK, LWORK, M,
17  $ maxcls, n, ndepth, needil, neediu, nzc, parity,
18  $ zoffset
19 
20  DOUBLE PRECISION PIVMIN, SCALE, WL, WU
21  LOGICAL VSTART, FINISH
22 
23 * ..
24 * .. Array Arguments ..
25  INTEGER ISUPPZ( * ), IWORK( * )
26  DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
27  DOUBLE PRECISION Z( LDZ, * )
28 * ..
29 *
30 * Purpose
31 * =======
32 *
33 * DSTEGR2B should only be called after a call to DSTEGR2A.
34 * From eigenvalues and initial representations computed by DSTEGR2A,
35 * DSTEGR2B computes the selected eigenvalues and eigenvectors
36 * of the real symmetric tridiagonal matrix in parallel
37 * on multiple processors. It is potentially invoked multiple times
38 * on a given processor because the locally relevant representation tree
39 * might depend on spectral information that is "owned" by other processors
40 * and might need to be communicated.
41 *
42 * Please note:
43 * 1. The calling sequence has two additional INTEGER parameters,
44 * DOL and DOU, that should satisfy M>=DOU>=DOL>=1.
45 * These parameters are only relevant for the case JOBZ = 'V'.
46 * DSTEGR2B ONLY computes the eigenVECTORS
47 * corresponding to eigenvalues DOL through DOU in W. (That is,
48 * instead of computing the eigenvectors belonging to W(1)
49 * through W(M), only the eigenvectors belonging to eigenvalues
50 * W(DOL) through W(DOU) are computed. In this case, only the
51 * eigenvalues DOL:DOU are guaranteed to be accurately refined
52 * to all figures by Rayleigh-Quotient iteration.
53 *
54 * 2. The additional arguments VSTART, FINISH, NDEPTH, PARITY, ZOFFSET
55 * are included as a thread-safe implementation equivalent to SAVE variables.
56 * These variables store details about the local representation tree which is
57 * computed layerwise. For scalability reasons, eigenvalues belonging to the
58 * locally relevant representation tree might be computed on other processors.
59 * These need to be communicated before the inspection of the RRRs can proceed
60 * on any given layer.
61 * Note that only when the variable FINISH is true, the computation has ended
62 * All eigenpairs between DOL and DOU have been computed. M is set = DOU - DOL + 1.
63 *
64 * 3. DSTEGR2B needs more workspace in Z than the sequential DSTEGR.
65 * It is used to store the conformal embedding of the local representation tree.
66 *
67 * Arguments
68 * =========
69 *
70 * JOBZ (input) CHARACTER*1
71 * = 'N': Compute eigenvalues only;
72 * = 'V': Compute eigenvalues and eigenvectors.
73 *
74 * N (input) INTEGER
75 * The order of the matrix. N >= 0.
76 *
77 * D (input/output) DOUBLE PRECISION array, dimension (N)
78 * On entry, the N diagonal elements of the tridiagonal matrix
79 * T. On exit, D is overwritten.
80 *
81 * E (input/output) DOUBLE PRECISION array, dimension (N)
82 * On entry, the (N-1) subdiagonal elements of the tridiagonal
83 * matrix T in elements 1 to N-1 of E. E(N) need not be set on
84 * input, but is used internally as workspace.
85 * On exit, E is overwritten.
86 *
87 * M (input) INTEGER
88 * The total number of eigenvalues found
89 * in DSTEGR2A. 0 <= M <= N.
90 *
91 * W (input) DOUBLE PRECISION array, dimension (N)
92 * The first M elements contain approximations to the selected
93 * eigenvalues in ascending order. Note that only the eigenvalues from
94 * the locally relevant part of the representation tree, that is
95 * all the clusters that include eigenvalues from DOL:DOU, are reliable
96 * on this processor. (It does not need to know about any others anyway.)
97 *
98 * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
99 * If JOBZ = 'V', and if INFO = 0, then
100 * a subset of the first M columns of Z
101 * contain the orthonormal eigenvectors of the matrix T
102 * corresponding to the selected eigenvalues, with the i-th
103 * column of Z holding the eigenvector associated with W(i).
104 * See DOL, DOU for more information.
105 *
106 * LDZ (input) INTEGER
107 * The leading dimension of the array Z. LDZ >= 1, and if
108 * JOBZ = 'V', then LDZ >= max(1,N).
109 *
110 * NZC (input) INTEGER
111 * The number of eigenvectors to be held in the array Z.
112 *
113 * ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
114 * The support of the eigenvectors in Z, i.e., the indices
115 * indicating the nonzero elements in Z. The i-th computed eigenvector
116 * is nonzero only in elements ISUPPZ( 2*i-1 ) through
117 * ISUPPZ( 2*i ). This is relevant in the case when the matrix
118 * is split. ISUPPZ is only set if N>2.
119 *
120 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
121 * On exit, if INFO = 0, WORK(1) returns the optimal
122 * (and minimal) LWORK.
123 *
124 * LWORK (input) INTEGER
125 * The dimension of the array WORK. LWORK >= max(1,18*N)
126 * if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
127 * If LWORK = -1, then a workspace query is assumed; the routine
128 * only calculates the optimal size of the WORK array, returns
129 * this value as the first entry of the WORK array, and no error
130 * message related to LWORK is issued.
131 *
132 * IWORK (workspace/output) INTEGER array, dimension (LIWORK)
133 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
134 *
135 * LIWORK (input) INTEGER
136 * The dimension of the array IWORK. LIWORK >= max(1,10*N)
137 * if the eigenvectors are desired, and LIWORK >= max(1,8*N)
138 * if only the eigenvalues are to be computed.
139 * If LIWORK = -1, then a workspace query is assumed; the
140 * routine only calculates the optimal size of the IWORK array,
141 * returns this value as the first entry of the IWORK array, and
142 * no error message related to LIWORK is issued.
143 *
144 * DOL (input) INTEGER
145 * DOU (input) INTEGER
146 * From the eigenvalues W(1:M), only eigenvectors
147 * Z(:,DOL) to Z(:,DOU) are computed.
148 * If DOL > 1, then Z(:,DOL-1-ZOFFSET) is used and overwritten.
149 * If DOU < M, then Z(:,DOU+1-ZOFFSET) is used and overwritten.
150 *
151 * NEEDIL (input/output) INTEGER
152 * NEEDIU (input/output) INTEGER
153 * Describes which are the left and right outermost eigenvalues
154 * still to be computed. Initially computed by DLARRE2A,
155 * modified in the course of the algorithm.
156 *
157 * INDWLC (output) DOUBLE PRECISION
158 * Pointer into the workspace, location where the local
159 * eigenvalue representations are stored. ("Local eigenvalues"
160 * are those relative to the individual shifts of the RRRs.)
161 *
162 * PIVMIN (input) DOUBLE PRECISION
163 * The minimum pivot in the sturm sequence for T.
164 *
165 * SCALE (input) DOUBLE PRECISION
166 * The scaling factor for T. Used for unscaling the eigenvalues
167 * at the very end of the algorithm.
168 *
169 * WL (input) DOUBLE PRECISION
170 * WU (input) DOUBLE PRECISION
171 * The interval (WL, WU] contains all the wanted eigenvalues.
172 *
173 * VSTART (input/output) LOGICAL
174 * .TRUE. on initialization, set to .FALSE. afterwards.
175 *
176 * FINISH (input/output) LOGICAL
177 * indicates whether all eigenpairs have been computed
178 *
179 * MAXCLS (input/output) INTEGER
180 * The largest cluster worked on by this processor in the
181 * representation tree.
182 *
183 * NDEPTH (input/output) INTEGER
184 * The current depth of the representation tree. Set to
185 * zero on initial pass, changed when the deeper levels of
186 * the representation tree are generated.
187 *
188 * PARITY (input/output) INTEGER
189 * An internal parameter needed for the storage of the
190 * clusters on the current level of the representation tree.
191 *
192 * ZOFFSET (input) INTEGER
193 * Offset for storing the eigenpairs when Z is distributed
194 * in 1D-cyclic fashion
195 *
196 * INFO (output) INTEGER
197 * On exit, INFO
198 * = 0: successful exit
199 * other:if INFO = -i, the i-th argument had an illegal value
200 * if INFO = 20X, internal error in DLARRV2.
201 * Here, the digit X = ABS( IINFO ) < 10, where IINFO is
202 * the nonzero error code returned by DLARRV2.
203 *
204 * .. Parameters ..
205  DOUBLE PRECISION ONE, FOUR, MINRGP
206  PARAMETER ( ONE = 1.0d0,
207  $ four = 4.0d0,
208  $ minrgp = 1.0d-3 )
209 * ..
210 * .. Local Scalars ..
211  LOGICAL LQUERY, WANTZ, ZQUERY
212  INTEGER IINDBL, IINDW, IINDWK, IINFO, IINSPL, INDERR,
213  $ INDGP, INDGRS, INDSDM, INDWRK, ITMP, J, LIWMIN,
214  $ LWMIN
215  DOUBLE PRECISION EPS, RTOL1, RTOL2
216 * ..
217 * .. External Functions ..
218  LOGICAL LSAME
219  DOUBLE PRECISION DLAMCH, DLANST
220  EXTERNAL lsame, dlamch, dlanst
221 * ..
222 * .. External Subroutines ..
223  EXTERNAL dlarrv2, dscal
224 * ..
225 * .. Intrinsic Functions ..
226  INTRINSIC dble, max, min, sqrt
227 * ..
228 * .. Executable Statements ..
229 *
230 * Test the input parameters.
231 *
232  wantz = lsame( jobz, 'V' )
233 *
234  lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
235  zquery = ( nzc.EQ.-1 )
236 
237 * DSTEGR2B needs WORK of size 6*N, IWORK of size 3*N.
238 * In addition, DLARRE2A needed WORK of size 6*N, IWORK of size 5*N.
239 * Workspace is kept consistent even though DLARRE2A is not called here.
240 * Furthermore, DLARRV2 needs WORK of size 12*N, IWORK of size 7*N.
241  IF( wantz ) THEN
242  lwmin = 18*n
243  liwmin = 10*n
244  ELSE
245 * need less workspace if only the eigenvalues are wanted
246  lwmin = 12*n
247  liwmin = 8*n
248  ENDIF
249 *
250  info = 0
251 *
252 * Get machine constants.
253 *
254  eps = dlamch( 'Precision' )
255 *
256  IF( (n.EQ.0).OR.(n.EQ.1) ) THEN
257  finish = .true.
258  RETURN
259  ENDIF
260 
261  IF(zquery.OR.lquery)
262  $ RETURN
263 *
264  indgrs = 1
265  inderr = 2*n + 1
266  indgp = 3*n + 1
267  indsdm = 4*n + 1
268  indwrk = 6*n + 1
269  indwlc = indwrk
270 *
271  iinspl = 1
272  iindbl = n + 1
273  iindw = 2*n + 1
274  iindwk = 3*n + 1
275 
276 * Set the tolerance parameters for bisection
277  rtol1 = four*sqrt(eps)
278  rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
279 
280 
281  IF( wantz ) THEN
282 *
283 * Compute the desired eigenvectors corresponding to the computed
284 * eigenvalues
285 *
286  CALL dlarrv2( n, wl, wu, d, e,
287  $ pivmin, iwork( iinspl ), m,
288  $ dol, dou, needil, neediu, minrgp, rtol1, rtol2,
289  $ w, work( inderr ), work( indgp ), iwork( iindbl ),
290  $ iwork( iindw ), work( indgrs ),
291  $ work( indsdm ), z, ldz,
292  $ isuppz, work( indwrk ), iwork( iindwk ),
293  $ vstart, finish,
294  $ maxcls, ndepth, parity, zoffset, iinfo )
295  IF( iinfo.NE.0 ) THEN
296  info = 200 + abs( iinfo )
297  RETURN
298  END IF
299 *
300  ELSE
301 * DLARRE2A computed eigenvalues of the (shifted) root representation
302 * DLARRV2 returns the eigenvalues of the unshifted matrix.
303 * However, if the eigenvectors are not desired by the user, we need
304 * to apply the corresponding shifts from DLARRE2A to obtain the
305 * eigenvalues of the original matrix.
306  DO 30 j = 1, m
307  itmp = iwork( iindbl+j-1 )
308  w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
309  30 CONTINUE
310 *
311  finish = .true.
312 *
313  END IF
314 *
315 
316  IF(finish) THEN
317 * All eigenpairs have been computed
318 
319 *
320 * If matrix was scaled, then rescale eigenvalues appropriately.
321 *
322  IF( scale.NE.one ) THEN
323  CALL dscal( m, one / scale, w, 1 )
324  END IF
325 *
326 * Correct M if needed
327 *
328  IF ( wantz ) THEN
329  IF( dol.NE.1 .OR. dou.NE.m ) THEN
330  m = dou - dol +1
331  ENDIF
332  ENDIF
333 *
334 * No sorting of eigenpairs is done here, done later in the
335 * calling subroutine
336 *
337  work( 1 ) = lwmin
338  iwork( 1 ) = liwmin
339  ENDIF
340 
341  RETURN
342 *
343 * End of DSTEGR2B
344 *
345  END
max
#define max(A, B)
Definition: pcgemr.c:180
dstegr2b
subroutine dstegr2b(JOBZ, N, D, E, M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK, LIWORK, DOL, DOU, NEEDIL, NEEDIU, INDWLC, PIVMIN, SCALE, WL, WU, VSTART, FINISH, MAXCLS, NDEPTH, PARITY, ZOFFSET, INFO)
Definition: dstegr2b.f:7
dlarrv2
subroutine dlarrv2(N, VL, VU, D, L, PIVMIN, ISPLIT, M, DOL, DOU, NEEDIL, NEEDIU, MINRGP, RTOL1, RTOL2, W, WERR, WGAP, IBLOCK, INDEXW, GERS, SDIAM, Z, LDZ, ISUPPZ, WORK, IWORK, VSTART, FINISH, MAXCLS, NDEPTH, PARITY, ZOFFSET, INFO)
Definition: dlarrv2.f:8
min
#define min(A, B)
Definition: pcgemr.c:181