SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pdstebz.f
Go to the documentation of this file.
1 SUBROUTINE pdstebz( ICTXT, RANGE, ORDER, N, VL, VU, IL, IU,
2 $ ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT,
3 $ WORK, LWORK, IWORK, LIWORK, INFO )
4*
5* -- ScaLAPACK routine (version 1.7) --
6* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7* and University of California, Berkeley.
8* November 15, 1997
9*
10* .. Scalar Arguments ..
11 CHARACTER ORDER, RANGE
12 INTEGER ICTXT, IL, INFO, IU, LIWORK, LWORK, M, N,
13 $ nsplit
14 DOUBLE PRECISION ABSTOL, VL, VU
15* ..
16* .. Array Arguments ..
17 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
18 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
19* ..
20*
21* Purpose
22* =======
23*
24* PDSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix in
25* parallel. The user may ask for all eigenvalues, all eigenvalues in
26* the interval [VL, VU], or the eigenvalues indexed IL through IU. A
27* static partitioning of work is done at the beginning of PDSTEBZ which
28* results in all processes finding an (almost) equal number of
29* eigenvalues.
30*
31* NOTE : It is assumed that the user is on an IEEE machine. If the user
32* is not on an IEEE mchine, set the compile time flag NO_IEEE
33* to 1 (in SLmake.inc). The features of IEEE arithmetic that
34* are needed for the "fast" Sturm Count are : (a) infinity
35* arithmetic (b) the sign bit of a single precision floating
36* point number is assumed be in the 32nd bit position
37* (c) the sign of negative zero.
38*
39* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
40* Matrix", Report CS41, Computer Science Dept., Stanford
41* University, July 21, 1966.
42*
43* Arguments
44* =========
45*
46* ICTXT (global input) INTEGER
47* The BLACS context handle.
48*
49* RANGE (global input) CHARACTER
50* Specifies which eigenvalues are to be found.
51* = 'A': ("All") all eigenvalues will be found.
52* = 'V': ("Value") all eigenvalues in the interval
53* [VL, VU] will be found.
54* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
55* entire matrix) will be found.
56*
57* ORDER (global input) CHARACTER
58* Specifies the order in which the eigenvalues and their block
59* numbers are stored in W and IBLOCK.
60* = 'B': ("By Block") the eigenvalues will be grouped by
61* split-off block (see IBLOCK, ISPLIT) and
62* ordered from smallest to largest within
63* the block.
64* = 'E': ("Entire matrix")
65* the eigenvalues for the entire matrix
66* will be ordered from smallest to largest.
67*
68* N (global input) INTEGER
69* The order of the tridiagonal matrix T. N >= 0.
70*
71* VL (global input) DOUBLE PRECISION
72* If RANGE='V', the lower bound of the interval to be searched
73* for eigenvalues. Eigenvalues less than VL will not be
74* returned. Not referenced if RANGE='A' or 'I'.
75*
76* VU (global input) DOUBLE PRECISION
77* If RANGE='V', the upper bound of the interval to be searched
78* for eigenvalues. Eigenvalues greater than VU will not be
79* returned. VU must be greater than VL. Not referenced if
80* RANGE='A' or 'I'.
81*
82* IL (global input) INTEGER
83* If RANGE='I', the index (from smallest to largest) of the
84* smallest eigenvalue to be returned. IL must be at least 1.
85* Not referenced if RANGE='A' or 'V'.
86*
87* IU (global input) INTEGER
88* If RANGE='I', the index (from smallest to largest) of the
89* largest eigenvalue to be returned. IU must be at least IL
90* and no greater than N. Not referenced if RANGE='A' or 'V'.
91*
92* ABSTOL (global input) DOUBLE PRECISION
93* The absolute tolerance for the eigenvalues. An eigenvalue
94* (or cluster) is considered to be located if it has been
95* determined to lie in an interval whose width is ABSTOL or
96* less. If ABSTOL is less than or equal to zero, then ULP*|T|
97* will be used, where |T| means the 1-norm of T.
98* Eigenvalues will be computed most accurately when ABSTOL is
99* set to the underflow threshold DLAMCH('U'), not zero.
100* Note : If eigenvectors are desired later by inverse iteration
101* ( PDSTEIN ), ABSTOL should be set to 2*PDLAMCH('S').
102*
103* D (global input) DOUBLE PRECISION array, dimension (N)
104* The n diagonal elements of the tridiagonal matrix T. To
105* avoid overflow, the matrix must be scaled so that its largest
106* entry is no greater than overflow**(1/2) * underflow**(1/4)
107* in absolute value, and for greatest accuracy, it should not
108* be much smaller than that.
109*
110* E (global input) DOUBLE PRECISION array, dimension (N-1)
111* The (n-1) off-diagonal elements of the tridiagonal matrix T.
112* To avoid overflow, the matrix must be scaled so that its
113* largest entry is no greater than overflow**(1/2) *
114* underflow**(1/4) in absolute value, and for greatest
115* accuracy, it should not be much smaller than that.
116*
117* M (global output) INTEGER
118* The actual number of eigenvalues found. 0 <= M <= N.
119* (See also the description of INFO=2)
120*
121* NSPLIT (global output) INTEGER
122* The number of diagonal blocks in the matrix T.
123* 1 <= NSPLIT <= N.
124*
125* W (global output) DOUBLE PRECISION array, dimension (N)
126* On exit, the first M elements of W contain the eigenvalues
127* on all processes.
128*
129* IBLOCK (global output) INTEGER array, dimension (N)
130* At each row/column j where E(j) is zero or small, the
131* matrix T is considered to split into a block diagonal
132* matrix. On exit IBLOCK(i) specifies which block (from 1
133* to the number of blocks) the eigenvalue W(i) belongs to.
134* NOTE: in the (theoretically impossible) event that bisection
135* does not converge for some or all eigenvalues, INFO is set
136* to 1 and the ones for which it did not are identified by a
137* negative block number.
138*
139* ISPLIT (global output) INTEGER array, dimension (N)
140* The splitting points, at which T breaks up into submatrices.
141* The first submatrix consists of rows/columns 1 to ISPLIT(1),
142* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
143* etc., and the NSPLIT-th consists of rows/columns
144* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
145* (Only the first NSPLIT elements will actually be used, but
146* since the user cannot know a priori what value NSPLIT will
147* have, N words must be reserved for ISPLIT.)
148*
149* WORK (local workspace) DOUBLE PRECISION array,
150* dimension ( MAX( 5*N, 7 ) )
151*
152* LWORK (local input) INTEGER
153* size of array WORK must be >= MAX( 5*N, 7 )
154* If LWORK = -1, then LWORK is global input and a workspace
155* query is assumed; the routine only calculates the minimum
156* and optimal size for all work arrays. Each of these
157* values is returned in the first entry of the corresponding
158* work array, and no error message is issued by PXERBLA.
159*
160* IWORK (local workspace) INTEGER array, dimension ( MAX( 4*N, 14 ) )
161*
162* LIWORK (local input) INTEGER
163* size of array IWORK must be >= MAX( 4*N, 14, NPROCS )
164* If LIWORK = -1, then LIWORK is global input and a workspace
165* query is assumed; the routine only calculates the minimum
166* and optimal size for all work arrays. Each of these
167* values is returned in the first entry of the corresponding
168* work array, and no error message is issued by PXERBLA.
169*
170* INFO (global output) INTEGER
171* = 0 : successful exit
172* < 0 : if INFO = -i, the i-th argument had an illegal value
173* > 0 : some or all of the eigenvalues failed to converge or
174* were not computed:
175* = 1 : Bisection failed to converge for some eigenvalues;
176* these eigenvalues are flagged by a negative block
177* number. The effect is that the eigenvalues may not
178* be as accurate as the absolute and relative
179* tolerances. This is generally caused by arithmetic
180* which is less accurate than PDLAMCH says.
181* = 2 : There is a mismatch between the number of
182* eigenvalues output and the number desired.
183* = 3 : RANGE='i', and the Gershgorin interval initially
184* used was incorrect. No eigenvalues were computed.
185* Probable cause: your machine has sloppy floating
186* point arithmetic.
187* Cure: Increase the PARAMETER "FUDGE", recompile,
188* and try again.
189*
190* Internal Parameters
191* ===================
192*
193* RELFAC DOUBLE PRECISION, default = 2.0
194* The relative tolerance. An interval [a,b] lies within
195* "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
196* where "ulp" is the machine precision (distance from 1 to
197* the next larger floating point number.)
198*
199* FUDGE DOUBLE PRECISION, default = 2.0
200* A "fudge factor" to widen the Gershgorin intervals. Ideally,
201* a value of 1 should work, but on machines with sloppy
202* arithmetic, this needs to be larger. The default for
203* publicly released versions should be large enough to handle
204* the worst machine around. Note that this has no effect
205* on the accuracy of the solution.
206*
207* =====================================================================
208*
209* .. Intrinsic Functions ..
210 INTRINSIC abs, dble, ichar, max, min, mod
211* ..
212* .. External Functions ..
213 LOGICAL LSAME
214 INTEGER BLACS_PNUM
215 DOUBLE PRECISION PDLAMCH
216 EXTERNAL lsame, blacs_pnum, pdlamch
217* ..
218* .. External Subroutines ..
219 EXTERNAL blacs_freebuff, blacs_get, blacs_gridexit,
220 $ blacs_gridinfo, blacs_gridmap, dgebr2d,
221 $ dgebs2d, dgerv2d, dgesd2d, dlasrt2, globchk,
222 $ igebr2d, igebs2d, igerv2d, igesd2d, igsum2d,
223 $ pdlaebz, pdlaiectb, pdlaiectl, pdlapdct,
224 $ pdlasnbt, pxerbla
225* ..
226* .. Parameters ..
227 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
228 $ MB_, NB_, RSRC_, CSRC_, LLD_
229 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
230 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
231 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
232 INTEGER BIGNUM, DESCMULT
233 PARAMETER ( BIGNUM = 10000, descmult = 100 )
234 DOUBLE PRECISION ZERO, ONE, TWO, FIVE, HALF
235 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
236 $ five = 5.0d+0, half = 1.0d+0 / two )
237 DOUBLE PRECISION FUDGE, RELFAC
238 PARAMETER ( FUDGE = 2.0d+0, relfac = 2.0d+0 )
239* ..
240* .. Local Scalars ..
241 LOGICAL LQUERY
242 INTEGER BLKNO, FOUND, I, IBEGIN, IEFLAG, IEND, IFRST,
243 $ iinfo, ilast, iload, im, imyload, in, indriw1,
244 $ indriw2, indrw1, indrw2, inxtload, ioff,
245 $ iorder, iout, irange, irecv, irem, itmp1,
246 $ itmp2, j, jb, k, last, lextra, lreq, mycol,
247 $ myrow, nalpha, nbeta, ncmp, neigint, next, ngl,
248 $ nglob, ngu, nint, npcol, nprow, offset,
249 $ onedcontext, p, prev, rextra, rreq, self
250 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*
881 END
882*
883 SUBROUTINE pdlaebz( IJOB, N, MMAX, MINP, ABSTOL, RELTOL, PIVMIN,
884 $ D, NVAL, INTVL, INTVLCT, MOUT, LSAVE, IEFLAG,
885 $ INFO )
886*
887* -- ScaLAPACK routine (version 1.7) --
888* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
889* and University of California, Berkeley.
890* November 15, 1997
891*
892*
893* .. Scalar Arguments ..
894 INTEGER IEFLAG, IJOB, INFO, MINP, MMAX, MOUT, N
895 DOUBLE PRECISION ABSTOL, LSAVE, PIVMIN, RELTOL
896* ..
897* .. Array Arguments ..
898 INTEGER INTVLCT( * ), NVAL( * )
899 DOUBLE PRECISION D( * ), INTVL( * )
900* ..
901*
902* Purpose
903* =======
904*
905* PDLAEBZ contains the iteration loop which computes the eigenvalues
906* contained in the input intervals [ INTVL(2*j-1), INTVL(2*j) ] where
907* j = 1,...,MINP. It uses and computes the function N(w), which is
908* the count of eigenvalues of a symmetric tridiagonal matrix less than
909* or equal to its argument w.
910*
911* This is a ScaLAPACK internal subroutine and arguments are not
912* checked for unreasonable values.
913*
914* Arguments
915* =========
916*
917* IJOB (input) INTEGER
918* Specifies the computation done by PDLAEBZ
919* = 0 : Find an interval with desired values of N(w) at the
920* endpoints of the interval.
921* = 1 : Find a floating point number contained in the initial
922* interval with a desired value of N(w).
923* = 2 : Perform bisection iteration to find eigenvalues of T.
924*
925* N (input) INTEGER
926* The order of the tridiagonal matrix T. N >= 1.
927*
928* MMAX (input) INTEGER
929* The maximum number of intervals that may be generated. If
930* more than MMAX intervals are generated, then PDLAEBZ will
931* quit with INFO = MMAX+1.
932*
933* MINP (input) INTEGER
934* The initial number of intervals. MINP <= MMAX.
935*
936* ABSTOL (input) DOUBLE PRECISION
937* The minimum (absolute) width of an interval. When an interval
938* is narrower than ABSTOL, or than RELTOL times the larger (in
939* magnitude) endpoint, then it is considered to be sufficiently
940* small, i.e., converged.
941* This must be at least zero.
942*
943* RELTOL (input) DOUBLE PRECISION
944* The minimum relative width of an interval. When an interval
945* is narrower than ABSTOL, or than RELTOL times the larger (in
946* magnitude) endpoint, then it is considered to be sufficiently
947* small, i.e., converged.
948* Note : This should be at least radix*machine epsilon.
949*
950* PIVMIN (input) DOUBLE PRECISION
951* The minimum absolute of a "pivot" in the "paranoid"
952* implementation of the Sturm sequence loop. This must be at
953* least max_j |e(j)^2| *safe_min, and at least safe_min, where
954* safe_min is at least the smallest number that can divide 1.0
955* without overflow.
956* See PDLAPDCT for the "paranoid" implementation of the Sturm
957* sequence loop.
958*
959* D (input) DOUBLE PRECISION array, dimension (2*N - 1)
960* Contains the diagonals and the squares of the off-diagonal
961* elements of the tridiagonal matrix T. These elements are
962* assumed to be interleaved in memory for better cache
963* performance. The diagonal entries of T are in the entries
964* D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
965* entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
966* matrix must be scaled so that its largest entry is no greater
967* than overflow**(1/2) * underflow**(1/4) in absolute value,
968* and for greatest accuracy, it should not be much smaller
969* than that.
970*
971* NVAL (input/output) INTEGER array, dimension (4)
972* If IJOB = 0, the desired values of N(w) are in NVAL(1) and
973* NVAL(2).
974* If IJOB = 1, NVAL(2) is the desired value of N(w).
975* If IJOB = 2, not referenced.
976* This array will, in general, be reordered on output.
977*
978* INTVL (input/output) DOUBLE PRECISION array, dimension (2*MMAX)
979* The endpoints of the intervals. INTVL(2*j-1) is the left
980* endpoint of the j-th interval, and INTVL(2*j) is the right
981* endpoint of the j-th interval. The input intervals will,
982* in general, be modified, split and reordered by the
983* calculation.
984* On input, INTVL contains the MINP input intervals.
985* On output, INTVL contains the converged intervals.
986*
987* INTVLCT (input/output) INTEGER array, dimension (2*MMAX)
988* The counts at the endpoints of the intervals. INTVLCT(2*j-1)
989* is the count at the left endpoint of the j-th interval, i.e.,
990* the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
991* count at the right endpoint of the j-th interval.
992* On input, INTVLCT contains the counts at the endpoints of
993* the MINP input intervals.
994* On output, INTVLCT contains the counts at the endpoints of
995* the converged intervals.
996*
997* MOUT (output) INTEGER
998* The number of intervals output.
999*
1000* LSAVE (output) DOUBLE PRECISION
1001* If IJOB = 0 or 2, not referenced.
1002* If IJOB = 1, this is the largest floating point number
1003* encountered which has count N(w) = NVAL(1).
1004*
1005* IEFLAG (input) INTEGER
1006* A flag which indicates whether N(w) should be speeded up by
1007* exploiting IEEE Arithmetic.
1008*
1009* INFO (output) INTEGER
1010* = 0 : All intervals converged.
1011* = 1 - MMAX : The last INFO intervals did not converge.
1012* = MMAX + 1 : More than MMAX intervals were generated.
1013*
1014* =====================================================================
1015*
1016* .. Intrinsic Functions ..
1017 INTRINSIC abs, int, log, max, min
1018* ..
1019* .. External Subroutines ..
1020 EXTERNAL pdlaecv, pdlaiectb, pdlaiectl, pdlapdct
1021* ..
1022* .. Parameters ..
1023 DOUBLE PRECISION ZERO, TWO, HALF
1024 parameter( zero = 0.0d+0, two = 2.0d+0,
1025 $ half = 1.0d+0 / two )
1026* ..
1027* .. Local Scalars ..
1028 INTEGER I, ITMAX, J, K, KF, KL, KLNEW, L, LCNT, LREQ,
1029 $ NALPHA, NBETA, NMID, RCNT, RREQ
1030 DOUBLE PRECISION ALPHA, BETA, MID
1031* ..
1032* .. Executable Statements ..
1033*
1034 kf = 1
1035 kl = minp + 1
1036 info = 0
1037 IF( intvl( 2 )-intvl( 1 ).LE.zero ) THEN
1038 info = minp
1039 mout = kf
1040 RETURN
1041 END IF
1042 IF( ijob.EQ.0 ) THEN
1043*
1044* Check if some input intervals have "converged"
1045*
1046 CALL pdlaecv( 0, kf, kl, intvl, intvlct, nval,
1047 $ max( abstol, pivmin ), reltol )
1048 IF( kf.GE.kl )
1049 $ GO TO 60
1050*
1051* Compute upper bound on number of iterations needed
1052*
1053 itmax = int( ( log( intvl( 2 )-intvl( 1 )+pivmin )-
1054 $ log( pivmin ) ) / log( two ) ) + 2
1055*
1056* Iteration Loop
1057*
1058 DO 20 i = 1, itmax
1059 klnew = kl
1060 DO 10 j = kf, kl - 1
1061 k = 2*j
1062*
1063* Bisect the interval and find the count at that point
1064*
1065 mid = half*( intvl( k-1 )+intvl( k ) )
1066 IF( ieflag.EQ.0 ) THEN
1067 CALL pdlapdct( mid, n, d, pivmin, nmid )
1068 ELSE IF( ieflag.EQ.1 ) THEN
1069 CALL pdlaiectb( mid, n, d, nmid )
1070 ELSE
1071 CALL pdlaiectl( mid, n, d, nmid )
1072 END IF
1073 lreq = nval( k-1 )
1074 rreq = nval( k )
1075 IF( kl.EQ.1 )
1076 $ nmid = min( intvlct( k ),
1077 $ max( intvlct( k-1 ), nmid ) )
1078 IF( nmid.LE.nval( k-1 ) ) THEN
1079 intvl( k-1 ) = mid
1080 intvlct( k-1 ) = nmid
1081 END IF
1082 IF( nmid.GE.nval( k ) ) THEN
1083 intvl( k ) = mid
1084 intvlct( k ) = nmid
1085 END IF
1086 IF( nmid.GT.lreq .AND. nmid.LT.rreq ) THEN
1087 l = 2*klnew
1088 intvl( l-1 ) = mid
1089 intvl( l ) = intvl( k )
1090 intvlct( l-1 ) = nval( k )
1091 intvlct( l ) = intvlct( k )
1092 intvl( k ) = mid
1093 intvlct( k ) = nval( k-1 )
1094 nval( l-1 ) = nval( k )
1095 nval( l ) = nval( l-1 )
1096 nval( k ) = nval( k-1 )
1097 klnew = klnew + 1
1098 END IF
1099 10 CONTINUE
1100 kl = klnew
1101 CALL pdlaecv( 0, kf, kl, intvl, intvlct, nval,
1102 $ max( abstol, pivmin ), reltol )
1103 IF( kf.GE.kl )
1104 $ GO TO 60
1105 20 CONTINUE
1106 ELSE IF( ijob.EQ.1 ) THEN
1107 alpha = intvl( 1 )
1108 beta = intvl( 2 )
1109 nalpha = intvlct( 1 )
1110 nbeta = intvlct( 2 )
1111 lsave = alpha
1112 lreq = nval( 1 )
1113 rreq = nval( 2 )
1114 30 CONTINUE
1115 IF( nbeta.NE.rreq .AND. beta-alpha.GT.
1116 $ max( abstol, reltol*max( abs( alpha ), abs( beta ) ) ) )
1117 $ THEN
1118*
1119* Bisect the interval and find the count at that point
1120*
1121 mid = half*( alpha+beta )
1122 IF( ieflag.EQ.0 ) THEN
1123 CALL pdlapdct( mid, n, d, pivmin, nmid )
1124 ELSE IF( ieflag.EQ.1 ) THEN
1125 CALL pdlaiectb( mid, n, d, nmid )
1126 ELSE
1127 CALL pdlaiectl( mid, n, d, nmid )
1128 END IF
1129 nmid = min( nbeta, max( nalpha, nmid ) )
1130 IF( nmid.GE.rreq ) THEN
1131 beta = mid
1132 nbeta = nmid
1133 ELSE
1134 alpha = mid
1135 nalpha = nmid
1136 IF( nmid.EQ.lreq )
1137 $ lsave = alpha
1138 END IF
1139 GO TO 30
1140 END IF
1141 kl = kf
1142 intvl( 1 ) = alpha
1143 intvl( 2 ) = beta
1144 intvlct( 1 ) = nalpha
1145 intvlct( 2 ) = nbeta
1146 ELSE IF( ijob.EQ.2 ) THEN
1147*
1148* Check if some input intervals have "converged"
1149*
1150 CALL pdlaecv( 1, kf, kl, intvl, intvlct, nval,
1151 $ max( abstol, pivmin ), reltol )
1152 IF( kf.GE.kl )
1153 $ GO TO 60
1154*
1155* Compute upper bound on number of iterations needed
1156*
1157 itmax = int( ( log( intvl( 2 )-intvl( 1 )+pivmin )-
1158 $ log( pivmin ) ) / log( two ) ) + 2
1159*
1160* Iteration Loop
1161*
1162 DO 50 i = 1, itmax
1163 klnew = kl
1164 DO 40 j = kf, kl - 1
1165 k = 2*j
1166 mid = half*( intvl( k-1 )+intvl( k ) )
1167 IF( ieflag.EQ.0 ) THEN
1168 CALL pdlapdct( mid, n, d, pivmin, nmid )
1169 ELSE IF( ieflag.EQ.1 ) THEN
1170 CALL pdlaiectb( mid, n, d, nmid )
1171 ELSE
1172 CALL pdlaiectl( mid, n, d, nmid )
1173 END IF
1174 lcnt = intvlct( k-1 )
1175 rcnt = intvlct( k )
1176 nmid = min( rcnt, max( lcnt, nmid ) )
1177*
1178* Form New Interval(s)
1179*
1180 IF( nmid.EQ.lcnt ) THEN
1181 intvl( k-1 ) = mid
1182 ELSE IF( nmid.EQ.rcnt ) THEN
1183 intvl( k ) = mid
1184 ELSE IF( klnew.LT.mmax+1 ) THEN
1185 l = 2*klnew
1186 intvl( l-1 ) = mid
1187 intvl( l ) = intvl( k )
1188 intvlct( l-1 ) = nmid
1189 intvlct( l ) = intvlct( k )
1190 intvl( k ) = mid
1191 intvlct( k ) = nmid
1192 klnew = klnew + 1
1193 ELSE
1194 info = mmax + 1
1195 RETURN
1196 END IF
1197 40 CONTINUE
1198 kl = klnew
1199 CALL pdlaecv( 1, kf, kl, intvl, intvlct, nval,
1200 $ max( abstol, pivmin ), reltol )
1201 IF( kf.GE.kl )
1202 $ GO TO 60
1203 50 CONTINUE
1204 END IF
1205 60 CONTINUE
1206 info = max( kl-kf, 0 )
1207 mout = kl - 1
1208 RETURN
1209*
1210* End of PDLAEBZ
1211*
1212 END
1213*
1214*
1215 SUBROUTINE pdlaecv( IJOB, KF, KL, INTVL, INTVLCT, NVAL, ABSTOL,
1216 $ RELTOL )
1217*
1218* -- ScaLAPACK routine (version 1.7) --
1219* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
1220* and University of California, Berkeley.
1221* November 15, 1997
1222*
1223*
1224* .. Scalar Arguments ..
1225 INTEGER IJOB, KF, KL
1226 DOUBLE PRECISION ABSTOL, RELTOL
1227* ..
1228* .. Array Arguments ..
1229 INTEGER INTVLCT( * ), NVAL( * )
1230 DOUBLE PRECISION INTVL( * )
1231* ..
1232*
1233* Purpose
1234* =======
1235*
1236* PDLAECV checks if the input intervals [ INTVL(2*i-1), INTVL(2*i) ],
1237* i = KF, ... , KL-1, have "converged".
1238* PDLAECV modifies KF to be the index of the last converged interval,
1239* i.e., on output, all intervals [ INTVL(2*i-1), INTVL(2*i) ], i < KF,
1240* have converged. Note that the input intervals may be reordered by
1241* PDLAECV.
1242*
1243* This is a SCALAPACK internal procedure and arguments are not checked
1244* for unreasonable values.
1245*
1246* Arguments
1247* =========
1248*
1249* IJOB (input) INTEGER
1250* Specifies the criterion for "convergence" of an interval.
1251* = 0 : When an interval is narrower than ABSTOL, or than
1252* RELTOL times the larger (in magnitude) endpoint, then
1253* it is considered to have "converged".
1254* = 1 : When an interval is narrower than ABSTOL, or than
1255* RELTOL times the larger (in magnitude) endpoint, or if
1256* the counts at the endpoints are identical to the counts
1257* specified by NVAL ( see NVAL ) then the interval is
1258* considered to have "converged".
1259*
1260* KF (input/output) INTEGER
1261* On input, the index of the first input interval is 2*KF-1.
1262* On output, the index of the last converged interval
1263* is 2*KF-3.
1264*
1265* KL (input) INTEGER
1266* The index of the last input interval is 2*KL-3.
1267*
1268* INTVL (input/output) DOUBLE PRECISION array, dimension (2*(KL-KF))
1269* The endpoints of the intervals. INTVL(2*j-1) is the left
1270* oendpoint f the j-th interval, and INTVL(2*j) is the right
1271* endpoint of the j-th interval. The input intervals will,
1272* in general, be reordered on output.
1273* On input, INTVL contains the KL-KF input intervals.
1274* On output, INTVL contains the converged intervals, 1 thru'
1275* KF-1, and the unconverged intervals, KF thru' KL-1.
1276*
1277* INTVLCT (input/output) INTEGER array, dimension (2*(KL-KF))
1278* The counts at the endpoints of the intervals. INTVLCT(2*j-1)
1279* is the count at the left endpoint of the j-th interval, i.e.,
1280* the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
1281* count at the right endpoint of the j-th interval. This array
1282* will, in general, be reordered on output.
1283* See the comments in PDLAEBZ for more on the function N(w).
1284*
1285* NVAL (input/output) INTEGER array, dimension (2*(KL-KF))
1286* The desired counts, N(w), at the endpoints of the
1287* corresponding intervals. This array will, in general,
1288* be reordered on output.
1289*
1290* ABSTOL (input) DOUBLE PRECISION
1291* The minimum (absolute) width of an interval. When an interval
1292* is narrower than ABSTOL, or than RELTOL times the larger (in
1293* magnitude) endpoint, then it is considered to be sufficiently
1294* small, i.e., converged.
1295* Note : This must be at least zero.
1296*
1297* RELTOL (input) DOUBLE PRECISION
1298* The minimum relative width of an interval. When an interval
1299* is narrower than ABSTOL, or than RELTOL times the larger (in
1300* magnitude) endpoint, then it is considered to be sufficiently
1301* small, i.e., converged.
1302* Note : This should be at least radix*machine epsilon.
1303*
1304* =====================================================================
1305*
1306* .. Intrinsic Functions ..
1307 INTRINSIC abs, max
1308* ..
1309* .. Local Scalars ..
1310 LOGICAL CONDN
1311 INTEGER I, ITMP1, ITMP2, J, K, KFNEW
1312 DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4
1313* ..
1314* .. Executable Statements ..
1315*
1316 kfnew = kf
1317 DO 10 i = kf, kl - 1
1318 k = 2*i
1319 tmp3 = intvl( k-1 )
1320 tmp4 = intvl( k )
1321 tmp1 = abs( tmp4-tmp3 )
1322 tmp2 = max( abs( tmp3 ), abs( tmp4 ) )
1323 condn = tmp1.LT.max( abstol, reltol*tmp2 )
1324 IF( ijob.EQ.0 )
1325 $ condn = condn .OR. ( ( intvlct( k-1 ).EQ.nval( k-1 ) ) .AND.
1326 $ intvlct( k ).EQ.nval( k ) )
1327 IF( condn ) THEN
1328 IF( i.GT.kfnew ) THEN
1329*
1330* Reorder Intervals
1331*
1332 j = 2*kfnew
1333 tmp1 = intvl( k-1 )
1334 tmp2 = intvl( k )
1335 itmp1 = intvlct( k-1 )
1336 itmp2 = intvlct( k )
1337 intvl( k-1 ) = intvl( j-1 )
1338 intvl( k ) = intvl( j )
1339 intvlct( k-1 ) = intvlct( j-1 )
1340 intvlct( k ) = intvlct( j )
1341 intvl( j-1 ) = tmp1
1342 intvl( j ) = tmp2
1343 intvlct( j-1 ) = itmp1
1344 intvlct( j ) = itmp2
1345 IF( ijob.EQ.0 ) THEN
1346 itmp1 = nval( k-1 )
1347 nval( k-1 ) = nval( j-1 )
1348 nval( j-1 ) = itmp1
1349 itmp1 = nval( k )
1350 nval( k ) = nval( j )
1351 nval( j ) = itmp1
1352 END IF
1353 END IF
1354 kfnew = kfnew + 1
1355 END IF
1356 10 CONTINUE
1357 kf = kfnew
1358 RETURN
1359*
1360* End of PDLAECV
1361*
1362 END
1363*
1364 SUBROUTINE pdlapdct( SIGMA, N, D, PIVMIN, COUNT )
1365*
1366* -- ScaLAPACK routine (version 1.7) --
1367* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
1368* and University of California, Berkeley.
1369* November 15, 1997
1370*
1371*
1372* .. Scalar Arguments ..
1373 INTEGER COUNT, N
1374 DOUBLE PRECISION PIVMIN, SIGMA
1375* ..
1376* .. Array Arguments ..
1377 DOUBLE PRECISION D( * )
1378* ..
1379*
1380* Purpose
1381* =======
1382*
1383* PDLAPDCT counts the number of negative eigenvalues of (T - SIGMA I).
1384* This implementation of the Sturm Sequence loop has conditionals in
1385* the innermost loop to avoid overflow and determine the sign of a
1386* floating point number. PDLAPDCT will be referred to as the "paranoid"
1387* implementation of the Sturm Sequence loop.
1388*
1389* This is a SCALAPACK internal procedure and arguments are not checked
1390* for unreasonable values.
1391*
1392* Arguments
1393* =========
1394*
1395* SIGMA (input) DOUBLE PRECISION
1396* The shift. PDLAPDCT finds the number of eigenvalues of T less
1397* than or equal to SIGMA.
1398*
1399* N (input) INTEGER
1400* The order of the tridiagonal matrix T. N >= 1.
1401*
1402* D (input) DOUBLE PRECISION array, dimension (2*N - 1)
1403* Contains the diagonals and the squares of the off-diagonal
1404* elements of the tridiagonal matrix T. These elements are
1405* assumed to be interleaved in memory for better cache
1406* performance. The diagonal entries of T are in the entries
1407* D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
1408* entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
1409* matrix must be scaled so that its largest entry is no greater
1410* than overflow**(1/2) * underflow**(1/4) in absolute value,
1411* and for greatest accuracy, it should not be much smaller
1412* than that.
1413*
1414* PIVMIN (input) DOUBLE PRECISION
1415* The minimum absolute of a "pivot" in this "paranoid"
1416* implementation of the Sturm sequence loop. This must be at
1417* least max_j |e(j)^2| *safe_min, and at least safe_min, where
1418* safe_min is at least the smallest number that can divide 1.0
1419* without overflow.
1420*
1421* COUNT (output) INTEGER
1422* The count of the number of eigenvalues of T less than or
1423* equal to SIGMA.
1424*
1425* =====================================================================
1426*
1427* .. Intrinsic Functions ..
1428 INTRINSIC abs
1429* ..
1430* .. Parameters ..
1431 DOUBLE PRECISION ZERO
1432 PARAMETER ( ZERO = 0.0d+0 )
1433* ..
1434* .. Local Scalars ..
1435 INTEGER I
1436 DOUBLE PRECISION TMP
1437* ..
1438* .. Executable Statements ..
1439*
1440 tmp = d( 1 ) - sigma
1441 IF( abs( tmp ).LE.pivmin )
1442 $ tmp = -pivmin
1443 count = 0
1444 IF( tmp.LE.zero )
1445 $ count = 1
1446 DO 10 i = 3, 2*n - 1, 2
1447 tmp = d( i ) - d( i-1 ) / tmp - sigma
1448 IF( abs( tmp ).LE.pivmin )
1449 $ tmp = -pivmin
1450 IF( tmp.LE.zero )
1451 $ count = count + 1
1452 10 CONTINUE
1453*
1454 RETURN
1455*
1456* End of PDLAPDCT
1457*
1458 END
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
pdlaebz
subroutine pdlaebz(ijob, n, mmax, minp, abstol, reltol, pivmin, d, nval, intvl, intvlct, mout, lsave, ieflag, info)
Definition pdstebz.f:886
pdlaecv
subroutine pdlaecv(ijob, kf, kl, intvl, intvlct, nval, abstol, reltol)
Definition pdstebz.f:1217
pdstebz
subroutine pdstebz(ictxt, range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, lwork, iwork, liwork, info)
Definition pdstebz.f:4
pdlapdct
subroutine pdlapdct(sigma, n, d, pivmin, count)
Definition pdstebz.f:1365
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2