SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ psstebz()

subroutine psstebz ( integer  ictxt,
character  range,
character  order,
integer  n,
real  vl,
real  vu,
integer  il,
integer  iu,
real  abstol,
real, dimension( * )  d,
real, dimension( * )  e,
integer  m,
integer  nsplit,
real, dimension( * )  w,
integer, dimension( * )  iblock,
integer, dimension( * )  isplit,
real, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer  liwork,
integer  info 
)

Definition at line 1 of file psstebz.f.

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*
pslamch
real function pslamch(ictxt, cmach)
Definition pcblastst.f:7455
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
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
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
Here is the call graph for this function:
Here is the caller graph for this function: