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

◆ pdstebz()

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

Definition at line 1 of file pdstebz.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 DOUBLE PRECISION ABSTOL, VL, VU
15* ..
16* .. Array Arguments ..
17 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
18 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
19* ..
20*
21* Purpose
22* =======
23*
24* PDSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix in
25* parallel. The user may ask for all eigenvalues, all eigenvalues in
26* the interval [VL, VU], or the eigenvalues indexed IL through IU. A
27* static partitioning of work is done at the beginning of PDSTEBZ which
28* results in all processes finding an (almost) equal number of
29* eigenvalues.
30*
31* NOTE : It is assumed that the user is on an IEEE machine. If the user
32* is not on an IEEE mchine, set the compile time flag NO_IEEE
33* to 1 (in SLmake.inc). The features of IEEE arithmetic that
34* are needed for the "fast" Sturm Count are : (a) infinity
35* arithmetic (b) the sign bit of a single precision floating
36* point number is assumed be in the 32nd bit position
37* (c) the sign of negative zero.
38*
39* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
40* Matrix", Report CS41, Computer Science Dept., Stanford
41* University, July 21, 1966.
42*
43* Arguments
44* =========
45*
46* ICTXT (global input) INTEGER
47* The BLACS context handle.
48*
49* RANGE (global input) CHARACTER
50* Specifies which eigenvalues are to be found.
51* = 'A': ("All") all eigenvalues will be found.
52* = 'V': ("Value") all eigenvalues in the interval
53* [VL, VU] will be found.
54* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
55* entire matrix) will be found.
56*
57* ORDER (global input) CHARACTER
58* Specifies the order in which the eigenvalues and their block
59* numbers are stored in W and IBLOCK.
60* = 'B': ("By Block") the eigenvalues will be grouped by
61* split-off block (see IBLOCK, ISPLIT) and
62* ordered from smallest to largest within
63* the block.
64* = 'E': ("Entire matrix")
65* the eigenvalues for the entire matrix
66* will be ordered from smallest to largest.
67*
68* N (global input) INTEGER
69* The order of the tridiagonal matrix T. N >= 0.
70*
71* VL (global input) DOUBLE PRECISION
72* If RANGE='V', the lower bound of the interval to be searched
73* for eigenvalues. Eigenvalues less than VL will not be
74* returned. Not referenced if RANGE='A' or 'I'.
75*
76* VU (global input) DOUBLE PRECISION
77* If RANGE='V', the upper bound of the interval to be searched
78* for eigenvalues. Eigenvalues greater than VU will not be
79* returned. VU must be greater than VL. Not referenced if
80* RANGE='A' or 'I'.
81*
82* IL (global input) INTEGER
83* If RANGE='I', the index (from smallest to largest) of the
84* smallest eigenvalue to be returned. IL must be at least 1.
85* Not referenced if RANGE='A' or 'V'.
86*
87* IU (global input) INTEGER
88* If RANGE='I', the index (from smallest to largest) of the
89* largest eigenvalue to be returned. IU must be at least IL
90* and no greater than N. Not referenced if RANGE='A' or 'V'.
91*
92* ABSTOL (global input) DOUBLE PRECISION
93* The absolute tolerance for the eigenvalues. An eigenvalue
94* (or cluster) is considered to be located if it has been
95* determined to lie in an interval whose width is ABSTOL or
96* less. If ABSTOL is less than or equal to zero, then ULP*|T|
97* will be used, where |T| means the 1-norm of T.
98* Eigenvalues will be computed most accurately when ABSTOL is
99* set to the underflow threshold DLAMCH('U'), not zero.
100* Note : If eigenvectors are desired later by inverse iteration
101* ( PDSTEIN ), ABSTOL should be set to 2*PDLAMCH('S').
102*
103* D (global input) DOUBLE PRECISION array, dimension (N)
104* The n diagonal elements of the tridiagonal matrix T. To
105* avoid overflow, the matrix must be scaled so that its largest
106* entry is no greater than overflow**(1/2) * underflow**(1/4)
107* in absolute value, and for greatest accuracy, it should not
108* be much smaller than that.
109*
110* E (global input) DOUBLE PRECISION array, dimension (N-1)
111* The (n-1) off-diagonal elements of the tridiagonal matrix T.
112* To avoid overflow, the matrix must be scaled so that its
113* largest entry is no greater than overflow**(1/2) *
114* underflow**(1/4) in absolute value, and for greatest
115* accuracy, it should not be much smaller than that.
116*
117* M (global output) INTEGER
118* The actual number of eigenvalues found. 0 <= M <= N.
119* (See also the description of INFO=2)
120*
121* NSPLIT (global output) INTEGER
122* The number of diagonal blocks in the matrix T.
123* 1 <= NSPLIT <= N.
124*
125* W (global output) DOUBLE PRECISION array, dimension (N)
126* On exit, the first M elements of W contain the eigenvalues
127* on all processes.
128*
129* IBLOCK (global output) INTEGER array, dimension (N)
130* At each row/column j where E(j) is zero or small, the
131* matrix T is considered to split into a block diagonal
132* matrix. On exit IBLOCK(i) specifies which block (from 1
133* to the number of blocks) the eigenvalue W(i) belongs to.
134* NOTE: in the (theoretically impossible) event that bisection
135* does not converge for some or all eigenvalues, INFO is set
136* to 1 and the ones for which it did not are identified by a
137* negative block number.
138*
139* ISPLIT (global output) INTEGER array, dimension (N)
140* The splitting points, at which T breaks up into submatrices.
141* The first submatrix consists of rows/columns 1 to ISPLIT(1),
142* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
143* etc., and the NSPLIT-th consists of rows/columns
144* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
145* (Only the first NSPLIT elements will actually be used, but
146* since the user cannot know a priori what value NSPLIT will
147* have, N words must be reserved for ISPLIT.)
148*
149* WORK (local workspace) DOUBLE PRECISION array,
150* dimension ( MAX( 5*N, 7 ) )
151*
152* LWORK (local input) INTEGER
153* size of array WORK must be >= MAX( 5*N, 7 )
154* If LWORK = -1, then LWORK is global input and a workspace
155* query is assumed; the routine only calculates the minimum
156* and optimal size for all work arrays. Each of these
157* values is returned in the first entry of the corresponding
158* work array, and no error message is issued by PXERBLA.
159*
160* IWORK (local workspace) INTEGER array, dimension ( MAX( 4*N, 14 ) )
161*
162* LIWORK (local input) INTEGER
163* size of array IWORK must be >= MAX( 4*N, 14, NPROCS )
164* If LIWORK = -1, then LIWORK is global input and a workspace
165* query is assumed; the routine only calculates the minimum
166* and optimal size for all work arrays. Each of these
167* values is returned in the first entry of the corresponding
168* work array, and no error message is issued by PXERBLA.
169*
170* INFO (global output) INTEGER
171* = 0 : successful exit
172* < 0 : if INFO = -i, the i-th argument had an illegal value
173* > 0 : some or all of the eigenvalues failed to converge or
174* were not computed:
175* = 1 : Bisection failed to converge for some eigenvalues;
176* these eigenvalues are flagged by a negative block
177* number. The effect is that the eigenvalues may not
178* be as accurate as the absolute and relative
179* tolerances. This is generally caused by arithmetic
180* which is less accurate than PDLAMCH says.
181* = 2 : There is a mismatch between the number of
182* eigenvalues output and the number desired.
183* = 3 : RANGE='i', and the Gershgorin interval initially
184* used was incorrect. No eigenvalues were computed.
185* Probable cause: your machine has sloppy floating
186* point arithmetic.
187* Cure: Increase the PARAMETER "FUDGE", recompile,
188* and try again.
189*
190* Internal Parameters
191* ===================
192*
193* RELFAC DOUBLE PRECISION, default = 2.0
194* The relative tolerance. An interval [a,b] lies within
195* "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
196* where "ulp" is the machine precision (distance from 1 to
197* the next larger floating point number.)
198*
199* FUDGE DOUBLE PRECISION, default = 2.0
200* A "fudge factor" to widen the Gershgorin intervals. Ideally,
201* a value of 1 should work, but on machines with sloppy
202* arithmetic, this needs to be larger. The default for
203* publicly released versions should be large enough to handle
204* the worst machine around. Note that this has no effect
205* on the accuracy of the solution.
206*
207* =====================================================================
208*
209* .. Intrinsic Functions ..
210 INTRINSIC abs, dble, ichar, max, min, mod
211* ..
212* .. External Functions ..
213 LOGICAL LSAME
214 INTEGER BLACS_PNUM
215 DOUBLE PRECISION PDLAMCH
216 EXTERNAL lsame, blacs_pnum, pdlamch
217* ..
218* .. External Subroutines ..
219 EXTERNAL blacs_freebuff, blacs_get, blacs_gridexit,
220 $ blacs_gridinfo, blacs_gridmap, dgebr2d,
221 $ dgebs2d, dgerv2d, dgesd2d, dlasrt2, globchk,
222 $ igebr2d, igebs2d, igerv2d, igesd2d, igsum2d,
223 $ pdlaebz, pdlaiectb, pdlaiectl, pdlapdct,
224 $ pdlasnbt, pxerbla
225* ..
226* .. Parameters ..
227 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
228 $ MB_, NB_, RSRC_, CSRC_, LLD_
229 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
230 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
231 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
232 INTEGER BIGNUM, DESCMULT
233 parameter( bignum = 10000, descmult = 100 )
234 DOUBLE PRECISION ZERO, ONE, TWO, FIVE, HALF
235 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
236 $ five = 5.0d+0, half = 1.0d+0 / two )
237 DOUBLE PRECISION FUDGE, RELFAC
238 parameter( fudge = 2.0d+0, relfac = 2.0d+0 )
239* ..
240* .. Local Scalars ..
241 LOGICAL LQUERY
242 INTEGER BLKNO, FOUND, I, IBEGIN, IEFLAG, IEND, IFRST,
243 $ IINFO, ILAST, ILOAD, IM, IMYLOAD, IN, INDRIW1,
244 $ INDRIW2, INDRW1, INDRW2, INXTLOAD, IOFF,
245 $ IORDER, IOUT, IRANGE, IRECV, IREM, ITMP1,
246 $ ITMP2, J, JB, K, LAST, LEXTRA, LREQ, MYCOL,
247 $ MYROW, NALPHA, NBETA, NCMP, NEIGINT, NEXT, NGL,
248 $ NGLOB, NGU, NINT, NPCOL, NPROW, OFFSET,
249 $ ONEDCONTEXT, P, PREV, REXTRA, RREQ, SELF
250 DOUBLE PRECISION ALPHA, ATOLI, BETA, BNORM, DRECV, DSEND, GL,
251 $ GU, INITVL, INITVU, LSAVE, MID, PIVMIN, RELTOL,
252 $ SAFEMN, TMP1, TMP2, TNORM, ULP
253* ..
254* .. Local Arrays ..
255 INTEGER IDUM( 5, 2 )
256 INTEGER TORECV( 1, 1 )
257* ..
258* .. Executable Statements ..
259* This is just to keep ftnchek happy
260 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
261 $ rsrc_.LT.0 )RETURN
262*
263* Set up process grid
264*
265 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
266*
267 info = 0
268 m = 0
269*
270* Decode RANGE
271*
272 IF( lsame( range, 'A' ) ) THEN
273 irange = 1
274 ELSE IF( lsame( range, 'V' ) ) THEN
275 irange = 2
276 ELSE IF( lsame( range, 'I' ) ) THEN
277 irange = 3
278 ELSE
279 irange = 0
280 END IF
281*
282* Decode ORDER
283*
284 IF( lsame( order, 'B' ) ) THEN
285 iorder = 2
286 ELSE IF( lsame( order, 'E' ) .OR. lsame( order, 'A' ) ) THEN
287 iorder = 1
288 ELSE
289 iorder = 0
290 END IF
291*
292* Check for Errors
293*
294 IF( nprow.EQ.-1 ) THEN
295 info = -1
296 ELSE
297*
298* Get machine constants
299*
300 safemn = pdlamch( ictxt, 'S' )
301 ulp = pdlamch( ictxt, 'P' )
302 reltol = ulp*relfac
303 idum( 1, 1 ) = ichar( range )
304 idum( 1, 2 ) = 2
305 idum( 2, 1 ) = ichar( order )
306 idum( 2, 2 ) = 3
307 idum( 3, 1 ) = n
308 idum( 3, 2 ) = 4
309 nglob = 5
310 IF( irange.EQ.3 ) THEN
311 idum( 4, 1 ) = il
312 idum( 4, 2 ) = 7
313 idum( 5, 1 ) = iu
314 idum( 5, 2 ) = 8
315 ELSE
316 idum( 4, 1 ) = 0
317 idum( 4, 2 ) = 0
318 idum( 5, 1 ) = 0
319 idum( 5, 2 ) = 0
320 END IF
321 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
322 work( 1 ) = abstol
323 IF( irange.EQ.2 ) THEN
324 work( 2 ) = vl
325 work( 3 ) = vu
326 ELSE
327 work( 2 ) = zero
328 work( 3 ) = zero
329 END IF
330 CALL dgebs2d( ictxt, 'ALL', ' ', 3, 1, work, 3 )
331 ELSE
332 CALL dgebr2d( ictxt, 'ALL', ' ', 3, 1, work, 3, 0, 0 )
333 END IF
334 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
335 IF( info.EQ.0 ) THEN
336 IF( irange.EQ.0 ) THEN
337 info = -2
338 ELSE IF( iorder.EQ.0 ) THEN
339 info = -3
340 ELSE IF( irange.EQ.2 .AND. vl.GE.vu ) THEN
341 info = -5
342 ELSE IF( irange.EQ.3 .AND. ( il.LT.1 .OR. il.GT.max( 1,
343 $ n ) ) ) THEN
344 info = -6
345 ELSE IF( irange.EQ.3 .AND. ( iu.LT.min( n,
346 $ il ) .OR. iu.GT.n ) ) THEN
347 info = -7
348 ELSE IF( lwork.LT.max( 5*n, 7 ) .AND. .NOT.lquery ) THEN
349 info = -18
350 ELSE IF( liwork.LT.max( 4*n, 14, nprow*npcol ) .AND. .NOT.
351 $ lquery ) THEN
352 info = -20
353 ELSE IF( irange.EQ.2 .AND. ( abs( work( 2 )-vl ).GT.five*
354 $ ulp*abs( vl ) ) ) THEN
355 info = -5
356 ELSE IF( irange.EQ.2 .AND. ( abs( work( 3 )-vu ).GT.five*
357 $ ulp*abs( vu ) ) ) THEN
358 info = -6
359 ELSE IF( abs( work( 1 )-abstol ).GT.five*ulp*abs( abstol ) )
360 $ THEN
361 info = -9
362 END IF
363 END IF
364 IF( info.EQ.0 )
365 $ info = bignum
366 CALL globchk( ictxt, nglob, idum, 5, iwork, info )
367 IF( info.EQ.bignum ) THEN
368 info = 0
369 ELSE IF( mod( info, descmult ).EQ.0 ) THEN
370 info = -info / descmult
371 ELSE
372 info = -info
373 END IF
374 END IF
375 work( 1 ) = dble( max( 5*n, 7 ) )
376 iwork( 1 ) = max( 4*n, 14, nprow*npcol )
377*
378 IF( info.NE.0 ) THEN
379 CALL pxerbla( ictxt, 'PDSTEBZ', -info )
380 RETURN
381 ELSE IF( lwork.EQ.-1 .AND. liwork.EQ.-1 ) THEN
382 RETURN
383 END IF
384*
385*
386* Quick return if possible
387*
388 IF( n.EQ.0 )
389 $ RETURN
390*
391 k = 1
392 DO 20 i = 0, nprow - 1
393 DO 10 j = 0, npcol - 1
394 iwork( k ) = blacs_pnum( ictxt, i, j )
395 k = k + 1
396 10 CONTINUE
397 20 CONTINUE
398*
399 p = nprow*npcol
400 nprow = 1
401 npcol = p
402*
403 CALL blacs_get( ictxt, 10, onedcontext )
404 CALL blacs_gridmap( onedcontext, iwork, nprow, nprow, npcol )
405 CALL blacs_gridinfo( onedcontext, i, j, k, self )
406*
407* Simplifications:
408*
409 IF( irange.EQ.3 .AND. il.EQ.1 .AND. iu.EQ.n )
410 $ irange = 1
411*
412 next = mod( self+1, p )
413 prev = mod( p+self-1, p )
414*
415* Compute squares of off-diagonals, splitting points and pivmin.
416* Interleave diagonals and off-diagonals.
417*
418 indrw1 = max( 2*n, 4 )
419 indrw2 = indrw1 + 2*n
420 indriw1 = max( 2*n, 8 )
421 nsplit = 1
422 work( indrw1+2*n ) = zero
423 pivmin = one
424*
425 DO 30 i = 1, n - 1
426 tmp1 = e( i )**2
427 j = 2*i
428 work( indrw1+j-1 ) = d( i )
429 IF( abs( d( i+1 )*d( i ) )*ulp**2+safemn.GT.tmp1 ) THEN
430 isplit( nsplit ) = i
431 nsplit = nsplit + 1
432 work( indrw1+j ) = zero
433 ELSE
434 work( indrw1+j ) = tmp1
435 pivmin = max( pivmin, tmp1 )
436 END IF
437 30 CONTINUE
438 work( indrw1+2*n-1 ) = d( n )
439 isplit( nsplit ) = n
440 pivmin = pivmin*safemn
441*
442* Compute Gershgorin interval [gl,gu] for entire matrix
443*
444 gu = d( 1 )
445 gl = d( 1 )
446 tmp1 = zero
447*
448 DO 40 i = 1, n - 1
449 tmp2 = abs( e( i ) )
450 gu = max( gu, d( i )+tmp1+tmp2 )
451 gl = min( gl, d( i )-tmp1-tmp2 )
452 tmp1 = tmp2
453 40 CONTINUE
454 gu = max( gu, d( n )+tmp1 )
455 gl = min( gl, d( n )-tmp1 )
456 tnorm = max( abs( gl ), abs( gu ) )
457 gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin
458 gu = gu + fudge*tnorm*ulp*n + fudge*pivmin
459*
460 IF( abstol.LE.zero ) THEN
461 atoli = ulp*tnorm
462 ELSE
463 atoli = abstol
464 END IF
465*
466* Find out if on an IEEE machine, the sign bit is the
467* 32nd bit (Big Endian) or the 64th bit (Little Endian)
468*
469 IF( irange.EQ.1 .OR. nsplit.EQ.1 ) THEN
470 CALL pdlasnbt( ieflag )
471 ELSE
472 ieflag = 0
473 END IF
474 lextra = 0
475 rextra = 0
476*
477* Form Initial Interval containing desired eigenvalues
478*
479 IF( irange.EQ.1 ) THEN
480 initvl = gl
481 initvu = gu
482 work( 1 ) = gl
483 work( 2 ) = gu
484 iwork( 1 ) = 0
485 iwork( 2 ) = n
486 ifrst = 1
487 ilast = n
488 ELSE IF( irange.EQ.2 ) THEN
489 IF( vl.GT.gl ) THEN
490 IF( ieflag.EQ.0 ) THEN
491 CALL pdlapdct( vl, n, work( indrw1+1 ), pivmin, ifrst )
492 ELSE IF( ieflag.EQ.1 ) THEN
493 CALL pdlaiectb( vl, n, work( indrw1+1 ), ifrst )
494 ELSE
495 CALL pdlaiectl( vl, n, work( indrw1+1 ), ifrst )
496 END IF
497 ifrst = ifrst + 1
498 initvl = vl
499 ELSE
500 initvl = gl
501 ifrst = 1
502 END IF
503 IF( vu.LT.gu ) THEN
504 IF( ieflag.EQ.0 ) THEN
505 CALL pdlapdct( vu, n, work( indrw1+1 ), pivmin, ilast )
506 ELSE IF( ieflag.EQ.1 ) THEN
507 CALL pdlaiectb( vu, n, work( indrw1+1 ), ilast )
508 ELSE
509 CALL pdlaiectl( vu, n, work( indrw1+1 ), ilast )
510 END IF
511 initvu = vu
512 ELSE
513 initvu = gu
514 ilast = n
515 END IF
516 work( 1 ) = initvl
517 work( 2 ) = initvu
518 iwork( 1 ) = ifrst - 1
519 iwork( 2 ) = ilast
520 ELSE IF( irange.EQ.3 ) THEN
521 work( 1 ) = gl
522 work( 2 ) = gu
523 iwork( 1 ) = 0
524 iwork( 2 ) = n
525 iwork( 5 ) = il - 1
526 iwork( 6 ) = iu
527 CALL pdlaebz( 0, n, 2, 1, atoli, reltol, pivmin,
528 $ work( indrw1+1 ), iwork( 5 ), work, iwork, nint,
529 $ lsave, ieflag, iinfo )
530 IF( iinfo.NE.0 ) THEN
531 info = 3
532 GO TO 230
533 END IF
534 IF( nint.GT.1 ) THEN
535 IF( iwork( 5 ).EQ.il-1 ) THEN
536 work( 2 ) = work( 4 )
537 iwork( 2 ) = iwork( 4 )
538 ELSE
539 work( 1 ) = work( 3 )
540 iwork( 1 ) = iwork( 3 )
541 END IF
542 IF( iwork( 1 ).LT.0 .OR. iwork( 1 ).GT.il-1 .OR.
543 $ iwork( 2 ).LE.min( iu-1, iwork( 1 ) ) .OR.
544 $ iwork( 2 ).GT.n ) THEN
545 info = 3
546 GO TO 230
547 END IF
548 END IF
549 lextra = il - 1 - iwork( 1 )
550 rextra = iwork( 2 ) - iu
551 initvl = work( 1 )
552 initvu = work( 2 )
553 ifrst = il
554 ilast = iu
555 END IF
556* NVL = IFRST - 1
557* NVU = ILAST
558 gl = initvl
559 gu = initvu
560 ngl = iwork( 1 )
561 ngu = iwork( 2 )
562 im = 0
563 found = 0
564 indriw2 = indriw1 + ngu - ngl
565 iend = 0
566 IF( ifrst.GT.ilast )
567 $ GO TO 100
568 IF( ifrst.EQ.1 .AND. ilast.EQ.n )
569 $ irange = 1
570*
571* Find Eigenvalues -- Loop Over Blocks
572*
573 DO 90 jb = 1, nsplit
574 ioff = iend
575 ibegin = ioff + 1
576 iend = isplit( jb )
577 in = iend - ioff
578 IF( jb.NE.1 ) THEN
579 IF( irange.NE.1 ) THEN
580 found = im
581*
582* Find total number of eigenvalues found thus far
583*
584 CALL igsum2d( onedcontext, 'All', ' ', 1, 1, found, 1,
585 $ -1, -1 )
586 ELSE
587 found = ioff
588 END IF
589 END IF
590* IF( SELF.GE.P )
591* $ GO TO 30
592 IF( in.NE.n ) THEN
593*
594* Compute Gershgorin interval [gl,gu] for split matrix
595*
596 gu = d( ibegin )
597 gl = d( ibegin )
598 tmp1 = zero
599*
600 DO 50 j = ibegin, iend - 1
601 tmp2 = abs( e( j ) )
602 gu = max( gu, d( j )+tmp1+tmp2 )
603 gl = min( gl, d( j )-tmp1-tmp2 )
604 tmp1 = tmp2
605 50 CONTINUE
606*
607 gu = max( gu, d( iend )+tmp1 )
608 gl = min( gl, d( iend )-tmp1 )
609 bnorm = max( abs( gl ), abs( gu ) )
610 gl = gl - fudge*bnorm*ulp*in - fudge*pivmin
611 gu = gu + fudge*bnorm*ulp*in + fudge*pivmin
612*
613* Compute ATOLI for the current submatrix
614*
615 IF( abstol.LE.zero ) THEN
616 atoli = ulp*bnorm
617 ELSE
618 atoli = abstol
619 END IF
620*
621 IF( gl.LT.initvl ) THEN
622 gl = initvl
623 IF( ieflag.EQ.0 ) THEN
624 CALL pdlapdct( gl, in, work( indrw1+2*ioff+1 ),
625 $ pivmin, ngl )
626 ELSE IF( ieflag.EQ.1 ) THEN
627 CALL pdlaiectb( gl, in, work( indrw1+2*ioff+1 ), ngl )
628 ELSE
629 CALL pdlaiectl( gl, in, work( indrw1+2*ioff+1 ), ngl )
630 END IF
631 ELSE
632 ngl = 0
633 END IF
634 IF( gu.GT.initvu ) THEN
635 gu = initvu
636 IF( ieflag.EQ.0 ) THEN
637 CALL pdlapdct( gu, in, work( indrw1+2*ioff+1 ),
638 $ pivmin, ngu )
639 ELSE IF( ieflag.EQ.1 ) THEN
640 CALL pdlaiectb( gu, in, work( indrw1+2*ioff+1 ), ngu )
641 ELSE
642 CALL pdlaiectl( gu, in, work( indrw1+2*ioff+1 ), ngu )
643 END IF
644 ELSE
645 ngu = in
646 END IF
647 IF( ngl.GE.ngu )
648 $ GO TO 90
649 work( 1 ) = gl
650 work( 2 ) = gu
651 iwork( 1 ) = ngl
652 iwork( 2 ) = ngu
653 END IF
654 offset = found - ngl
655 blkno = jb
656*
657* Do a static partitioning of work so that each process
658* has to find an (almost) equal number of eigenvalues
659*
660 ncmp = ngu - ngl
661 iload = ncmp / p
662 irem = ncmp - iload*p
663 itmp1 = mod( self-found, p )
664 IF( itmp1.LT.0 )
665 $ itmp1 = itmp1 + p
666 IF( itmp1.LT.irem ) THEN
667 imyload = iload + 1
668 ELSE
669 imyload = iload
670 END IF
671 IF( imyload.EQ.0 ) THEN
672 GO TO 90
673 ELSE IF( in.EQ.1 ) THEN
674 work( indrw2+im+1 ) = work( indrw1+2*ioff+1 )
675 iwork( indriw1+im+1 ) = blkno
676 iwork( indriw2+im+1 ) = offset + 1
677 im = im + 1
678 GO TO 90
679 ELSE
680 inxtload = iload
681 itmp2 = mod( self+1-found, p )
682 IF( itmp2.LT.0 )
683 $ itmp2 = itmp2 + p
684 IF( itmp2.LT.irem )
685 $ inxtload = inxtload + 1
686 lreq = ngl + itmp1*iload + min( irem, itmp1 )
687 rreq = lreq + imyload
688 iwork( 5 ) = lreq
689 iwork( 6 ) = rreq
690 tmp1 = work( 1 )
691 itmp1 = iwork( 1 )
692 CALL pdlaebz( 1, in, 1, 1, atoli, reltol, pivmin,
693 $ work( indrw1+2*ioff+1 ), iwork( 5 ), work,
694 $ iwork, nint, lsave, ieflag, iinfo )
695 alpha = work( 1 )
696 beta = work( 2 )
697 nalpha = iwork( 1 )
698 nbeta = iwork( 2 )
699 dsend = beta
700 IF( nbeta.GT.rreq+inxtload ) THEN
701 nbeta = rreq
702 dsend = alpha
703 END IF
704 last = mod( found+min( ngu-ngl, p )-1, p )
705 IF( last.LT.0 )
706 $ last = last + p
707 IF( self.NE.last ) THEN
708 CALL dgesd2d( onedcontext, 1, 1, dsend, 1, 0, next )
709 CALL igesd2d( onedcontext, 1, 1, nbeta, 1, 0, next )
710 END IF
711 IF( self.NE.mod( found, p ) ) THEN
712 CALL dgerv2d( onedcontext, 1, 1, drecv, 1, 0, prev )
713 CALL igerv2d( onedcontext, 1, 1, irecv, 1, 0, prev )
714 ELSE
715 drecv = tmp1
716 irecv = itmp1
717 END IF
718 work( 1 ) = max( lsave, drecv )
719 iwork( 1 ) = irecv
720 alpha = max( alpha, work( 1 ) )
721 nalpha = max( nalpha, irecv )
722 IF( beta-alpha.LE.max( atoli, reltol*max( abs( alpha ),
723 $ abs( beta ) ) ) ) THEN
724 mid = half*( alpha+beta )
725 DO 60 j = offset + nalpha + 1, offset + nbeta
726 work( indrw2+im+1 ) = mid
727 iwork( indriw1+im+1 ) = blkno
728 iwork( indriw2+im+1 ) = j
729 im = im + 1
730 60 CONTINUE
731 work( 2 ) = alpha
732 iwork( 2 ) = nalpha
733 END IF
734 END IF
735 neigint = iwork( 2 ) - iwork( 1 )
736 IF( neigint.LE.0 )
737 $ GO TO 90
738*
739* Call the main computational routine
740*
741 CALL pdlaebz( 2, in, neigint, 1, atoli, reltol, pivmin,
742 $ work( indrw1+2*ioff+1 ), iwork, work, iwork,
743 $ iout, lsave, ieflag, iinfo )
744 IF( iinfo.NE.0 ) THEN
745 info = 1
746 END IF
747 DO 80 i = 1, iout
748 mid = half*( work( 2*i-1 )+work( 2*i ) )
749 IF( i.GT.iout-iinfo )
750 $ blkno = -blkno
751 DO 70 j = offset + iwork( 2*i-1 ) + 1,
752 $ offset + iwork( 2*i )
753 work( indrw2+im+1 ) = mid
754 iwork( indriw1+im+1 ) = blkno
755 iwork( indriw2+im+1 ) = j
756 im = im + 1
757 70 CONTINUE
758 80 CONTINUE
759 90 CONTINUE
760*
761* Find out total number of eigenvalues computed
762*
763 100 CONTINUE
764 m = im
765 CALL igsum2d( onedcontext, 'ALL', ' ', 1, 1, m, 1, -1, -1 )
766*
767* Move the eigenvalues found to their final destinations
768*
769 DO 130 i = 1, p
770 IF( self.EQ.i-1 ) THEN
771 CALL igebs2d( onedcontext, 'ALL', ' ', 1, 1, im, 1 )
772 IF( im.NE.0 ) THEN
773 CALL igebs2d( onedcontext, 'ALL', ' ', im, 1,
774 $ iwork( indriw2+1 ), im )
775 CALL dgebs2d( onedcontext, 'ALL', ' ', im, 1,
776 $ work( indrw2+1 ), im )
777 CALL igebs2d( onedcontext, 'ALL', ' ', im, 1,
778 $ iwork( indriw1+1 ), im )
779 DO 110 j = 1, im
780 w( iwork( indriw2+j ) ) = work( indrw2+j )
781 iblock( iwork( indriw2+j ) ) = iwork( indriw1+j )
782 110 CONTINUE
783 END IF
784 ELSE
785 CALL igebr2d( onedcontext, 'ALL', ' ', 1, 1, torecv, 1, 0,
786 $ i-1 )
787 IF( torecv( 1, 1 ).NE.0 ) THEN
788 CALL igebr2d( onedcontext, 'ALL', ' ', torecv( 1, 1 ), 1,
789 $ iwork, torecv( 1, 1 ), 0, i-1 )
790 CALL dgebr2d( onedcontext, 'ALL', ' ', torecv( 1, 1 ), 1,
791 $ work, torecv( 1, 1 ), 0, i-1 )
792 CALL igebr2d( onedcontext, 'ALL', ' ', torecv( 1, 1 ), 1,
793 $ iwork( n+1 ), torecv( 1, 1 ), 0, i-1 )
794 DO 120 j = 1, torecv( 1, 1 )
795 w( iwork( j ) ) = work( j )
796 iblock( iwork( j ) ) = iwork( n+j )
797 120 CONTINUE
798 END IF
799 END IF
800 130 CONTINUE
801 IF( nsplit.GT.1 .AND. iorder.EQ.1 ) THEN
802*
803* Sort the eigenvalues
804*
805*
806 DO 140 i = 1, m
807 iwork( m+i ) = i
808 140 CONTINUE
809 CALL dlasrt2( 'I', m, w, iwork( m+1 ), iinfo )
810 DO 150 i = 1, m
811 iwork( i ) = iblock( i )
812 150 CONTINUE
813 DO 160 i = 1, m
814 iblock( i ) = iwork( iwork( m+i ) )
815 160 CONTINUE
816 END IF
817 IF( irange.EQ.3 .AND. ( lextra.GT.0 .OR. rextra.GT.0 ) ) THEN
818*
819* Discard unwanted eigenvalues (occurs only when RANGE = 'I',
820* and eigenvalues IL, and/or IU are in a cluster)
821*
822 DO 170 i = 1, m
823 work( i ) = w( i )
824 iwork( i ) = i
825 iwork( m+i ) = i
826 170 CONTINUE
827 DO 190 i = 1, lextra
828 itmp1 = i
829 DO 180 j = i + 1, m
830 IF( work( j ).LT.work( itmp1 ) ) THEN
831 itmp1 = j
832 END IF
833 180 CONTINUE
834 tmp1 = work( i )
835 work( i ) = work( itmp1 )
836 work( itmp1 ) = tmp1
837 iwork( iwork( m+itmp1 ) ) = i
838 iwork( iwork( m+i ) ) = itmp1
839 itmp2 = iwork( m+i )
840 iwork( m+i ) = iwork( m+itmp1 )
841 iwork( m+itmp1 ) = itmp2
842 190 CONTINUE
843 DO 210 i = 1, rextra
844 itmp1 = m - i + 1
845 DO 200 j = m - i, lextra + 1, -1
846 IF( work( j ).GT.work( itmp1 ) ) THEN
847 itmp1 = j
848 END IF
849 200 CONTINUE
850 tmp1 = work( m-i+1 )
851 work( m-i+1 ) = work( itmp1 )
852 work( itmp1 ) = tmp1
853 iwork( iwork( m+itmp1 ) ) = m - i + 1
854 iwork( iwork( 2*m-i+1 ) ) = itmp1
855 itmp2 = iwork( 2*m-i+1 )
856 iwork( 2*m-i+1 ) = iwork( m+itmp1 )
857 iwork( m+itmp1 ) = itmp2
858* IWORK( ITMP1 ) = 1
859 210 CONTINUE
860 j = 0
861 DO 220 i = 1, m
862 IF( iwork( i ).GT.lextra .AND. iwork( i ).LE.m-rextra ) THEN
863 j = j + 1
864 w( j ) = work( iwork( i ) )
865 iblock( j ) = iblock( i )
866 END IF
867 220 CONTINUE
868 m = m - lextra - rextra
869 END IF
870 IF( m.NE.ilast-ifrst+1 ) THEN
871 info = 2
872 END IF
873*
874 230 CONTINUE
875 CALL blacs_freebuff( onedcontext, 1 )
876 CALL blacs_gridexit( onedcontext )
877 RETURN
878*
879* End of PDSTEBZ
880*
dlasrt2
subroutine dlasrt2(id, n, d, key, info)
Definition dlasrt2.f:4
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
pdlamch
double precision function pdlamch(ictxt, cmach)
Definition pdblastst.f:6769
pdlaebz
subroutine pdlaebz(ijob, n, mmax, minp, abstol, reltol, pivmin, d, nval, intvl, intvlct, mout, lsave, ieflag, info)
Definition pdstebz.f:886
pdlapdct
subroutine pdlapdct(sigma, n, d, pivmin, count)
Definition pdstebz.f:1365
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2
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: