ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdstebz.f
Go to the documentation of this file.
1  SUBROUTINE pdstebz( ICTXT, RANGE, ORDER, N, VL, VU, IL, IU,
2  $ ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT,
3  $ WORK, LWORK, IWORK, LIWORK, INFO )
4 *
5 * -- ScaLAPACK routine (version 1.7) --
6 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7 * and University of California, Berkeley.
8 * November 15, 1997
9 *
10 * .. Scalar Arguments ..
11  CHARACTER ORDER, RANGE
12  INTEGER ICTXT, IL, INFO, IU, LIWORK, LWORK, M, N,
13  $ nsplit
14  DOUBLE PRECISION ABSTOL, VL, VU
15 * ..
16 * .. Array Arguments ..
17  INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
18  DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * PDSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix in
25 * parallel. The user may ask for all eigenvalues, all eigenvalues in
26 * the interval [VL, VU], or the eigenvalues indexed IL through IU. A
27 * static partitioning of work is done at the beginning of PDSTEBZ which
28 * results in all processes finding an (almost) equal number of
29 * eigenvalues.
30 *
31 * NOTE : It is assumed that the user is on an IEEE machine. If the user
32 * is not on an IEEE mchine, set the compile time flag NO_IEEE
33 * to 1 (in SLmake.inc). The features of IEEE arithmetic that
34 * are needed for the "fast" Sturm Count are : (a) infinity
35 * arithmetic (b) the sign bit of a single precision floating
36 * point number is assumed be in the 32nd bit position
37 * (c) the sign of negative zero.
38 *
39 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
40 * Matrix", Report CS41, Computer Science Dept., Stanford
41 * University, July 21, 1966.
42 *
43 * Arguments
44 * =========
45 *
46 * ICTXT (global input) INTEGER
47 * The BLACS context handle.
48 *
49 * RANGE (global input) CHARACTER
50 * Specifies which eigenvalues are to be found.
51 * = 'A': ("All") all eigenvalues will be found.
52 * = 'V': ("Value") all eigenvalues in the interval
53 * [VL, VU] will be found.
54 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
55 * entire matrix) will be found.
56 *
57 * ORDER (global input) CHARACTER
58 * Specifies the order in which the eigenvalues and their block
59 * numbers are stored in W and IBLOCK.
60 * = 'B': ("By Block") the eigenvalues will be grouped by
61 * split-off block (see IBLOCK, ISPLIT) and
62 * ordered from smallest to largest within
63 * the block.
64 * = 'E': ("Entire matrix")
65 * the eigenvalues for the entire matrix
66 * will be ordered from smallest to largest.
67 *
68 * N (global input) INTEGER
69 * The order of the tridiagonal matrix T. N >= 0.
70 *
71 * VL (global input) DOUBLE PRECISION
72 * If RANGE='V', the lower bound of the interval to be searched
73 * for eigenvalues. Eigenvalues less than VL will not be
74 * returned. Not referenced if RANGE='A' or 'I'.
75 *
76 * VU (global input) DOUBLE PRECISION
77 * If RANGE='V', the upper bound of the interval to be searched
78 * for eigenvalues. Eigenvalues greater than VU will not be
79 * returned. VU must be greater than VL. Not referenced if
80 * RANGE='A' or 'I'.
81 *
82 * IL (global input) INTEGER
83 * If RANGE='I', the index (from smallest to largest) of the
84 * smallest eigenvalue to be returned. IL must be at least 1.
85 * Not referenced if RANGE='A' or 'V'.
86 *
87 * IU (global input) INTEGER
88 * If RANGE='I', the index (from smallest to largest) of the
89 * largest eigenvalue to be returned. IU must be at least IL
90 * and no greater than N. Not referenced if RANGE='A' or 'V'.
91 *
92 * ABSTOL (global input) DOUBLE PRECISION
93 * The absolute tolerance for the eigenvalues. An eigenvalue
94 * (or cluster) is considered to be located if it has been
95 * determined to lie in an interval whose width is ABSTOL or
96 * less. If ABSTOL is less than or equal to zero, then ULP*|T|
97 * will be used, where |T| means the 1-norm of T.
98 * Eigenvalues will be computed most accurately when ABSTOL is
99 * set to the underflow threshold DLAMCH('U'), not zero.
100 * Note : If eigenvectors are desired later by inverse iteration
101 * ( PDSTEIN ), ABSTOL should be set to 2*PDLAMCH('S').
102 *
103 * D (global input) DOUBLE PRECISION array, dimension (N)
104 * The n diagonal elements of the tridiagonal matrix T. To
105 * avoid overflow, the matrix must be scaled so that its largest
106 * entry is no greater than overflow**(1/2) * underflow**(1/4)
107 * in absolute value, and for greatest accuracy, it should not
108 * be much smaller than that.
109 *
110 * E (global input) DOUBLE PRECISION array, dimension (N-1)
111 * The (n-1) off-diagonal elements of the tridiagonal matrix T.
112 * To avoid overflow, the matrix must be scaled so that its
113 * largest entry is no greater than overflow**(1/2) *
114 * underflow**(1/4) in absolute value, and for greatest
115 * accuracy, it should not be much smaller than that.
116 *
117 * M (global output) INTEGER
118 * The actual number of eigenvalues found. 0 <= M <= N.
119 * (See also the description of INFO=2)
120 *
121 * NSPLIT (global output) INTEGER
122 * The number of diagonal blocks in the matrix T.
123 * 1 <= NSPLIT <= N.
124 *
125 * W (global output) DOUBLE PRECISION array, dimension (N)
126 * On exit, the first M elements of W contain the eigenvalues
127 * on all processes.
128 *
129 * IBLOCK (global output) INTEGER array, dimension (N)
130 * At each row/column j where E(j) is zero or small, the
131 * matrix T is considered to split into a block diagonal
132 * matrix. On exit IBLOCK(i) specifies which block (from 1
133 * to the number of blocks) the eigenvalue W(i) belongs to.
134 * NOTE: in the (theoretically impossible) event that bisection
135 * does not converge for some or all eigenvalues, INFO is set
136 * to 1 and the ones for which it did not are identified by a
137 * negative block number.
138 *
139 * ISPLIT (global output) INTEGER array, dimension (N)
140 * The splitting points, at which T breaks up into submatrices.
141 * The first submatrix consists of rows/columns 1 to ISPLIT(1),
142 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
143 * etc., and the NSPLIT-th consists of rows/columns
144 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
145 * (Only the first NSPLIT elements will actually be used, but
146 * since the user cannot know a priori what value NSPLIT will
147 * have, N words must be reserved for ISPLIT.)
148 *
149 * WORK (local workspace) DOUBLE PRECISION array,
150 * dimension ( MAX( 5*N, 7 ) )
151 *
152 * LWORK (local input) INTEGER
153 * size of array WORK must be >= MAX( 5*N, 7 )
154 * If LWORK = -1, then LWORK is global input and a workspace
155 * query is assumed; the routine only calculates the minimum
156 * and optimal size for all work arrays. Each of these
157 * values is returned in the first entry of the corresponding
158 * work array, and no error message is issued by PXERBLA.
159 *
160 * IWORK (local workspace) INTEGER array, dimension ( MAX( 4*N, 14 ) )
161 *
162 * LIWORK (local input) INTEGER
163 * size of array IWORK must be >= MAX( 4*N, 14, NPROCS )
164 * If LIWORK = -1, then LIWORK is global input and a workspace
165 * query is assumed; the routine only calculates the minimum
166 * and optimal size for all work arrays. Each of these
167 * values is returned in the first entry of the corresponding
168 * work array, and no error message is issued by PXERBLA.
169 *
170 * INFO (global output) INTEGER
171 * = 0 : successful exit
172 * < 0 : if INFO = -i, the i-th argument had an illegal value
173 * > 0 : some or all of the eigenvalues failed to converge or
174 * were not computed:
175 * = 1 : Bisection failed to converge for some eigenvalues;
176 * these eigenvalues are flagged by a negative block
177 * number. The effect is that the eigenvalues may not
178 * be as accurate as the absolute and relative
179 * tolerances. This is generally caused by arithmetic
180 * which is less accurate than PDLAMCH says.
181 * = 2 : There is a mismatch between the number of
182 * eigenvalues output and the number desired.
183 * = 3 : RANGE='i', and the Gershgorin interval initially
184 * used was incorrect. No eigenvalues were computed.
185 * Probable cause: your machine has sloppy floating
186 * point arithmetic.
187 * Cure: Increase the PARAMETER "FUDGE", recompile,
188 * and try again.
189 *
190 * Internal Parameters
191 * ===================
192 *
193 * RELFAC DOUBLE PRECISION, default = 2.0
194 * The relative tolerance. An interval [a,b] lies within
195 * "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
196 * where "ulp" is the machine precision (distance from 1 to
197 * the next larger floating point number.)
198 *
199 * FUDGE DOUBLE PRECISION, default = 2.0
200 * A "fudge factor" to widen the Gershgorin intervals. Ideally,
201 * a value of 1 should work, but on machines with sloppy
202 * arithmetic, this needs to be larger. The default for
203 * publicly released versions should be large enough to handle
204 * the worst machine around. Note that this has no effect
205 * on the accuracy of the solution.
206 *
207 * =====================================================================
208 *
209 * .. Intrinsic Functions ..
210  INTRINSIC abs, dble, ichar, max, min, mod
211 * ..
212 * .. External Functions ..
213  LOGICAL LSAME
214  INTEGER BLACS_PNUM
215  DOUBLE PRECISION PDLAMCH
216  EXTERNAL lsame, blacs_pnum, pdlamch
217 * ..
218 * .. External Subroutines ..
219  EXTERNAL blacs_freebuff, blacs_get, blacs_gridexit,
220  $ blacs_gridinfo, blacs_gridmap, dgebr2d,
221  $ dgebs2d, dgerv2d, dgesd2d, dlasrt2, globchk,
222  $ igebr2d, igebs2d, igerv2d, igesd2d, igsum2d,
223  $ pdlaebz, pdlaiectb, pdlaiectl, pdlapdct,
224  $ pdlasnbt, pxerbla
225 * ..
226 * .. Parameters ..
227  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
228  $ MB_, NB_, RSRC_, CSRC_, LLD_
229  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
230  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
231  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
232  INTEGER BIGNUM, DESCMULT
233  PARAMETER ( BIGNUM = 10000, descmult = 100 )
234  DOUBLE PRECISION ZERO, ONE, TWO, FIVE, HALF
235  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
236  $ five = 5.0d+0, half = 1.0d+0 / two )
237  DOUBLE PRECISION FUDGE, RELFAC
238  PARAMETER ( FUDGE = 2.0d+0, relfac = 2.0d+0 )
239 * ..
240 * .. Local Scalars ..
241  LOGICAL LQUERY
242  INTEGER BLKNO, FOUND, I, IBEGIN, IEFLAG, IEND, IFRST,
243  $ iinfo, ilast, iload, im, imyload, in, indriw1,
244  $ indriw2, indrw1, indrw2, inxtload, ioff,
245  $ iorder, iout, irange, irecv, irem, itmp1,
246  $ itmp2, j, jb, k, last, lextra, lreq, mycol,
247  $ myrow, nalpha, nbeta, ncmp, neigint, next, ngl,
248  $ nglob, ngu, nint, npcol, nprow, offset,
249  $ onedcontext, p, prev, rextra, rreq, self,
250  $ torecv
251  DOUBLE PRECISION ALPHA, ATOLI, BETA, BNORM, DRECV, DSEND, GL,
252  $ GU, INITVL, INITVU, LSAVE, MID, PIVMIN, RELTOL,
253  $ safemn, tmp1, tmp2, tnorm, ulp
254 * ..
255 * .. Local Arrays ..
256  INTEGER IDUM( 5, 2 )
257 * ..
258 * .. Executable Statements ..
259 * This is just to keep ftnchek happy
260  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
261  $ rsrc_.LT.0 )RETURN
262 *
263 * Set up process grid
264 *
265  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
266 *
267  info = 0
268  m = 0
269 *
270 * Decode RANGE
271 *
272  IF( lsame( range, 'A' ) ) THEN
273  irange = 1
274  ELSE IF( lsame( range, 'V' ) ) THEN
275  irange = 2
276  ELSE IF( lsame( range, 'I' ) ) THEN
277  irange = 3
278  ELSE
279  irange = 0
280  END IF
281 *
282 * Decode ORDER
283 *
284  IF( lsame( order, 'B' ) ) THEN
285  iorder = 2
286  ELSE IF( lsame( order, 'E' ) .OR. lsame( order, 'A' ) ) THEN
287  iorder = 1
288  ELSE
289  iorder = 0
290  END IF
291 *
292 * Check for Errors
293 *
294  IF( nprow.EQ.-1 ) THEN
295  info = -1
296  ELSE
297 *
298 * Get machine constants
299 *
300  safemn = pdlamch( ictxt, 'S' )
301  ulp = pdlamch( ictxt, 'P' )
302  reltol = ulp*relfac
303  idum( 1, 1 ) = ichar( range )
304  idum( 1, 2 ) = 2
305  idum( 2, 1 ) = ichar( order )
306  idum( 2, 2 ) = 3
307  idum( 3, 1 ) = n
308  idum( 3, 2 ) = 4
309  nglob = 5
310  IF( irange.EQ.3 ) THEN
311  idum( 4, 1 ) = il
312  idum( 4, 2 ) = 7
313  idum( 5, 1 ) = iu
314  idum( 5, 2 ) = 8
315  ELSE
316  idum( 4, 1 ) = 0
317  idum( 4, 2 ) = 0
318  idum( 5, 1 ) = 0
319  idum( 5, 2 ) = 0
320  END IF
321  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
322  work( 1 ) = abstol
323  IF( irange.EQ.2 ) THEN
324  work( 2 ) = vl
325  work( 3 ) = vu
326  ELSE
327  work( 2 ) = zero
328  work( 3 ) = zero
329  END IF
330  CALL dgebs2d( ictxt, 'ALL', ' ', 3, 1, work, 3 )
331  ELSE
332  CALL dgebr2d( ictxt, 'ALL', ' ', 3, 1, work, 3, 0, 0 )
333  END IF
334  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
335  IF( info.EQ.0 ) THEN
336  IF( irange.EQ.0 ) THEN
337  info = -2
338  ELSE IF( iorder.EQ.0 ) THEN
339  info = -3
340  ELSE IF( irange.EQ.2 .AND. vl.GE.vu ) THEN
341  info = -5
342  ELSE IF( irange.EQ.3 .AND. ( il.LT.1 .OR. il.GT.max( 1,
343  $ n ) ) ) THEN
344  info = -6
345  ELSE IF( irange.EQ.3 .AND. ( iu.LT.min( n,
346  $ il ) .OR. iu.GT.n ) ) THEN
347  info = -7
348  ELSE IF( lwork.LT.max( 5*n, 7 ) .AND. .NOT.lquery ) THEN
349  info = -18
350  ELSE IF( liwork.LT.max( 4*n, 14, nprow*npcol ) .AND. .NOT.
351  $ lquery ) THEN
352  info = -20
353  ELSE IF( irange.EQ.2 .AND. ( abs( work( 2 )-vl ).GT.five*
354  $ ulp*abs( vl ) ) ) THEN
355  info = -5
356  ELSE IF( irange.EQ.2 .AND. ( abs( work( 3 )-vu ).GT.five*
357  $ ulp*abs( vu ) ) ) THEN
358  info = -6
359  ELSE IF( abs( work( 1 )-abstol ).GT.five*ulp*abs( abstol ) )
360  $ THEN
361  info = -9
362  END IF
363  END IF
364  IF( info.EQ.0 )
365  $ info = bignum
366  CALL globchk( ictxt, nglob, idum, 5, iwork, info )
367  IF( info.EQ.bignum ) THEN
368  info = 0
369  ELSE IF( mod( info, descmult ).EQ.0 ) THEN
370  info = -info / descmult
371  ELSE
372  info = -info
373  END IF
374  END IF
375  work( 1 ) = dble( max( 5*n, 7 ) )
376  iwork( 1 ) = max( 4*n, 14, nprow*npcol )
377 *
378  IF( info.NE.0 ) THEN
379  CALL pxerbla( ictxt, 'PDSTEBZ', -info )
380  RETURN
381  ELSE IF( lwork.EQ.-1 .AND. liwork.EQ.-1 ) THEN
382  RETURN
383  END IF
384 *
385 *
386 * Quick return if possible
387 *
388  IF( n.EQ.0 )
389  $ RETURN
390 *
391  k = 1
392  DO 20 i = 0, nprow - 1
393  DO 10 j = 0, npcol - 1
394  iwork( k ) = blacs_pnum( ictxt, i, j )
395  k = k + 1
396  10 CONTINUE
397  20 CONTINUE
398 *
399  p = nprow*npcol
400  nprow = 1
401  npcol = p
402 *
403  CALL blacs_get( ictxt, 10, onedcontext )
404  CALL blacs_gridmap( onedcontext, iwork, nprow, nprow, npcol )
405  CALL blacs_gridinfo( onedcontext, i, j, k, self )
406 *
407 * Simplifications:
408 *
409  IF( irange.EQ.3 .AND. il.EQ.1 .AND. iu.EQ.n )
410  $ irange = 1
411 *
412  next = mod( self+1, p )
413  prev = mod( p+self-1, p )
414 *
415 * Compute squares of off-diagonals, splitting points and pivmin.
416 * Interleave diagonals and off-diagonals.
417 *
418  indrw1 = max( 2*n, 4 )
419  indrw2 = indrw1 + 2*n
420  indriw1 = max( 2*n, 8 )
421  nsplit = 1
422  work( indrw1+2*n ) = zero
423  pivmin = one
424 *
425  DO 30 i = 1, n - 1
426  tmp1 = e( i )**2
427  j = 2*i
428  work( indrw1+j-1 ) = d( i )
429  IF( abs( d( i+1 )*d( i ) )*ulp**2+safemn.GT.tmp1 ) THEN
430  isplit( nsplit ) = i
431  nsplit = nsplit + 1
432  work( indrw1+j ) = zero
433  ELSE
434  work( indrw1+j ) = tmp1
435  pivmin = max( pivmin, tmp1 )
436  END IF
437  30 CONTINUE
438  work( indrw1+2*n-1 ) = d( n )
439  isplit( nsplit ) = n
440  pivmin = pivmin*safemn
441 *
442 * Compute Gershgorin interval [gl,gu] for entire matrix
443 *
444  gu = d( 1 )
445  gl = d( 1 )
446  tmp1 = zero
447 *
448  DO 40 i = 1, n - 1
449  tmp2 = abs( e( i ) )
450  gu = max( gu, d( i )+tmp1+tmp2 )
451  gl = min( gl, d( i )-tmp1-tmp2 )
452  tmp1 = tmp2
453  40 CONTINUE
454  gu = max( gu, d( n )+tmp1 )
455  gl = min( gl, d( n )-tmp1 )
456  tnorm = max( abs( gl ), abs( gu ) )
457  gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin
458  gu = gu + fudge*tnorm*ulp*n + fudge*pivmin
459 *
460  IF( abstol.LE.zero ) THEN
461  atoli = ulp*tnorm
462  ELSE
463  atoli = abstol
464  END IF
465 *
466 * Find out if on an IEEE machine, the sign bit is the
467 * 32nd bit (Big Endian) or the 64th bit (Little Endian)
468 *
469  IF( irange.EQ.1 .OR. nsplit.EQ.1 ) THEN
470  CALL pdlasnbt( ieflag )
471  ELSE
472  ieflag = 0
473  END IF
474  lextra = 0
475  rextra = 0
476 *
477 * Form Initial Interval containing desired eigenvalues
478 *
479  IF( irange.EQ.1 ) THEN
480  initvl = gl
481  initvu = gu
482  work( 1 ) = gl
483  work( 2 ) = gu
484  iwork( 1 ) = 0
485  iwork( 2 ) = n
486  ifrst = 1
487  ilast = n
488  ELSE IF( irange.EQ.2 ) THEN
489  IF( vl.GT.gl ) THEN
490  IF( ieflag.EQ.0 ) THEN
491  CALL pdlapdct( vl, n, work( indrw1+1 ), pivmin, ifrst )
492  ELSE IF( ieflag.EQ.1 ) THEN
493  CALL pdlaiectb( vl, n, work( indrw1+1 ), ifrst )
494  ELSE
495  CALL pdlaiectl( vl, n, work( indrw1+1 ), ifrst )
496  END IF
497  ifrst = ifrst + 1
498  initvl = vl
499  ELSE
500  initvl = gl
501  ifrst = 1
502  END IF
503  IF( vu.LT.gu ) THEN
504  IF( ieflag.EQ.0 ) THEN
505  CALL pdlapdct( vu, n, work( indrw1+1 ), pivmin, ilast )
506  ELSE IF( ieflag.EQ.1 ) THEN
507  CALL pdlaiectb( vu, n, work( indrw1+1 ), ilast )
508  ELSE
509  CALL pdlaiectl( vu, n, work( indrw1+1 ), ilast )
510  END IF
511  initvu = vu
512  ELSE
513  initvu = gu
514  ilast = n
515  END IF
516  work( 1 ) = initvl
517  work( 2 ) = initvu
518  iwork( 1 ) = ifrst - 1
519  iwork( 2 ) = ilast
520  ELSE IF( irange.EQ.3 ) THEN
521  work( 1 ) = gl
522  work( 2 ) = gu
523  iwork( 1 ) = 0
524  iwork( 2 ) = n
525  iwork( 5 ) = il - 1
526  iwork( 6 ) = iu
527  CALL pdlaebz( 0, n, 2, 1, atoli, reltol, pivmin,
528  $ work( indrw1+1 ), iwork( 5 ), work, iwork, nint,
529  $ lsave, ieflag, iinfo )
530  IF( iinfo.NE.0 ) THEN
531  info = 3
532  GO TO 230
533  END IF
534  IF( nint.GT.1 ) THEN
535  IF( iwork( 5 ).EQ.il-1 ) THEN
536  work( 2 ) = work( 4 )
537  iwork( 2 ) = iwork( 4 )
538  ELSE
539  work( 1 ) = work( 3 )
540  iwork( 1 ) = iwork( 3 )
541  END IF
542  IF( iwork( 1 ).LT.0 .OR. iwork( 1 ).GT.il-1 .OR.
543  $ iwork( 2 ).LE.min( iu-1, iwork( 1 ) ) .OR.
544  $ iwork( 2 ).GT.n ) THEN
545  info = 3
546  GO TO 230
547  END IF
548  END IF
549  lextra = il - 1 - iwork( 1 )
550  rextra = iwork( 2 ) - iu
551  initvl = work( 1 )
552  initvu = work( 2 )
553  ifrst = il
554  ilast = iu
555  END IF
556 * NVL = IFRST - 1
557 * NVU = ILAST
558  gl = initvl
559  gu = initvu
560  ngl = iwork( 1 )
561  ngu = iwork( 2 )
562  im = 0
563  found = 0
564  indriw2 = indriw1 + ngu - ngl
565  iend = 0
566  IF( ifrst.GT.ilast )
567  $ GO TO 100
568  IF( ifrst.EQ.1 .AND. ilast.EQ.n )
569  $ irange = 1
570 *
571 * Find Eigenvalues -- Loop Over Blocks
572 *
573  DO 90 jb = 1, nsplit
574  ioff = iend
575  ibegin = ioff + 1
576  iend = isplit( jb )
577  in = iend - ioff
578  IF( jb.NE.1 ) THEN
579  IF( irange.NE.1 ) THEN
580  found = im
581 *
582 * Find total number of eigenvalues found thus far
583 *
584  CALL igsum2d( onedcontext, 'All', ' ', 1, 1, found, 1,
585  $ -1, -1 )
586  ELSE
587  found = ioff
588  END IF
589  END IF
590 * IF( SELF.GE.P )
591 * $ GO TO 30
592  IF( in.NE.n ) THEN
593 *
594 * Compute Gershgorin interval [gl,gu] for split matrix
595 *
596  gu = d( ibegin )
597  gl = d( ibegin )
598  tmp1 = zero
599 *
600  DO 50 j = ibegin, iend - 1
601  tmp2 = abs( e( j ) )
602  gu = max( gu, d( j )+tmp1+tmp2 )
603  gl = min( gl, d( j )-tmp1-tmp2 )
604  tmp1 = tmp2
605  50 CONTINUE
606 *
607  gu = max( gu, d( iend )+tmp1 )
608  gl = min( gl, d( iend )-tmp1 )
609  bnorm = max( abs( gl ), abs( gu ) )
610  gl = gl - fudge*bnorm*ulp*in - fudge*pivmin
611  gu = gu + fudge*bnorm*ulp*in + fudge*pivmin
612 *
613 * Compute ATOLI for the current submatrix
614 *
615  IF( abstol.LE.zero ) THEN
616  atoli = ulp*bnorm
617  ELSE
618  atoli = abstol
619  END IF
620 *
621  IF( gl.LT.initvl ) THEN
622  gl = initvl
623  IF( ieflag.EQ.0 ) THEN
624  CALL pdlapdct( gl, in, work( indrw1+2*ioff+1 ),
625  $ pivmin, ngl )
626  ELSE IF( ieflag.EQ.1 ) THEN
627  CALL pdlaiectb( gl, in, work( indrw1+2*ioff+1 ), ngl )
628  ELSE
629  CALL pdlaiectl( gl, in, work( indrw1+2*ioff+1 ), ngl )
630  END IF
631  ELSE
632  ngl = 0
633  END IF
634  IF( gu.GT.initvu ) THEN
635  gu = initvu
636  IF( ieflag.EQ.0 ) THEN
637  CALL pdlapdct( gu, in, work( indrw1+2*ioff+1 ),
638  $ pivmin, ngu )
639  ELSE IF( ieflag.EQ.1 ) THEN
640  CALL pdlaiectb( gu, in, work( indrw1+2*ioff+1 ), ngu )
641  ELSE
642  CALL pdlaiectl( gu, in, work( indrw1+2*ioff+1 ), ngu )
643  END IF
644  ELSE
645  ngu = in
646  END IF
647  IF( ngl.GE.ngu )
648  $ GO TO 90
649  work( 1 ) = gl
650  work( 2 ) = gu
651  iwork( 1 ) = ngl
652  iwork( 2 ) = ngu
653  END IF
654  offset = found - ngl
655  blkno = jb
656 *
657 * Do a static partitioning of work so that each process
658 * has to find an (almost) equal number of eigenvalues
659 *
660  ncmp = ngu - ngl
661  iload = ncmp / p
662  irem = ncmp - iload*p
663  itmp1 = mod( self-found, p )
664  IF( itmp1.LT.0 )
665  $ itmp1 = itmp1 + p
666  IF( itmp1.LT.irem ) THEN
667  imyload = iload + 1
668  ELSE
669  imyload = iload
670  END IF
671  IF( imyload.EQ.0 ) THEN
672  GO TO 90
673  ELSE IF( in.EQ.1 ) THEN
674  work( indrw2+im+1 ) = work( indrw1+2*ioff+1 )
675  iwork( indriw1+im+1 ) = blkno
676  iwork( indriw2+im+1 ) = offset + 1
677  im = im + 1
678  GO TO 90
679  ELSE
680  inxtload = iload
681  itmp2 = mod( self+1-found, p )
682  IF( itmp2.LT.0 )
683  $ itmp2 = itmp2 + p
684  IF( itmp2.LT.irem )
685  $ inxtload = inxtload + 1
686  lreq = ngl + itmp1*iload + min( irem, itmp1 )
687  rreq = lreq + imyload
688  iwork( 5 ) = lreq
689  iwork( 6 ) = rreq
690  tmp1 = work( 1 )
691  itmp1 = iwork( 1 )
692  CALL pdlaebz( 1, in, 1, 1, atoli, reltol, pivmin,
693  $ work( indrw1+2*ioff+1 ), iwork( 5 ), work,
694  $ iwork, nint, lsave, ieflag, iinfo )
695  alpha = work( 1 )
696  beta = work( 2 )
697  nalpha = iwork( 1 )
698  nbeta = iwork( 2 )
699  dsend = beta
700  IF( nbeta.GT.rreq+inxtload ) THEN
701  nbeta = rreq
702  dsend = alpha
703  END IF
704  last = mod( found+min( ngu-ngl, p )-1, p )
705  IF( last.LT.0 )
706  $ last = last + p
707  IF( self.NE.last ) THEN
708  CALL dgesd2d( onedcontext, 1, 1, dsend, 1, 0, next )
709  CALL igesd2d( onedcontext, 1, 1, nbeta, 1, 0, next )
710  END IF
711  IF( self.NE.mod( found, p ) ) THEN
712  CALL dgerv2d( onedcontext, 1, 1, drecv, 1, 0, prev )
713  CALL igerv2d( onedcontext, 1, 1, irecv, 1, 0, prev )
714  ELSE
715  drecv = tmp1
716  irecv = itmp1
717  END IF
718  work( 1 ) = max( lsave, drecv )
719  iwork( 1 ) = irecv
720  alpha = max( alpha, work( 1 ) )
721  nalpha = max( nalpha, irecv )
722  IF( beta-alpha.LE.max( atoli, reltol*max( abs( alpha ),
723  $ abs( beta ) ) ) ) THEN
724  mid = half*( alpha+beta )
725  DO 60 j = offset + nalpha + 1, offset + nbeta
726  work( indrw2+im+1 ) = mid
727  iwork( indriw1+im+1 ) = blkno
728  iwork( indriw2+im+1 ) = j
729  im = im + 1
730  60 CONTINUE
731  work( 2 ) = alpha
732  iwork( 2 ) = nalpha
733  END IF
734  END IF
735  neigint = iwork( 2 ) - iwork( 1 )
736  IF( neigint.LE.0 )
737  $ GO TO 90
738 *
739 * Call the main computational routine
740 *
741  CALL pdlaebz( 2, in, neigint, 1, atoli, reltol, pivmin,
742  $ work( indrw1+2*ioff+1 ), iwork, work, iwork,
743  $ iout, lsave, ieflag, iinfo )
744  IF( iinfo.NE.0 ) THEN
745  info = 1
746  END IF
747  DO 80 i = 1, iout
748  mid = half*( work( 2*i-1 )+work( 2*i ) )
749  IF( i.GT.iout-iinfo )
750  $ blkno = -blkno
751  DO 70 j = offset + iwork( 2*i-1 ) + 1,
752  $ offset + iwork( 2*i )
753  work( indrw2+im+1 ) = mid
754  iwork( indriw1+im+1 ) = blkno
755  iwork( indriw2+im+1 ) = j
756  im = im + 1
757  70 CONTINUE
758  80 CONTINUE
759  90 CONTINUE
760 *
761 * Find out total number of eigenvalues computed
762 *
763  100 CONTINUE
764  m = im
765  CALL igsum2d( onedcontext, 'ALL', ' ', 1, 1, m, 1, -1, -1 )
766 *
767 * Move the eigenvalues found to their final destinations
768 *
769  DO 130 i = 1, p
770  IF( self.EQ.i-1 ) THEN
771  CALL igebs2d( onedcontext, 'ALL', ' ', 1, 1, im, 1 )
772  IF( im.NE.0 ) THEN
773  CALL igebs2d( onedcontext, 'ALL', ' ', im, 1,
774  $ iwork( indriw2+1 ), im )
775  CALL dgebs2d( onedcontext, 'ALL', ' ', im, 1,
776  $ work( indrw2+1 ), im )
777  CALL igebs2d( onedcontext, 'ALL', ' ', im, 1,
778  $ iwork( indriw1+1 ), im )
779  DO 110 j = 1, im
780  w( iwork( indriw2+j ) ) = work( indrw2+j )
781  iblock( iwork( indriw2+j ) ) = iwork( indriw1+j )
782  110 CONTINUE
783  END IF
784  ELSE
785  CALL igebr2d( onedcontext, 'ALL', ' ', 1, 1, torecv, 1, 0,
786  $ i-1 )
787  IF( torecv.NE.0 ) THEN
788  CALL igebr2d( onedcontext, 'ALL', ' ', torecv, 1, iwork,
789  $ torecv, 0, i-1 )
790  CALL dgebr2d( onedcontext, 'ALL', ' ', torecv, 1, work,
791  $ torecv, 0, i-1 )
792  CALL igebr2d( onedcontext, 'ALL', ' ', torecv, 1,
793  $ iwork( n+1 ), torecv, 0, i-1 )
794  DO 120 j = 1, torecv
795  w( iwork( j ) ) = work( j )
796  iblock( iwork( j ) ) = iwork( n+j )
797  120 CONTINUE
798  END IF
799  END IF
800  130 CONTINUE
801  IF( nsplit.GT.1 .AND. iorder.EQ.1 ) THEN
802 *
803 * Sort the eigenvalues
804 *
805 *
806  DO 140 i = 1, m
807  iwork( m+i ) = i
808  140 CONTINUE
809  CALL dlasrt2( 'I', m, w, iwork( m+1 ), iinfo )
810  DO 150 i = 1, m
811  iwork( i ) = iblock( i )
812  150 CONTINUE
813  DO 160 i = 1, m
814  iblock( i ) = iwork( iwork( m+i ) )
815  160 CONTINUE
816  END IF
817  IF( irange.EQ.3 .AND. ( lextra.GT.0 .OR. rextra.GT.0 ) ) THEN
818 *
819 * Discard unwanted eigenvalues (occurs only when RANGE = 'I',
820 * and eigenvalues IL, and/or IU are in a cluster)
821 *
822  DO 170 i = 1, m
823  work( i ) = w( i )
824  iwork( i ) = i
825  iwork( m+i ) = i
826  170 CONTINUE
827  DO 190 i = 1, lextra
828  itmp1 = i
829  DO 180 j = i + 1, m
830  IF( work( j ).LT.work( itmp1 ) ) THEN
831  itmp1 = j
832  END IF
833  180 CONTINUE
834  tmp1 = work( i )
835  work( i ) = work( itmp1 )
836  work( itmp1 ) = tmp1
837  iwork( iwork( m+itmp1 ) ) = i
838  iwork( iwork( m+i ) ) = itmp1
839  itmp2 = iwork( m+i )
840  iwork( m+i ) = iwork( m+itmp1 )
841  iwork( m+itmp1 ) = itmp2
842  190 CONTINUE
843  DO 210 i = 1, rextra
844  itmp1 = m - i + 1
845  DO 200 j = m - i, lextra + 1, -1
846  IF( work( j ).GT.work( itmp1 ) ) THEN
847  itmp1 = j
848  END IF
849  200 CONTINUE
850  tmp1 = work( m-i+1 )
851  work( m-i+1 ) = work( itmp1 )
852  work( itmp1 ) = tmp1
853  iwork( iwork( m+itmp1 ) ) = m - i + 1
854  iwork( iwork( 2*m-i+1 ) ) = itmp1
855  itmp2 = iwork( 2*m-i+1 )
856  iwork( 2*m-i+1 ) = iwork( m+itmp1 )
857  iwork( m+itmp1 ) = itmp2
858 * IWORK( ITMP1 ) = 1
859  210 CONTINUE
860  j = 0
861  DO 220 i = 1, m
862  IF( iwork( i ).GT.lextra .AND. iwork( i ).LE.m-rextra ) THEN
863  j = j + 1
864  w( j ) = work( iwork( i ) )
865  iblock( j ) = iblock( i )
866  END IF
867  220 CONTINUE
868  m = m - lextra - rextra
869  END IF
870  IF( m.NE.ilast-ifrst+1 ) THEN
871  info = 2
872  END IF
873 *
874  230 CONTINUE
875  CALL blacs_freebuff( onedcontext, 1 )
876  CALL blacs_gridexit( onedcontext )
877  RETURN
878 *
879 * End of PDSTEBZ
880 *
881  END
882 *
883  SUBROUTINE pdlaebz( IJOB, N, MMAX, MINP, ABSTOL, RELTOL, PIVMIN,
884  $ D, NVAL, INTVL, INTVLCT, MOUT, LSAVE, IEFLAG,
885  $ INFO )
886 *
887 * -- ScaLAPACK routine (version 1.7) --
888 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
889 * and University of California, Berkeley.
890 * November 15, 1997
891 *
892 *
893 * .. Scalar Arguments ..
894  INTEGER IEFLAG, IJOB, INFO, MINP, MMAX, MOUT, N
895  DOUBLE PRECISION ABSTOL, LSAVE, PIVMIN, RELTOL
896 * ..
897 * .. Array Arguments ..
898  INTEGER INTVLCT( * ), NVAL( * )
899  DOUBLE PRECISION D( * ), INTVL( * )
900 * ..
901 *
902 * Purpose
903 * =======
904 *
905 * PDLAEBZ contains the iteration loop which computes the eigenvalues
906 * contained in the input intervals [ INTVL(2*j-1), INTVL(2*j) ] where
907 * j = 1,...,MINP. It uses and computes the function N(w), which is
908 * the count of eigenvalues of a symmetric tridiagonal matrix less than
909 * or equal to its argument w.
910 *
911 * This is a ScaLAPACK internal subroutine and arguments are not
912 * checked for unreasonable values.
913 *
914 * Arguments
915 * =========
916 *
917 * IJOB (input) INTEGER
918 * Specifies the computation done by PDLAEBZ
919 * = 0 : Find an interval with desired values of N(w) at the
920 * endpoints of the interval.
921 * = 1 : Find a floating point number contained in the initial
922 * interval with a desired value of N(w).
923 * = 2 : Perform bisection iteration to find eigenvalues of T.
924 *
925 * N (input) INTEGER
926 * The order of the tridiagonal matrix T. N >= 1.
927 *
928 * MMAX (input) INTEGER
929 * The maximum number of intervals that may be generated. If
930 * more than MMAX intervals are generated, then PDLAEBZ will
931 * quit with INFO = MMAX+1.
932 *
933 * MINP (input) INTEGER
934 * The initial number of intervals. MINP <= MMAX.
935 *
936 * ABSTOL (input) DOUBLE PRECISION
937 * The minimum (absolute) width of an interval. When an interval
938 * is narrower than ABSTOL, or than RELTOL times the larger (in
939 * magnitude) endpoint, then it is considered to be sufficiently
940 * small, i.e., converged.
941 * This must be at least zero.
942 *
943 * RELTOL (input) DOUBLE PRECISION
944 * The minimum relative width of an interval. When an interval
945 * is narrower than ABSTOL, or than RELTOL times the larger (in
946 * magnitude) endpoint, then it is considered to be sufficiently
947 * small, i.e., converged.
948 * Note : This should be at least radix*machine epsilon.
949 *
950 * PIVMIN (input) DOUBLE PRECISION
951 * The minimum absolute of a "pivot" in the "paranoid"
952 * implementation of the Sturm sequence loop. This must be at
953 * least max_j |e(j)^2| *safe_min, and at least safe_min, where
954 * safe_min is at least the smallest number that can divide 1.0
955 * without overflow.
956 * See PDLAPDCT for the "paranoid" implementation of the Sturm
957 * sequence loop.
958 *
959 * D (input) DOUBLE PRECISION array, dimension (2*N - 1)
960 * Contains the diagonals and the squares of the off-diagonal
961 * elements of the tridiagonal matrix T. These elements are
962 * assumed to be interleaved in memory for better cache
963 * performance. The diagonal entries of T are in the entries
964 * D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
965 * entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
966 * matrix must be scaled so that its largest entry is no greater
967 * than overflow**(1/2) * underflow**(1/4) in absolute value,
968 * and for greatest accuracy, it should not be much smaller
969 * than that.
970 *
971 * NVAL (input/output) INTEGER array, dimension (4)
972 * If IJOB = 0, the desired values of N(w) are in NVAL(1) and
973 * NVAL(2).
974 * If IJOB = 1, NVAL(2) is the desired value of N(w).
975 * If IJOB = 2, not referenced.
976 * This array will, in general, be reordered on output.
977 *
978 * INTVL (input/output) DOUBLE PRECISION array, dimension (2*MMAX)
979 * The endpoints of the intervals. INTVL(2*j-1) is the left
980 * endpoint of the j-th interval, and INTVL(2*j) is the right
981 * endpoint of the j-th interval. The input intervals will,
982 * in general, be modified, split and reordered by the
983 * calculation.
984 * On input, INTVL contains the MINP input intervals.
985 * On output, INTVL contains the converged intervals.
986 *
987 * INTVLCT (input/output) INTEGER array, dimension (2*MMAX)
988 * The counts at the endpoints of the intervals. INTVLCT(2*j-1)
989 * is the count at the left endpoint of the j-th interval, i.e.,
990 * the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
991 * count at the right endpoint of the j-th interval.
992 * On input, INTVLCT contains the counts at the endpoints of
993 * the MINP input intervals.
994 * On output, INTVLCT contains the counts at the endpoints of
995 * the converged intervals.
996 *
997 * MOUT (output) INTEGER
998 * The number of intervals output.
999 *
1000 * LSAVE (output) DOUBLE PRECISION
1001 * If IJOB = 0 or 2, not referenced.
1002 * If IJOB = 1, this is the largest floating point number
1003 * encountered which has count N(w) = NVAL(1).
1004 *
1005 * IEFLAG (input) INTEGER
1006 * A flag which indicates whether N(w) should be speeded up by
1007 * exploiting IEEE Arithmetic.
1008 *
1009 * INFO (output) INTEGER
1010 * = 0 : All intervals converged.
1011 * = 1 - MMAX : The last INFO intervals did not converge.
1012 * = MMAX + 1 : More than MMAX intervals were generated.
1013 *
1014 * =====================================================================
1015 *
1016 * .. Intrinsic Functions ..
1017  INTRINSIC abs, int, log, max, min
1018 * ..
1019 * .. External Subroutines ..
1020  EXTERNAL pdlaecv, pdlaiectb, pdlaiectl, pdlapdct
1021 * ..
1022 * .. Parameters ..
1023  DOUBLE PRECISION ZERO, TWO, HALF
1024  parameter( zero = 0.0d+0, two = 2.0d+0,
1025  $ half = 1.0d+0 / two )
1026 * ..
1027 * .. Local Scalars ..
1028  INTEGER I, ITMAX, J, K, KF, KL, KLNEW, L, LCNT, LREQ,
1029  $ NALPHA, NBETA, NMID, RCNT, RREQ
1030  DOUBLE PRECISION ALPHA, BETA, MID
1031 * ..
1032 * .. Executable Statements ..
1033 *
1034  kf = 1
1035  kl = minp + 1
1036  info = 0
1037  IF( intvl( 2 )-intvl( 1 ).LE.zero ) THEN
1038  info = minp
1039  mout = kf
1040  RETURN
1041  END IF
1042  IF( ijob.EQ.0 ) THEN
1043 *
1044 * Check if some input intervals have "converged"
1045 *
1046  CALL pdlaecv( 0, kf, kl, intvl, intvlct, nval,
1047  $ max( abstol, pivmin ), reltol )
1048  IF( kf.GE.kl )
1049  $ GO TO 60
1050 *
1051 * Compute upper bound on number of iterations needed
1052 *
1053  itmax = int( ( log( intvl( 2 )-intvl( 1 )+pivmin )-
1054  $ log( pivmin ) ) / log( two ) ) + 2
1055 *
1056 * Iteration Loop
1057 *
1058  DO 20 i = 1, itmax
1059  klnew = kl
1060  DO 10 j = kf, kl - 1
1061  k = 2*j
1062 *
1063 * Bisect the interval and find the count at that point
1064 *
1065  mid = half*( intvl( k-1 )+intvl( k ) )
1066  IF( ieflag.EQ.0 ) THEN
1067  CALL pdlapdct( mid, n, d, pivmin, nmid )
1068  ELSE IF( ieflag.EQ.1 ) THEN
1069  CALL pdlaiectb( mid, n, d, nmid )
1070  ELSE
1071  CALL pdlaiectl( mid, n, d, nmid )
1072  END IF
1073  lreq = nval( k-1 )
1074  rreq = nval( k )
1075  IF( kl.EQ.1 )
1076  $ nmid = min( intvlct( k ),
1077  $ max( intvlct( k-1 ), nmid ) )
1078  IF( nmid.LE.nval( k-1 ) ) THEN
1079  intvl( k-1 ) = mid
1080  intvlct( k-1 ) = nmid
1081  END IF
1082  IF( nmid.GE.nval( k ) ) THEN
1083  intvl( k ) = mid
1084  intvlct( k ) = nmid
1085  END IF
1086  IF( nmid.GT.lreq .AND. nmid.LT.rreq ) THEN
1087  l = 2*klnew
1088  intvl( l-1 ) = mid
1089  intvl( l ) = intvl( k )
1090  intvlct( l-1 ) = nval( k )
1091  intvlct( l ) = intvlct( k )
1092  intvl( k ) = mid
1093  intvlct( k ) = nval( k-1 )
1094  nval( l-1 ) = nval( k )
1095  nval( l ) = nval( l-1 )
1096  nval( k ) = nval( k-1 )
1097  klnew = klnew + 1
1098  END IF
1099  10 CONTINUE
1100  kl = klnew
1101  CALL pdlaecv( 0, kf, kl, intvl, intvlct, nval,
1102  $ max( abstol, pivmin ), reltol )
1103  IF( kf.GE.kl )
1104  $ GO TO 60
1105  20 CONTINUE
1106  ELSE IF( ijob.EQ.1 ) THEN
1107  alpha = intvl( 1 )
1108  beta = intvl( 2 )
1109  nalpha = intvlct( 1 )
1110  nbeta = intvlct( 2 )
1111  lsave = alpha
1112  lreq = nval( 1 )
1113  rreq = nval( 2 )
1114  30 CONTINUE
1115  IF( nbeta.NE.rreq .AND. beta-alpha.GT.
1116  $ max( abstol, reltol*max( abs( alpha ), abs( beta ) ) ) )
1117  $ THEN
1118 *
1119 * Bisect the interval and find the count at that point
1120 *
1121  mid = half*( alpha+beta )
1122  IF( ieflag.EQ.0 ) THEN
1123  CALL pdlapdct( mid, n, d, pivmin, nmid )
1124  ELSE IF( ieflag.EQ.1 ) THEN
1125  CALL pdlaiectb( mid, n, d, nmid )
1126  ELSE
1127  CALL pdlaiectl( mid, n, d, nmid )
1128  END IF
1129  nmid = min( nbeta, max( nalpha, nmid ) )
1130  IF( nmid.GE.rreq ) THEN
1131  beta = mid
1132  nbeta = nmid
1133  ELSE
1134  alpha = mid
1135  nalpha = nmid
1136  IF( nmid.EQ.lreq )
1137  $ lsave = alpha
1138  END IF
1139  GO TO 30
1140  END IF
1141  kl = kf
1142  intvl( 1 ) = alpha
1143  intvl( 2 ) = beta
1144  intvlct( 1 ) = nalpha
1145  intvlct( 2 ) = nbeta
1146  ELSE IF( ijob.EQ.2 ) THEN
1147 *
1148 * Check if some input intervals have "converged"
1149 *
1150  CALL pdlaecv( 1, kf, kl, intvl, intvlct, nval,
1151  $ max( abstol, pivmin ), reltol )
1152  IF( kf.GE.kl )
1153  $ GO TO 60
1154 *
1155 * Compute upper bound on number of iterations needed
1156 *
1157  itmax = int( ( log( intvl( 2 )-intvl( 1 )+pivmin )-
1158  $ log( pivmin ) ) / log( two ) ) + 2
1159 *
1160 * Iteration Loop
1161 *
1162  DO 50 i = 1, itmax
1163  klnew = kl
1164  DO 40 j = kf, kl - 1
1165  k = 2*j
1166  mid = half*( intvl( k-1 )+intvl( k ) )
1167  IF( ieflag.EQ.0 ) THEN
1168  CALL pdlapdct( mid, n, d, pivmin, nmid )
1169  ELSE IF( ieflag.EQ.1 ) THEN
1170  CALL pdlaiectb( mid, n, d, nmid )
1171  ELSE
1172  CALL pdlaiectl( mid, n, d, nmid )
1173  END IF
1174  lcnt = intvlct( k-1 )
1175  rcnt = intvlct( k )
1176  nmid = min( rcnt, max( lcnt, nmid ) )
1177 *
1178 * Form New Interval(s)
1179 *
1180  IF( nmid.EQ.lcnt ) THEN
1181  intvl( k-1 ) = mid
1182  ELSE IF( nmid.EQ.rcnt ) THEN
1183  intvl( k ) = mid
1184  ELSE IF( klnew.LT.mmax+1 ) THEN
1185  l = 2*klnew
1186  intvl( l-1 ) = mid
1187  intvl( l ) = intvl( k )
1188  intvlct( l-1 ) = nmid
1189  intvlct( l ) = intvlct( k )
1190  intvl( k ) = mid
1191  intvlct( k ) = nmid
1192  klnew = klnew + 1
1193  ELSE
1194  info = mmax + 1
1195  RETURN
1196  END IF
1197  40 CONTINUE
1198  kl = klnew
1199  CALL pdlaecv( 1, kf, kl, intvl, intvlct, nval,
1200  $ max( abstol, pivmin ), reltol )
1201  IF( kf.GE.kl )
1202  $ GO TO 60
1203  50 CONTINUE
1204  END IF
1205  60 CONTINUE
1206  info = max( kl-kf, 0 )
1207  mout = kl - 1
1208  RETURN
1209 *
1210 * End of PDLAEBZ
1211 *
1212  END
1213 *
1214 *
1215  SUBROUTINE pdlaecv( IJOB, KF, KL, INTVL, INTVLCT, NVAL, ABSTOL,
1216  $ RELTOL )
1217 *
1218 * -- ScaLAPACK routine (version 1.7) --
1219 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
1220 * and University of California, Berkeley.
1221 * November 15, 1997
1222 *
1223 *
1224 * .. Scalar Arguments ..
1225  INTEGER IJOB, KF, KL
1226  DOUBLE PRECISION ABSTOL, RELTOL
1227 * ..
1228 * .. Array Arguments ..
1229  INTEGER INTVLCT( * ), NVAL( * )
1230  DOUBLE PRECISION INTVL( * )
1231 * ..
1232 *
1233 * Purpose
1234 * =======
1235 *
1236 * PDLAECV checks if the input intervals [ INTVL(2*i-1), INTVL(2*i) ],
1237 * i = KF, ... , KL-1, have "converged".
1238 * PDLAECV modifies KF to be the index of the last converged interval,
1239 * i.e., on output, all intervals [ INTVL(2*i-1), INTVL(2*i) ], i < KF,
1240 * have converged. Note that the input intervals may be reordered by
1241 * PDLAECV.
1242 *
1243 * This is a SCALAPACK internal procedure and arguments are not checked
1244 * for unreasonable values.
1245 *
1246 * Arguments
1247 * =========
1248 *
1249 * IJOB (input) INTEGER
1250 * Specifies the criterion for "convergence" of an interval.
1251 * = 0 : When an interval is narrower than ABSTOL, or than
1252 * RELTOL times the larger (in magnitude) endpoint, then
1253 * it is considered to have "converged".
1254 * = 1 : When an interval is narrower than ABSTOL, or than
1255 * RELTOL times the larger (in magnitude) endpoint, or if
1256 * the counts at the endpoints are identical to the counts
1257 * specified by NVAL ( see NVAL ) then the interval is
1258 * considered to have "converged".
1259 *
1260 * KF (input/output) INTEGER
1261 * On input, the index of the first input interval is 2*KF-1.
1262 * On output, the index of the last converged interval
1263 * is 2*KF-3.
1264 *
1265 * KL (input) INTEGER
1266 * The index of the last input interval is 2*KL-3.
1267 *
1268 * INTVL (input/output) DOUBLE PRECISION array, dimension (2*(KL-KF))
1269 * The endpoints of the intervals. INTVL(2*j-1) is the left
1270 * oendpoint f the j-th interval, and INTVL(2*j) is the right
1271 * endpoint of the j-th interval. The input intervals will,
1272 * in general, be reordered on output.
1273 * On input, INTVL contains the KL-KF input intervals.
1274 * On output, INTVL contains the converged intervals, 1 thru'
1275 * KF-1, and the unconverged intervals, KF thru' KL-1.
1276 *
1277 * INTVLCT (input/output) INTEGER array, dimension (2*(KL-KF))
1278 * The counts at the endpoints of the intervals. INTVLCT(2*j-1)
1279 * is the count at the left endpoint of the j-th interval, i.e.,
1280 * the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
1281 * count at the right endpoint of the j-th interval. This array
1282 * will, in general, be reordered on output.
1283 * See the comments in PDLAEBZ for more on the function N(w).
1284 *
1285 * NVAL (input/output) INTEGER array, dimension (2*(KL-KF))
1286 * The desired counts, N(w), at the endpoints of the
1287 * corresponding intervals. This array will, in general,
1288 * be reordered on output.
1289 *
1290 * ABSTOL (input) DOUBLE PRECISION
1291 * The minimum (absolute) width of an interval. When an interval
1292 * is narrower than ABSTOL, or than RELTOL times the larger (in
1293 * magnitude) endpoint, then it is considered to be sufficiently
1294 * small, i.e., converged.
1295 * Note : This must be at least zero.
1296 *
1297 * RELTOL (input) DOUBLE PRECISION
1298 * The minimum relative width of an interval. When an interval
1299 * is narrower than ABSTOL, or than RELTOL times the larger (in
1300 * magnitude) endpoint, then it is considered to be sufficiently
1301 * small, i.e., converged.
1302 * Note : This should be at least radix*machine epsilon.
1303 *
1304 * =====================================================================
1305 *
1306 * .. Intrinsic Functions ..
1307  INTRINSIC abs, max
1308 * ..
1309 * .. Local Scalars ..
1310  LOGICAL CONDN
1311  INTEGER I, ITMP1, ITMP2, J, K, KFNEW
1312  DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4
1313 * ..
1314 * .. Executable Statements ..
1315 *
1316  kfnew = kf
1317  DO 10 i = kf, kl - 1
1318  k = 2*i
1319  tmp3 = intvl( k-1 )
1320  tmp4 = intvl( k )
1321  tmp1 = abs( tmp4-tmp3 )
1322  tmp2 = max( abs( tmp3 ), abs( tmp4 ) )
1323  condn = tmp1.LT.max( abstol, reltol*tmp2 )
1324  IF( ijob.EQ.0 )
1325  $ condn = condn .OR. ( ( intvlct( k-1 ).EQ.nval( k-1 ) ) .AND.
1326  $ intvlct( k ).EQ.nval( k ) )
1327  IF( condn ) THEN
1328  IF( i.GT.kfnew ) THEN
1329 *
1330 * Reorder Intervals
1331 *
1332  j = 2*kfnew
1333  tmp1 = intvl( k-1 )
1334  tmp2 = intvl( k )
1335  itmp1 = intvlct( k-1 )
1336  itmp2 = intvlct( k )
1337  intvl( k-1 ) = intvl( j-1 )
1338  intvl( k ) = intvl( j )
1339  intvlct( k-1 ) = intvlct( j-1 )
1340  intvlct( k ) = intvlct( j )
1341  intvl( j-1 ) = tmp1
1342  intvl( j ) = tmp2
1343  intvlct( j-1 ) = itmp1
1344  intvlct( j ) = itmp2
1345  IF( ijob.EQ.0 ) THEN
1346  itmp1 = nval( k-1 )
1347  nval( k-1 ) = nval( j-1 )
1348  nval( j-1 ) = itmp1
1349  itmp1 = nval( k )
1350  nval( k ) = nval( j )
1351  nval( j ) = itmp1
1352  END IF
1353  END IF
1354  kfnew = kfnew + 1
1355  END IF
1356  10 CONTINUE
1357  kf = kfnew
1358  RETURN
1359 *
1360 * End of PDLAECV
1361 *
1362  END
1363 *
1364  SUBROUTINE pdlapdct( SIGMA, N, D, PIVMIN, COUNT )
1365 *
1366 * -- ScaLAPACK routine (version 1.7) --
1367 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
1368 * and University of California, Berkeley.
1369 * November 15, 1997
1370 *
1371 *
1372 * .. Scalar Arguments ..
1373  INTEGER COUNT, N
1374  DOUBLE PRECISION PIVMIN, SIGMA
1375 * ..
1376 * .. Array Arguments ..
1377  DOUBLE PRECISION D( * )
1378 * ..
1379 *
1380 * Purpose
1381 * =======
1382 *
1383 * PDLAPDCT counts the number of negative eigenvalues of (T - SIGMA I).
1384 * This implementation of the Sturm Sequence loop has conditionals in
1385 * the innermost loop to avoid overflow and determine the sign of a
1386 * floating point number. PDLAPDCT will be referred to as the "paranoid"
1387 * implementation of the Sturm Sequence loop.
1388 *
1389 * This is a SCALAPACK internal procedure and arguments are not checked
1390 * for unreasonable values.
1391 *
1392 * Arguments
1393 * =========
1394 *
1395 * SIGMA (input) DOUBLE PRECISION
1396 * The shift. PDLAPDCT finds the number of eigenvalues of T less
1397 * than or equal to SIGMA.
1398 *
1399 * N (input) INTEGER
1400 * The order of the tridiagonal matrix T. N >= 1.
1401 *
1402 * D (input) DOUBLE PRECISION array, dimension (2*N - 1)
1403 * Contains the diagonals and the squares of the off-diagonal
1404 * elements of the tridiagonal matrix T. These elements are
1405 * assumed to be interleaved in memory for better cache
1406 * performance. The diagonal entries of T are in the entries
1407 * D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
1408 * entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
1409 * matrix must be scaled so that its largest entry is no greater
1410 * than overflow**(1/2) * underflow**(1/4) in absolute value,
1411 * and for greatest accuracy, it should not be much smaller
1412 * than that.
1413 *
1414 * PIVMIN (input) DOUBLE PRECISION
1415 * The minimum absolute of a "pivot" in this "paranoid"
1416 * implementation of the Sturm sequence loop. This must be at
1417 * least max_j |e(j)^2| *safe_min, and at least safe_min, where
1418 * safe_min is at least the smallest number that can divide 1.0
1419 * without overflow.
1420 *
1421 * COUNT (output) INTEGER
1422 * The count of the number of eigenvalues of T less than or
1423 * equal to SIGMA.
1424 *
1425 * =====================================================================
1426 *
1427 * .. Intrinsic Functions ..
1428  INTRINSIC abs
1429 * ..
1430 * .. Parameters ..
1431  DOUBLE PRECISION ZERO
1432  PARAMETER ( ZERO = 0.0d+0 )
1433 * ..
1434 * .. Local Scalars ..
1435  INTEGER I
1436  DOUBLE PRECISION TMP
1437 * ..
1438 * .. Executable Statements ..
1439 *
1440  tmp = d( 1 ) - sigma
1441  IF( abs( tmp ).LE.pivmin )
1442  $ tmp = -pivmin
1443  count = 0
1444  IF( tmp.LE.zero )
1445  $ count = 1
1446  DO 10 i = 3, 2*n - 1, 2
1447  tmp = d( i ) - d( i-1 ) / tmp - sigma
1448  IF( abs( tmp ).LE.pivmin )
1449  $ tmp = -pivmin
1450  IF( tmp.LE.zero )
1451  $ count = count + 1
1452  10 CONTINUE
1453 *
1454  RETURN
1455 *
1456 * End of PDLAPDCT
1457 *
1458  END
globchk
subroutine globchk(ICTXT, N, X, LDX, IWORK, INFO)
Definition: pchkxmat.f:403
max
#define max(A, B)
Definition: pcgemr.c:180
dlasrt2
subroutine dlasrt2(ID, N, D, KEY, INFO)
Definition: dlasrt2.f:4
pdlapdct
subroutine pdlapdct(SIGMA, N, D, PIVMIN, COUNT)
Definition: pdstebz.f:1365
pdlaecv
subroutine pdlaecv(IJOB, KF, KL, INTVL, INTVLCT, NVAL, ABSTOL, RELTOL)
Definition: pdstebz.f:1217
pdstebz
subroutine pdstebz(ICTXT, RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdstebz.f:4
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pdlaebz
subroutine pdlaebz(IJOB, N, MMAX, MINP, ABSTOL, RELTOL, PIVMIN, D, NVAL, INTVL, INTVLCT, MOUT, LSAVE, IEFLAG, INFO)
Definition: pdstebz.f:886
min
#define min(A, B)
Definition: pcgemr.c:181