ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
dstegr2.f
Go to the documentation of this file.
1  SUBROUTINE dstegr2( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
2  $ M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK,
3  $ LIWORK, DOL, DOU, ZOFFSET, INFO )
4 *
5 * -- ScaLAPACK auxiliary routine (version 2.0) --
6 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
7 * July 4, 2010
8 *
9 * .. Scalar Arguments ..
10  CHARACTER JOBZ, RANGE
11  INTEGER DOL, DOU, IL, INFO, IU,
12  $ ldz, nzc, liwork, lwork, m, n, zoffset
13  DOUBLE PRECISION VL, VU
14 
15 * ..
16 * .. Array Arguments ..
17  INTEGER ISUPPZ( * ), IWORK( * )
18  DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
19  DOUBLE PRECISION Z( LDZ, * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DSTEGR2 computes selected eigenvalues and, optionally, eigenvectors
26 * of a real symmetric tridiagonal matrix T. It is invoked in the
27 * ScaLAPACK MRRR driver PDSYEVR and the corresponding Hermitian
28 * version either when only eigenvalues are to be computed, or when only
29 * a single processor is used (the sequential-like case).
30 *
31 * DSTEGR2 has been adapted from LAPACK's DSTEGR. Please note the
32 * following crucial changes.
33 *
34 * 1. The calling sequence has two additional INTEGER parameters,
35 * DOL and DOU, that should satisfy M>=DOU>=DOL>=1.
36 * DSTEGR2 ONLY computes the eigenpairs
37 * corresponding to eigenvalues DOL through DOU in W. (That is,
38 * instead of computing the eigenpairs belonging to W(1)
39 * through W(M), only the eigenvectors belonging to eigenvalues
40 * W(DOL) through W(DOU) are computed. In this case, only the
41 * eigenvalues DOL:DOU are guaranteed to be fully accurate.
42 *
43 * 2. M is NOT the number of eigenvalues specified by RANGE, but is
44 * M = DOU - DOL + 1. This concerns the case where only eigenvalues
45 * are computed, but on more than one processor. Thus, in this case
46 * M refers to the number of eigenvalues computed on this processor.
47 *
48 * 3. The arrays W and Z might not contain all the wanted eigenpairs
49 * locally, instead this information is distributed over other
50 * processors.
51 *
52 * Arguments
53 * =========
54 *
55 * JOBZ (input) CHARACTER*1
56 * = 'N': Compute eigenvalues only;
57 * = 'V': Compute eigenvalues and eigenvectors.
58 *
59 * RANGE (input) CHARACTER*1
60 * = 'A': all eigenvalues will be found.
61 * = 'V': all eigenvalues in the half-open interval (VL,VU]
62 * will be found.
63 * = 'I': the IL-th through IU-th eigenvalues will be found.
64 *
65 * N (input) INTEGER
66 * The order of the matrix. N >= 0.
67 *
68 * D (input/output) DOUBLE PRECISION array, dimension (N)
69 * On entry, the N diagonal elements of the tridiagonal matrix
70 * T. On exit, D is overwritten.
71 *
72 * E (input/output) DOUBLE PRECISION array, dimension (N)
73 * On entry, the (N-1) subdiagonal elements of the tridiagonal
74 * matrix T in elements 1 to N-1 of E. E(N) need not be set on
75 * input, but is used internally as workspace.
76 * On exit, E is overwritten.
77 *
78 * VL (input) DOUBLE PRECISION
79 * VU (input) DOUBLE PRECISION
80 * If RANGE='V', the lower and upper bounds of the interval to
81 * be searched for eigenvalues. VL < VU.
82 * Not referenced if RANGE = 'A' or 'I'.
83 *
84 * IL (input) INTEGER
85 * IU (input) INTEGER
86 * If RANGE='I', the indices (in ascending order) of the
87 * smallest and largest eigenvalues to be returned.
88 * 1 <= IL <= IU <= N, if N > 0.
89 * Not referenced if RANGE = 'A' or 'V'.
90 *
91 * M (output) INTEGER
92 * Globally summed over all processors, M equals
93 * the total number of eigenvalues found. 0 <= M <= N.
94 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
95 * The local output equals M = DOU - DOL + 1.
96 *
97 * W (output) DOUBLE PRECISION array, dimension (N)
98 * The first M elements contain the selected eigenvalues in
99 * ascending order. Note that immediately after exiting this
100 * routine, only the eigenvalues from
101 * position DOL:DOU are to reliable on this processor
102 * because the eigenvalue computation is done in parallel.
103 * Other processors will hold reliable information on other
104 * parts of the W array. This information is communicated in
105 * the ScaLAPACK driver.
106 *
107 * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
108 * If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
109 * contain some of the orthonormal eigenvectors of the matrix T
110 * corresponding to the selected eigenvalues, with the i-th
111 * column of Z holding the eigenvector associated with W(i).
112 * If JOBZ = 'N', then Z is not referenced.
113 * Note: the user must ensure that at least max(1,M) columns are
114 * supplied in the array Z; if RANGE = 'V', the exact value of M
115 * is not known in advance and can be computed with a workspace
116 * query by setting NZC = -1, see below.
117 *
118 * LDZ (input) INTEGER
119 * The leading dimension of the array Z. LDZ >= 1, and if
120 * JOBZ = 'V', then LDZ >= max(1,N).
121 *
122 * NZC (input) INTEGER
123 * The number of eigenvectors to be held in the array Z.
124 * If RANGE = 'A', then NZC >= max(1,N).
125 * If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
126 * If RANGE = 'I', then NZC >= IU-IL+1.
127 * If NZC = -1, then a workspace query is assumed; the
128 * routine calculates the number of columns of the array Z that
129 * are needed to hold the eigenvectors.
130 * This value is returned as the first entry of the Z array, and
131 * no error message related to NZC is issued.
132 *
133 * ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
134 * The support of the eigenvectors in Z, i.e., the indices
135 * indicating the nonzero elements in Z. The i-th computed eigenvector
136 * is nonzero only in elements ISUPPZ( 2*i-1 ) through
137 * ISUPPZ( 2*i ). This is relevant in the case when the matrix
138 * is split. ISUPPZ is only set if N>2.
139 *
140 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
141 * On exit, if INFO = 0, WORK(1) returns the optimal
142 * (and minimal) LWORK.
143 *
144 * LWORK (input) INTEGER
145 * The dimension of the array WORK. LWORK >= max(1,18*N)
146 * if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
147 * If LWORK = -1, then a workspace query is assumed; the routine
148 * only calculates the optimal size of the WORK array, returns
149 * this value as the first entry of the WORK array, and no error
150 * message related to LWORK is issued.
151 *
152 * IWORK (workspace/output) INTEGER array, dimension (LIWORK)
153 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
154 *
155 * LIWORK (input) INTEGER
156 * The dimension of the array IWORK. LIWORK >= max(1,10*N)
157 * if the eigenvectors are desired, and LIWORK >= max(1,8*N)
158 * if only the eigenvalues are to be computed.
159 * If LIWORK = -1, then a workspace query is assumed; the
160 * routine only calculates the optimal size of the IWORK array,
161 * returns this value as the first entry of the IWORK array, and
162 * no error message related to LIWORK is issued.
163 *
164 * DOL (input) INTEGER
165 * DOU (input) INTEGER
166 * From the eigenvalues W(1:M), only eigenvectors
167 * Z(:,DOL) to Z(:,DOU) are computed.
168 * If DOL > 1, then Z(:,DOL-1-ZOFFSET) is used and overwritten.
169 * If DOU < M, then Z(:,DOU+1-ZOFFSET) is used and overwritten.
170 *
171 * ZOFFSET (input) INTEGER
172 * Offset for storing the eigenpairs when Z is distributed
173 * in 1D-cyclic fashion
174 *
175 * INFO (output) INTEGER
176 * On exit, INFO
177 * = 0: successful exit
178 * other:if INFO = -i, the i-th argument had an illegal value
179 * if INFO = 10X, internal error in DLARRE2,
180 * if INFO = 20X, internal error in DLARRV.
181 * Here, the digit X = ABS( IINFO ) < 10, where IINFO is
182 * the nonzero error code returned by DLARRE2 or
183 * DLARRV, respectively.
184 *
185 * =====================================================================
186 *
187 * .. Parameters ..
188  DOUBLE PRECISION ZERO, ONE, FOUR, MINRGP
189  PARAMETER ( ZERO = 0.0d0, one = 1.0d0,
190  $ four = 4.0d0,
191  $ minrgp = 1.0d-3 )
192 * ..
193 * .. Local Scalars ..
194  LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
195  INTEGER I, IIL, IINDBL, IINDW, IINDWK, IINFO, IINSPL,
196  $ iiu, inde2, inderr, indgp, indgrs, indwrk,
197  $ itmp, itmp2, j, jj, liwmin, lwmin, nsplit,
198  $ nzcmin
199  DOUBLE PRECISION BIGNUM, EPS, PIVMIN, RMAX, RMIN, RTOL1, RTOL2,
200  $ SAFMIN, SCALE, SMLNUM, THRESH, TMP, TNRM, WL,
201  $ wu
202 * ..
203 * .. External Functions ..
204  LOGICAL LSAME
205  DOUBLE PRECISION DLAMCH, DLANST
206  EXTERNAL lsame, dlamch, dlanst
207 * ..
208 * .. External Subroutines ..
209  EXTERNAL dcopy, dlae2, dlaev2, dlarrc, dlarre2,
210  $ dlarrv, dlasrt, dscal, dswap
211 * ..
212 * .. Intrinsic Functions ..
213  INTRINSIC dble, max, min, sqrt
214 * ..
215 * .. Executable Statements ..
216 *
217 * Test the input parameters.
218 *
219  wantz = lsame( jobz, 'V' )
220  alleig = lsame( range, 'A' )
221  valeig = lsame( range, 'V' )
222  indeig = lsame( range, 'I' )
223 *
224  lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
225  zquery = ( nzc.EQ.-1 )
226 
227 * DSTEGR2 needs WORK of size 6*N, IWORK of size 3*N.
228 * In addition, DLARRE2 needs WORK of size 6*N, IWORK of size 5*N.
229 * Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N.
230  IF( wantz ) THEN
231  lwmin = 18*n
232  liwmin = 10*n
233  ELSE
234 * need less workspace if only the eigenvalues are wanted
235  lwmin = 12*n
236  liwmin = 8*n
237  ENDIF
238 
239  wl = zero
240  wu = zero
241  iil = 0
242  iiu = 0
243 
244  IF( valeig ) THEN
245 * We do not reference VL, VU in the cases RANGE = 'I','A'
246 * The interval (WL, WU] contains all the wanted eigenvalues.
247 * It is either given by the user or computed in DLARRE2.
248  wl = vl
249  wu = vu
250  ELSEIF( indeig ) THEN
251 * We do not reference IL, IU in the cases RANGE = 'V','A'
252  iil = il
253  iiu = iu
254  ENDIF
255 *
256  info = 0
257  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
258  info = -1
259  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
260  info = -2
261  ELSE IF( n.LT.0 ) THEN
262  info = -3
263  ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
264  info = -7
265  ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
266  info = -8
267  ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
268  info = -9
269  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
270  info = -13
271  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
272  info = -17
273  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
274  info = -19
275  END IF
276 *
277 * Get machine constants.
278 *
279  safmin = dlamch( 'Safe minimum' )
280  eps = dlamch( 'Precision' )
281  smlnum = safmin / eps
282  bignum = one / smlnum
283  rmin = sqrt( smlnum )
284  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
285 *
286  IF( info.EQ.0 ) THEN
287  work( 1 ) = lwmin
288  iwork( 1 ) = liwmin
289 *
290  IF( wantz .AND. alleig ) THEN
291  nzcmin = n
292  iil = 1
293  iiu = n
294  ELSE IF( wantz .AND. valeig ) THEN
295  CALL dlarrc( 'T', n, vl, vu, d, e, safmin,
296  $ nzcmin, itmp, itmp2, info )
297  iil = itmp+1
298  iiu = itmp2
299  ELSE IF( wantz .AND. indeig ) THEN
300  nzcmin = iiu-iil+1
301  ELSE
302 * WANTZ .EQ. FALSE.
303  nzcmin = 0
304  ENDIF
305  IF( zquery .AND. info.EQ.0 ) THEN
306  z( 1,1 ) = nzcmin
307  ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
308  info = -14
309  END IF
310  END IF
311 
312  IF ( wantz ) THEN
313  IF ( dol.LT.1 .OR. dol.GT.nzcmin ) THEN
314  info = -20
315  ENDIF
316  IF ( dou.LT.1 .OR. dou.GT.nzcmin .OR. dou.LT.dol) THEN
317  info = -21
318  ENDIF
319  ENDIF
320 
321  IF( info.NE.0 ) THEN
322 *
323 C Disable sequential error handler
324 C for parallel case
325 C CALL XERBLA( 'DSTEGR2', -INFO )
326 *
327  RETURN
328  ELSE IF( lquery .OR. zquery ) THEN
329  RETURN
330  END IF
331 *
332 * Quick return if possible
333 *
334  m = 0
335  IF( n.EQ.0 )
336  $ RETURN
337 *
338  IF( n.EQ.1 ) THEN
339  IF( alleig .OR. indeig ) THEN
340  m = 1
341  w( 1 ) = d( 1 )
342  ELSE
343  IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
344  m = 1
345  w( 1 ) = d( 1 )
346  END IF
347  END IF
348  IF( wantz )
349  $ z( 1, 1 ) = one
350  RETURN
351  END IF
352 *
353  indgrs = 1
354  inderr = 2*n + 1
355  indgp = 3*n + 1
356  inde2 = 5*n + 1
357  indwrk = 6*n + 1
358 *
359  iinspl = 1
360  iindbl = n + 1
361  iindw = 2*n + 1
362  iindwk = 3*n + 1
363 *
364 * Scale matrix to allowable range, if necessary.
365 *
366  scale = one
367  tnrm = dlanst( 'M', n, d, e )
368  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
369  scale = rmin / tnrm
370  ELSE IF( tnrm.GT.rmax ) THEN
371  scale = rmax / tnrm
372  END IF
373  IF( scale.NE.one ) THEN
374  CALL dscal( n, scale, d, 1 )
375  CALL dscal( n-1, scale, e, 1 )
376  tnrm = tnrm*scale
377  IF( valeig ) THEN
378 * If eigenvalues in interval have to be found,
379 * scale (WL, WU] accordingly
380  wl = wl*scale
381  wu = wu*scale
382  ENDIF
383  END IF
384 *
385 * Compute the desired eigenvalues of the tridiagonal after splitting
386 * into smaller subblocks if the corresponding off-diagonal elements
387 * are small
388 * THRESH is the splitting parameter for DLARRE2
389 * A negative THRESH forces the old splitting criterion based on the
390 * size of the off-diagonal. A positive THRESH switches to splitting
391 * which preserves relative accuracy.
392 *
393  iinfo = -1
394 * Set the splitting criterion
395  IF (iinfo.EQ.0) THEN
396  thresh = eps
397  ELSE
398  thresh = -eps
399  ENDIF
400 *
401 * Store the squares of the offdiagonal values of T
402  DO 5 j = 1, n-1
403  work( inde2+j-1 ) = e(j)**2
404  5 CONTINUE
405 
406 * Set the tolerance parameters for bisection
407  IF( .NOT.wantz ) THEN
408 * DLARRE2 computes the eigenvalues to full precision.
409  rtol1 = four * eps
410  rtol2 = four * eps
411  ELSE
412 * DLARRE2 computes the eigenvalues to less than full precision.
413 * DLARRV will refine the eigenvalue approximations, and we can
414 * need less accurate initial bisection in DLARRE2.
415 * Note: these settings do only affect the subset case and DLARRE2
416  rtol1 = sqrt(eps)
417  rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
418  ENDIF
419  CALL dlarre2( range, n, wl, wu, iil, iiu, d, e,
420  $ work(inde2), rtol1, rtol2, thresh, nsplit,
421  $ iwork( iinspl ), m, dol, dou,
422  $ w, work( inderr ),
423  $ work( indgp ), iwork( iindbl ),
424  $ iwork( iindw ), work( indgrs ), pivmin,
425  $ work( indwrk ), iwork( iindwk ), iinfo )
426  IF( iinfo.NE.0 ) THEN
427  info = 100 + abs( iinfo )
428  RETURN
429  END IF
430 * Note that if RANGE .NE. 'V', DLARRE2 computes bounds on the desired
431 * part of the spectrum. All desired eigenvalues are contained in
432 * (WL,WU]
433 
434 
435  IF( wantz ) THEN
436 *
437 * Compute the desired eigenvectors corresponding to the computed
438 * eigenvalues
439 *
440  CALL dlarrv( n, wl, wu, d, e,
441  $ pivmin, iwork( iinspl ), m,
442  $ dol, dou, minrgp, rtol1, rtol2,
443  $ w, work( inderr ), work( indgp ), iwork( iindbl ),
444  $ iwork( iindw ), work( indgrs ), z, ldz,
445  $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
446  IF( iinfo.NE.0 ) THEN
447  info = 200 + abs( iinfo )
448  RETURN
449  END IF
450  ELSE
451 * DLARRE2 computes eigenvalues of the (shifted) root representation
452 * DLARRV returns the eigenvalues of the unshifted matrix.
453 * However, if the eigenvectors are not desired by the user, we need
454 * to apply the corresponding shifts from DLARRE2 to obtain the
455 * eigenvalues of the original matrix.
456  DO 20 j = 1, m
457  itmp = iwork( iindbl+j-1 )
458  w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
459  20 CONTINUE
460  END IF
461 *
462 
463 *
464 * If matrix was scaled, then rescale eigenvalues appropriately.
465 *
466  IF( scale.NE.one ) THEN
467  CALL dscal( m, one / scale, w, 1 )
468  END IF
469 *
470 * Correct M if needed
471 *
472  IF ( wantz ) THEN
473  IF( dol.NE.1 .OR. dou.NE.m ) THEN
474  m = dou - dol +1
475  ENDIF
476  ENDIF
477 *
478 * If eigenvalues are not in increasing order, then sort them,
479 * possibly along with eigenvectors.
480 *
481  IF( nsplit.GT.1 ) THEN
482  IF( .NOT. wantz ) THEN
483  CALL dlasrt( 'I', dou - dol +1, w(dol), iinfo )
484  IF( iinfo.NE.0 ) THEN
485  info = 3
486  RETURN
487  END IF
488  ELSE
489  DO 60 j = dol, dou - 1
490  i = 0
491  tmp = w( j )
492  DO 50 jj = j + 1, m
493  IF( w( jj ).LT.tmp ) THEN
494  i = jj
495  tmp = w( jj )
496  END IF
497  50 CONTINUE
498  IF( i.NE.0 ) THEN
499  w( i ) = w( j )
500  w( j ) = tmp
501  IF( wantz ) THEN
502  CALL dswap( n, z( 1, i-zoffset ),
503  $ 1, z( 1, j-zoffset ), 1 )
504  itmp = isuppz( 2*i-1 )
505  isuppz( 2*i-1 ) = isuppz( 2*j-1 )
506  isuppz( 2*j-1 ) = itmp
507  itmp = isuppz( 2*i )
508  isuppz( 2*i ) = isuppz( 2*j )
509  isuppz( 2*j ) = itmp
510  END IF
511  END IF
512  60 CONTINUE
513  END IF
514  ENDIF
515 *
516  work( 1 ) = lwmin
517  iwork( 1 ) = liwmin
518  RETURN
519 *
520 * End of DSTEGR2
521 *
522  END
max
#define max(A, B)
Definition: pcgemr.c:180
dlarre2
subroutine dlarre2(RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, DOL, DOU, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO)
Definition: dlarre2.f:6
dstegr2
subroutine dstegr2(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK, LIWORK, DOL, DOU, ZOFFSET, INFO)
Definition: dstegr2.f:4
min
#define min(A, B)
Definition: pcgemr.c:181