ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
sstegr2a.f
Go to the documentation of this file.
1  SUBROUTINE sstegr2a( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
2  $ M, W, Z, LDZ, NZC, WORK, LWORK, IWORK,
3  $ LIWORK, DOL, DOU, NEEDIL, NEEDIU,
4  $ INDERR, NSPLIT, PIVMIN, SCALE, WL, WU,
5  $ INFO )
6 *
7 * -- ScaLAPACK auxiliary routine (version 2.0) --
8 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
9 * July 4, 2010
10 *
11  IMPLICIT NONE
12 *
13 * .. Scalar Arguments ..
14  CHARACTER JOBZ, RANGE
15  INTEGER DOL, DOU, IL, INDERR, INFO, IU, LDZ, LIWORK,
16  $ LWORK, M, N, NEEDIL, NEEDIU, NSPLIT, NZC
17  REAL PIVMIN, SCALE, VL, VU, WL, WU
18 
19 * ..
20 * .. Array Arguments ..
21  INTEGER IWORK( * )
22  REAL D( * ), E( * ), W( * ), WORK( * )
23  REAL Z( LDZ, * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * SSTEGR2A computes selected eigenvalues and initial representations.
30 * needed for eigenvector computations in SSTEGR2B. It is invoked in the
31 * ScaLAPACK MRRR driver PSSYEVR and the corresponding Hermitian
32 * version when both eigenvalues and eigenvectors are computed in parallel.
33 * on multiple processors. For this case, SSTEGR2A implements the FIRST
34 * part of the MRRR algorithm, parallel eigenvalue computation and finding
35 * the root RRR. At the end of SSTEGR2A,
36 * other processors might have a part of the spectrum that is needed to
37 * continue the computation locally. Once this eigenvalue information has
38 * been received by the processor, the computation can then proceed by calling
39 * the SECOND part of the parallel MRRR algorithm, SSTEGR2B.
40 *
41 * Please note:
42 * 1. The calling sequence has two additional INTEGER parameters,
43 * (compared to LAPACK's SSTEGR), these are
44 * DOL and DOU and should satisfy M>=DOU>=DOL>=1.
45 * These parameters are only relevant for the case JOBZ = 'V'.
46 *
47 * Globally invoked over all processors, SSTEGR2A computes
48 * ALL the eigenVALUES specified by RANGE.
49 * RANGE= 'A': all eigenvalues will be found.
50 * = 'V': all eigenvalues in (VL,VU] will be found.
51 * = 'I': the IL-th through IU-th eigenvalues will be found.
52 *
53 * SSTEGR2A LOCALLY only computes the eigenvalues
54 * corresponding to eigenvalues DOL through DOU in W. (That is,
55 * instead of computing the eigenvectors belonging to W(1)
56 * through W(M), only the eigenvectors belonging to eigenvalues
57 * W(DOL) through W(DOU) are computed. In this case, only the
58 * eigenvalues DOL:DOU are guaranteed to be fully accurate.
59 *
60 * 2. M is NOT the number of eigenvalues specified by RANGE, but it is
61 * M = DOU - DOL + 1. Instead, M refers to the number of eigenvalues computed on
62 * this processor.
63 *
64 * 3. While no eigenvectors are computed in SSTEGR2A itself (this is
65 * done later in SSTEGR2B), the interface
66 * If JOBZ = 'V' then, depending on RANGE and DOL, DOU, SSTEGR2A
67 * might need more workspace in Z then the original SSTEGR.
68 * In particular, the arrays W and Z might not contain all the wanted eigenpairs
69 * locally, instead this information is distributed over other
70 * processors.
71 *
72 * Arguments
73 * =========
74 *
75 * JOBZ (input) CHARACTER*1
76 * = 'N': Compute eigenvalues only;
77 * = 'V': Compute eigenvalues and eigenvectors.
78 *
79 * RANGE (input) CHARACTER*1
80 * = 'A': all eigenvalues will be found.
81 * = 'V': all eigenvalues in the half-open interval (VL,VU]
82 * will be found.
83 * = 'I': the IL-th through IU-th eigenvalues will be found.
84 *
85 * N (input) INTEGER
86 * The order of the matrix. N >= 0.
87 *
88 * D (input/output) REAL array, dimension (N)
89 * On entry, the N diagonal elements of the tridiagonal matrix
90 * T. On exit, D is overwritten.
91 *
92 * E (input/output) REAL array, dimension (N)
93 * On entry, the (N-1) subdiagonal elements of the tridiagonal
94 * matrix T in elements 1 to N-1 of E. E(N) need not be set on
95 * input, but is used internally as workspace.
96 * On exit, E is overwritten.
97 *
98 * VL (input) REAL
99 * VU (input) REAL
100 * If RANGE='V', the lower and upper bounds of the interval to
101 * be searched for eigenvalues. VL < VU.
102 * Not referenced if RANGE = 'A' or 'I'.
103 *
104 * IL (input) INTEGER
105 * IU (input) INTEGER
106 * If RANGE='I', the indices (in ascending order) of the
107 * smallest and largest eigenvalues to be returned.
108 * 1 <= IL <= IU <= N, if N > 0.
109 * Not referenced if RANGE = 'A' or 'V'.
110 *
111 * M (output) INTEGER
112 * Globally summed over all processors, M equals
113 * the total number of eigenvalues found. 0 <= M <= N.
114 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
115 * The local output equals M = DOU - DOL + 1.
116 *
117 * W (output) REAL array, dimension (N)
118 * The first M elements contain approximations to the selected
119 * eigenvalues in ascending order. Note that immediately after
120 * exiting this routine, only the eigenvalues from
121 * position DOL:DOU are to reliable on this processor
122 * because the eigenvalue computation is done in parallel.
123 * The other entries outside DOL:DOU are very crude preliminary
124 * approximations. Other processors hold reliable information on
125 * these other parts of the W array.
126 * This information is communicated in the ScaLAPACK driver.
127 *
128 * Z (output) REAL array, dimension (LDZ, max(1,M) )
129 * SSTEGR2A does not compute eigenvectors, this is done
130 * in SSTEGR2B. The argument Z as well as all related
131 * other arguments only appear to keep the interface consistent
132 * and to signal to the user that this subroutine is meant to
133 * be used when eigenvectors are computed.
134 *
135 * LDZ (input) INTEGER
136 * The leading dimension of the array Z. LDZ >= 1, and if
137 * JOBZ = 'V', then LDZ >= max(1,N).
138 *
139 * NZC (input) INTEGER
140 * The number of eigenvectors to be held in the array Z.
141 * If RANGE = 'A', then NZC >= max(1,N).
142 * If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
143 * If RANGE = 'I', then NZC >= IU-IL+1.
144 * If NZC = -1, then a workspace query is assumed; the
145 * routine calculates the number of columns of the array Z that
146 * are needed to hold the eigenvectors.
147 * This value is returned as the first entry of the Z array, and
148 * no error message related to NZC is issued.
149 *
150 * WORK (workspace/output) REAL array, dimension (LWORK)
151 * On exit, if INFO = 0, WORK(1) returns the optimal
152 * (and minimal) LWORK.
153 *
154 * LWORK (input) INTEGER
155 * The dimension of the array WORK. LWORK >= max(1,18*N)
156 * if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
157 * If LWORK = -1, then a workspace query is assumed; the routine
158 * only calculates the optimal size of the WORK array, returns
159 * this value as the first entry of the WORK array, and no error
160 * message related to LWORK is issued.
161 *
162 * IWORK (workspace/output) INTEGER array, dimension (LIWORK)
163 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
164 *
165 * LIWORK (input) INTEGER
166 * The dimension of the array IWORK. LIWORK >= max(1,10*N)
167 * if the eigenvectors are desired, and LIWORK >= max(1,8*N)
168 * if only the eigenvalues are to be computed.
169 * If LIWORK = -1, then a workspace query is assumed; the
170 * routine only calculates the optimal size of the IWORK array,
171 * returns this value as the first entry of the IWORK array, and
172 * no error message related to LIWORK is issued.
173 *
174 * DOL (input) INTEGER
175 * DOU (input) INTEGER
176 * From all the eigenvalues W(1:M), only eigenvalues
177 * W(DOL:DOU) are computed.
178 *
179 * NEEDIL (output) INTEGER
180 * NEEDIU (output) INTEGER
181 * The indices of the leftmost and rightmost eigenvalues
182 * needed to accurately compute the relevant part of the
183 * representation tree. This information can be used to
184 * find out which processors have the relevant eigenvalue
185 * information needed so that it can be communicated.
186 *
187 * INDERR (output) INTEGER
188 * INDERR points to the place in the work space where
189 * the eigenvalue uncertainties (errors) are stored.
190 *
191 * NSPLIT (output) INTEGER
192 * The number of blocks T splits into. 1 <= NSPLIT <= N.
193 *
194 * PIVMIN (output) REAL
195 * The minimum pivot in the sturm sequence for T.
196 *
197 * SCALE (output) REAL
198 * The scaling factor for the tridiagonal T.
199 *
200 * WL (output) REAL
201 * WU (output) REAL
202 * The interval (WL, WU] contains all the wanted eigenvalues.
203 * It is either given by the user or computed in SLARRE2A.
204 *
205 * INFO (output) INTEGER
206 * On exit, INFO
207 * = 0: successful exit
208 * other:if INFO = -i, the i-th argument had an illegal value
209 * if INFO = 10X, internal error in SLARRE2A,
210 * Here, the digit X = ABS( IINFO ) < 10, where IINFO is
211 * the nonzero error code returned by SLARRE2A.
212 *
213 * =====================================================================
214 *
215 * .. Parameters ..
216  REAL ZERO, ONE, FOUR, MINRGP
217  PARAMETER ( ZERO = 0.0e0, one = 1.0e0,
218  $ four = 4.0e0,
219  $ minrgp = 3.0e-3 )
220 * ..
221 * .. Local Scalars ..
222  LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
223  INTEGER IIL, IINDBL, IINDW, IINDWK, IINFO, IINSPL, IIU,
224  $ INDE2, INDGP, INDGRS, INDSDM, INDWRK, ITMP,
225  $ ITMP2, J, LIWMIN, LWMIN, NZCMIN
226  REAL BIGNUM, EPS, RMAX, RMIN, RTOL1, RTOL2, SAFMIN,
227  $ smlnum, thresh, tnrm
228 * ..
229 * .. External Functions ..
230  LOGICAL LSAME
231  REAL SLAMCH, SLANST
232  EXTERNAL LSAME, SLAMCH, SLANST
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL slarrc, slarre2a, sscal
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC max, min, real, sqrt
239 * ..
240 * .. Executable Statements ..
241 *
242 * Test the input parameters.
243 *
244  wantz = lsame( jobz, 'V' )
245  alleig = lsame( range, 'A' )
246  valeig = lsame( range, 'V' )
247  indeig = lsame( range, 'I' )
248 *
249  lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
250  zquery = ( nzc.EQ.-1 )
251 
252 * SSTEGR2A needs WORK of size 6*N, IWORK of size 3*N.
253 * In addition, SLARRE2A needs WORK of size 6*N, IWORK of size 5*N.
254 * Furthermore, SLARRV2 needs WORK of size 12*N, IWORK of size 7*N.
255 * Workspace is kept consistent with SSTEGR2B even though
256 * SLARRV2 is not called here.
257  IF( wantz ) THEN
258  lwmin = 18*n
259  liwmin = 10*n
260  ELSE
261 * need less workspace if only the eigenvalues are wanted
262  lwmin = 12*n
263  liwmin = 8*n
264  ENDIF
265 
266  wl = zero
267  wu = zero
268  iil = 0
269  iiu = 0
270 
271  IF( valeig ) THEN
272 * We do not reference VL, VU in the cases RANGE = 'I','A'
273 * The interval (WL, WU] contains all the wanted eigenvalues.
274 * It is either given by the user or computed in SLARRE2A.
275  wl = vl
276  wu = vu
277  ELSEIF( indeig ) THEN
278 * We do not reference IL, IU in the cases RANGE = 'V','A'
279  iil = il
280  iiu = iu
281  ENDIF
282 *
283  info = 0
284  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
285  info = -1
286  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
287  info = -2
288  ELSE IF( n.LT.0 ) THEN
289  info = -3
290  ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
291  info = -7
292  ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
293  info = -8
294  ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
295  info = -9
296  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
297  info = -13
298  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
299  info = -17
300  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
301  info = -19
302  END IF
303 *
304 * Get machine constants.
305 *
306  safmin = slamch( 'Safe minimum' )
307  eps = slamch( 'Precision' )
308  smlnum = safmin / eps
309  bignum = one / smlnum
310  rmin = sqrt( smlnum )
311  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
312 *
313  IF( info.EQ.0 ) THEN
314  work( 1 ) = lwmin
315  iwork( 1 ) = liwmin
316 *
317  IF( wantz .AND. alleig ) THEN
318  nzcmin = n
319  iil = 1
320  iiu = n
321  ELSE IF( wantz .AND. valeig ) THEN
322  CALL slarrc( 'T', n, vl, vu, d, e, safmin,
323  $ nzcmin, itmp, itmp2, info )
324  iil = itmp+1
325  iiu = itmp2
326  ELSE IF( wantz .AND. indeig ) THEN
327  nzcmin = iiu-iil+1
328  ELSE
329 * WANTZ .EQ. FALSE.
330  nzcmin = 0
331  ENDIF
332  IF( zquery .AND. info.EQ.0 ) THEN
333  z( 1,1 ) = nzcmin
334  ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
335  info = -14
336  END IF
337  END IF
338 
339  IF ( wantz ) THEN
340  IF ( dol.LT.1 .OR. dol.GT.nzcmin ) THEN
341  info = -20
342  ENDIF
343  IF ( dou.LT.1 .OR. dou.GT.nzcmin .OR. dou.LT.dol) THEN
344  info = -21
345  ENDIF
346  ENDIF
347 
348  IF( info.NE.0 ) THEN
349 *
350 C Disable sequential error handler
351 C for parallel case
352 C CALL XERBLA( 'SSTEGR2A', -INFO )
353 *
354  RETURN
355  ELSE IF( lquery .OR. zquery ) THEN
356  RETURN
357  END IF
358 
359 * Initialize NEEDIL and NEEDIU, these values are changed in SLARRE2A
360  needil = dou
361  neediu = dol
362 *
363 * Quick return if possible
364 *
365  m = 0
366  IF( n.EQ.0 )
367  $ RETURN
368 *
369  IF( n.EQ.1 ) THEN
370  IF( alleig .OR. indeig ) THEN
371  m = 1
372  w( 1 ) = d( 1 )
373  ELSE
374  IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
375  m = 1
376  w( 1 ) = d( 1 )
377  END IF
378  END IF
379  IF( wantz )
380  $ z( 1, 1 ) = one
381  RETURN
382  END IF
383 *
384  indgrs = 1
385  inderr = 2*n + 1
386  indgp = 3*n + 1
387  indsdm = 4*n + 1
388  inde2 = 5*n + 1
389  indwrk = 6*n + 1
390 *
391  iinspl = 1
392  iindbl = n + 1
393  iindw = 2*n + 1
394  iindwk = 3*n + 1
395 *
396 * Scale matrix to allowable range, if necessary.
397 *
398  scale = one
399  tnrm = slanst( 'M', n, d, e )
400  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
401  scale = rmin / tnrm
402  ELSE IF( tnrm.GT.rmax ) THEN
403  scale = rmax / tnrm
404  END IF
405  IF( scale.NE.one ) THEN
406  CALL sscal( n, scale, d, 1 )
407  CALL sscal( n-1, scale, e, 1 )
408  tnrm = tnrm*scale
409  IF( valeig ) THEN
410 * If eigenvalues in interval have to be found,
411 * scale (WL, WU] accordingly
412  wl = wl*scale
413  wu = wu*scale
414  ENDIF
415  END IF
416 *
417 * Compute the desired eigenvalues of the tridiagonal after splitting
418 * into smaller subblocks if the corresponding off-diagonal elements
419 * are small
420 * THRESH is the splitting parameter for SLARRA in SLARRE2A
421 * A negative THRESH forces the old splitting criterion based on the
422 * size of the off-diagonal.
423  thresh = -eps
424  iinfo = 0
425 
426 * Store the squares of the offdiagonal values of T
427  DO 5 j = 1, n-1
428  work( inde2+j-1 ) = e(j)**2
429  5 CONTINUE
430 
431 * Set the tolerance parameters for bisection
432  IF( .NOT.wantz ) THEN
433 * SLARRE2A computes the eigenvalues to full precision.
434  rtol1 = four * eps
435  rtol2 = four * eps
436  ELSE
437 * SLARRE2A computes the eigenvalues to less than full precision.
438 * SLARRV2 will refine the eigenvalue approximations, and we can
439 * need less accurate initial bisection in SLARRE2A.
440  rtol1 = four*sqrt(eps)
441  rtol2 = max( sqrt(eps)*5.0e-3, four * eps )
442  ENDIF
443  CALL slarre2a( range, n, wl, wu, iil, iiu, d, e,
444  $ work(inde2), rtol1, rtol2, thresh, nsplit,
445  $ iwork( iinspl ), m, dol, dou, needil, neediu,
446  $ w, work( inderr ),
447  $ work( indgp ), iwork( iindbl ),
448  $ iwork( iindw ), work( indgrs ),
449  $ work( indsdm ), pivmin,
450  $ work( indwrk ), iwork( iindwk ),
451  $ minrgp, iinfo )
452  IF( iinfo.NE.0 ) THEN
453  info = 100 + abs( iinfo )
454  RETURN
455  END IF
456 * Note that if RANGE .NE. 'V', SLARRE2A computes bounds on the desired
457 * part of the spectrum. All desired eigenvalues are contained in
458 * (WL,WU]
459 
460 
461  RETURN
462 *
463 * End of SSTEGR2A
464 *
465  END
max
#define max(A, B)
Definition: pcgemr.c:180
sstegr2a
subroutine sstegr2a(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, WORK, LWORK, IWORK, LIWORK, DOL, DOU, NEEDIL, NEEDIU, INDERR, NSPLIT, PIVMIN, SCALE, WL, WU, INFO)
Definition: sstegr2a.f:6
slarre2a
subroutine slarre2a(RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, DOL, DOU, NEEDIL, NEEDIU, W, WERR, WGAP, IBLOCK, INDEXW, GERS, SDIAM, PIVMIN, WORK, IWORK, MINRGP, INFO)
Definition: slarre2a.f:6
min
#define min(A, B)
Definition: pcgemr.c:181