SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 REAL ALPHA, ATOLI, BETA, BNORM, DRECV, DSEND, GL,
249 $ GU, INITVL, INITVU, LSAVE, MID, PIVMIN, RELTOL,
250 $ safemn, tmp1, tmp2, tnorm, ulp
251* ..
252* .. Local Arrays ..
253 INTEGER IDUM( 5, 2 )
254 INTEGER TORECV( 1, 1 )
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( 1, 1 ).NE.0 ) THEN
778 CALL igebr2d( onedcontext, 'ALL', ' ', torecv( 1, 1 ), 1,
779 $ iwork, torecv( 1, 1 ), 0, i-1 )
780 CALL sgebr2d( onedcontext, 'ALL', ' ', torecv( 1, 1 ), 1,
781 $ work, torecv( 1, 1 ), 0, i-1 )
782 CALL igebr2d( onedcontext, 'ALL', ' ', torecv( 1, 1 ), 1,
783 $ iwork( n+1 ), torecv( 1, 1 ), 0, i-1 )
784 DO 120 j = 1, torecv( 1, 1 )
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
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
globchk
subroutine globchk(ictxt, n, x, ldx, iwork, info)
Definition pchkxmat.f:403
pslaecv
subroutine pslaecv(ijob, kf, kl, intvl, intvlct, nval, abstol, reltol)
Definition psstebz.f:1201
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
pslaebz
subroutine pslaebz(ijob, n, mmax, minp, abstol, reltol, pivmin, d, nval, intvl, intvlct, mout, lsave, ieflag, info)
Definition psstebz.f:876
pslapdct
subroutine pslapdct(sigma, n, d, pivmin, count)
Definition psstebz.f:1349
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2
slasrt2
subroutine slasrt2(id, n, d, key, info)
Definition slasrt2.f:4