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