SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slarrd2.f
Go to the documentation of this file.
1 SUBROUTINE slarrd2( RANGE, ORDER, N, VL, VU, IL, IU, GERS,
2 $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
3 $ M, W, WERR, WL, WU, IBLOCK, INDEXW,
4 $ WORK, IWORK, DOL, DOU, INFO )
5*
6* -- ScaLAPACK auxiliary routine (version 2.0) --
7* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
8* July 4, 2010
9*
10* .. Scalar Arguments ..
11 CHARACTER ORDER, RANGE
12 INTEGER DOL, DOU, IL, INFO, IU, M, N, NSPLIT
13 REAL PIVMIN, RELTOL, VL, VU, WL, WU
14* ..
15* .. Array Arguments ..
16 INTEGER IBLOCK( * ), INDEXW( * ),
17 $ ISPLIT( * ), IWORK( * )
18 REAL D( * ), E( * ), E2( * ),
19 $ gers( * ), w( * ), werr( * ), work( * )
20* ..
21*
22* Purpose
23* =======
24*
25* SLARRD2 computes the eigenvalues of a symmetric tridiagonal
26* matrix T to limited initial accuracy. This is an auxiliary code to be
27* called from SLARRE2A.
28*
29* SLARRD2 has been created using the LAPACK code SLARRD
30* which itself stems from SSTEBZ. The motivation for creating
31* SLARRD2 is efficiency: When computing eigenvalues in parallel
32* and the input tridiagonal matrix splits into blocks, SLARRD2
33* can skip over blocks which contain none of the eigenvalues from
34* DOL to DOU for which the processor responsible. In extreme cases (such
35* as large matrices consisting of many blocks of small size, e.g. 2x2,
36* the gain can be substantial.
37*
38* Arguments
39* =========
40*
41* RANGE (input) CHARACTER
42* = 'A': ("All") all eigenvalues will be found.
43* = 'V': ("Value") all eigenvalues in the half-open interval
44* (VL, VU] will be found.
45* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
46* entire matrix) will be found.
47*
48* ORDER (input) CHARACTER
49* = 'B': ("By Block") the eigenvalues will be grouped by
50* split-off block (see IBLOCK, ISPLIT) and
51* ordered from smallest to largest within
52* the block.
53* = 'E': ("Entire matrix")
54* the eigenvalues for the entire matrix
55* will be ordered from smallest to
56* largest.
57*
58* N (input) INTEGER
59* The order of the tridiagonal matrix T. N >= 0.
60*
61* VL (input) REAL
62* VU (input) REAL
63* If RANGE='V', the lower and upper bounds of the interval to
64* be searched for eigenvalues. Eigenvalues less than or equal
65* to VL, or greater than VU, will not be returned. VL < VU.
66* Not referenced if RANGE = 'A' or 'I'.
67*
68* IL (input) INTEGER
69* IU (input) INTEGER
70* If RANGE='I', the indices (in ascending order) of the
71* smallest and largest eigenvalues to be returned.
72* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
73* Not referenced if RANGE = 'A' or 'V'.
74*
75* GERS (input) REAL array, dimension (2*N)
76* The N Gerschgorin intervals (the i-th Gerschgorin interval
77* is (GERS(2*i-1), GERS(2*i)).
78*
79* RELTOL (input) REAL
80* The minimum relative width of an interval. When an interval
81* is narrower than RELTOL times the larger (in
82* magnitude) endpoint, then it is considered to be
83* sufficiently small, i.e., converged. Note: this should
84* always be at least radix*machine epsilon.
85*
86* D (input) REAL array, dimension (N)
87* The n diagonal elements of the tridiagonal matrix T.
88*
89* E (input) REAL array, dimension (N-1)
90* The (n-1) off-diagonal elements of the tridiagonal matrix T.
91*
92* E2 (input) REAL array, dimension (N-1)
93* The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
94*
95* PIVMIN (input) REAL
96* The minimum pivot allowed in the sturm sequence for T.
97*
98* NSPLIT (input) INTEGER
99* The number of diagonal blocks in the matrix T.
100* 1 <= NSPLIT <= N.
101*
102* ISPLIT (input) INTEGER array, dimension (N)
103* The splitting points, at which T breaks up into submatrices.
104* The first submatrix consists of rows/columns 1 to ISPLIT(1),
105* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
106* etc., and the NSPLIT-th consists of rows/columns
107* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
108* (Only the first NSPLIT elements will actually be used, but
109* since the user cannot know a priori what value NSPLIT will
110* have, N words must be reserved for ISPLIT.)
111*
112* M (output) INTEGER
113* The actual number of eigenvalues found. 0 <= M <= N.
114* (See also the description of INFO=2,3.)
115*
116* W (output) REAL array, dimension (N)
117* On exit, the first M elements of W will contain the
118* eigenvalue approximations. SLARRD2 computes an interval
119* I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
120* approximation is given as the interval midpoint
121* W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
122* WERR(j) = abs( a_j - b_j)/2
123*
124* WERR (output) REAL array, dimension (N)
125* The error bound on the corresponding eigenvalue approximation
126* in W.
127*
128* WL (output) REAL
129* WU (output) REAL
130* The interval (WL, WU] contains all the wanted eigenvalues.
131* If RANGE='V', then WL=VL and WU=VU.
132* If RANGE='A', then WL and WU are the global Gerschgorin bounds
133* on the spectrum.
134* If RANGE='I', then WL and WU are computed by SLAEBZ from the
135* index range specified.
136*
137* IBLOCK (output) INTEGER array, dimension (N)
138* At each row/column j where E(j) is zero or small, the
139* matrix T is considered to split into a block diagonal
140* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
141* block (from 1 to the number of blocks) the eigenvalue W(i)
142* belongs. (SLARRD2 may use the remaining N-M elements as
143* workspace.)
144*
145* INDEXW (output) INTEGER array, dimension (N)
146* The indices of the eigenvalues within each block (submatrix);
147* for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
148* i-th eigenvalue W(i) is the j-th eigenvalue in block k.
149*
150* WORK (workspace) REAL array, dimension (4*N)
151*
152* IWORK (workspace) INTEGER array, dimension (3*N)
153*
154* DOL (input) INTEGER
155* DOU (input) INTEGER
156* If the user wants to work on only a selected part of the
157* representation tree, he can specify an index range DOL:DOU.
158* Otherwise, the setting DOL=1, DOU=N should be applied.
159* Note that DOL and DOU refer to the order in which the eigenvalues
160* are stored in W.
161*
162* INFO (output) INTEGER
163* = 0: successful exit
164* < 0: if INFO = -i, the i-th argument had an illegal value
165* > 0: some or all of the eigenvalues failed to converge or
166* were not computed:
167* =1 or 3: Bisection failed to converge for some
168* eigenvalues; these eigenvalues are flagged by a
169* negative block number. The effect is that the
170* eigenvalues may not be as accurate as the
171* absolute and relative tolerances. This is
172* generally caused by unexpectedly inaccurate
173* arithmetic.
174* =2 or 3: RANGE='I' only: Not all of the eigenvalues
175* IL:IU were found.
176* Effect: M < IU+1-IL
177* Cause: non-monotonic arithmetic, causing the
178* Sturm sequence to be non-monotonic.
179* Cure: recalculate, using RANGE='A', and pick
180* out eigenvalues IL:IU. In some cases,
181* increasing the PARAMETER "FUDGE" may
182* make things work.
183* = 4: RANGE='I', and the Gershgorin interval
184* initially used was too small. No eigenvalues
185* were computed.
186* Probable cause: your machine has sloppy
187* floating-point arithmetic.
188* Cure: Increase the PARAMETER "FUDGE",
189* recompile, and try again.
190*
191* Internal Parameters
192* ===================
193*
194* FUDGE REAL , default = 2 originally, increased to 10.
195* A "fudge factor" to widen the Gershgorin intervals. Ideally,
196* a value of 1 should work, but on machines with sloppy
197* arithmetic, this needs to be larger. The default for
198* publicly released versions should be large enough to handle
199* the worst machine around. Note that this has no effect
200* on accuracy of the solution.
201*
202* =====================================================================
203*
204* .. Parameters ..
205 REAL ZERO, ONE, TWO, HALF, FUDGE
206 PARAMETER ( ZERO = 0.0e0, one = 1.0e0,
207 $ two = 2.0e0, half = one/two,
208 $ fudge = 10.0e0 )
209
210* ..
211* .. Local Scalars ..
212 LOGICAL NCNVRG, TOOFEW
213 INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
214 $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
215 $ itmp1, itmp2, iw, iwoff, j, jblk, jdisc, je,
216 $ jee, nb, nwl, nwu
217 REAL ATOLI, EPS, GL, GU, RTOLI, SPDIAM, TMP1, TMP2,
218 $ TNORM, UFLOW, WKILL, WLU, WUL
219
220* ..
221* .. Local Arrays ..
222 INTEGER IDUMMA( 1 )
223* ..
224* .. External Functions ..
225 LOGICAL LSAME
226 INTEGER ILAENV
227 REAL SLAMCH
228 EXTERNAL lsame, ilaenv, slamch
229* ..
230* .. External Subroutines ..
231 EXTERNAL slaebz
232* ..
233* .. Intrinsic Functions ..
234 INTRINSIC abs, int, log, max, min, sqrt
235* ..
236* .. Executable Statements ..
237*
238 info = 0
239*
240* Decode RANGE
241*
242 IF( lsame( range, 'A' ) ) THEN
243 irange = 1
244 ELSE IF( lsame( range, 'V' ) ) THEN
245 irange = 2
246 ELSE IF( lsame( range, 'I' ) ) THEN
247 irange = 3
248 ELSE
249 irange = 0
250 END IF
251*
252* Decode ORDER
253*
254 IF( lsame( order, 'B' ) ) THEN
255 iorder = 2
256 ELSE IF( lsame( order, 'E' ) ) THEN
257 iorder = 1
258 ELSE
259 iorder = 0
260 END IF
261*
262* Check for Errors
263*
264 IF( irange.LE.0 ) THEN
265 info = -1
266 ELSE IF( iorder.LE.0 ) THEN
267 info = -2
268 ELSE IF( n.LT.0 ) THEN
269 info = -3
270 ELSE IF( irange.EQ.2 ) THEN
271 IF( vl.GE.vu )
272 $ info = -5
273 ELSE IF( irange.EQ.3 .AND. ( il.LT.1 .OR. il.GT.max( 1, n ) ) )
274 $ THEN
275 info = -6
276 ELSE IF( irange.EQ.3 .AND. ( iu.LT.min( n, il ) .OR. iu.GT.n ) )
277 $ THEN
278 info = -7
279 END IF
280*
281 IF( info.NE.0 ) THEN
282 RETURN
283 END IF
284
285* Initialize error flags
286 info = 0
287 ncnvrg = .false.
288 toofew = .false.
289
290* Quick return if possible
291 m = 0
292 IF( n.EQ.0 ) RETURN
293
294* Simplification:
295 IF( irange.EQ.3 .AND. il.EQ.1 .AND. iu.EQ.n ) irange = 1
296
297* Get machine constants
298 eps = slamch( 'P' )
299 uflow = slamch( 'U' )
300
301
302* Special Case when N=1
303* Treat case of 1x1 matrix for quick return
304 IF( n.EQ.1 ) THEN
305 IF( (irange.EQ.1).OR.
306 $ ((irange.EQ.2).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
307 $ ((irange.EQ.3).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
308 m = 1
309 w(1) = d(1)
310* The computation error of the eigenvalue is zero
311 werr(1) = zero
312 iblock( 1 ) = 1
313 indexw( 1 ) = 1
314 ENDIF
315 RETURN
316 END IF
317
318* NB is the minimum vector length for vector bisection, or 0
319* if only scalar is to be done.
320 nb = ilaenv( 1, 'SSTEBZ', ' ', n, -1, -1, -1 )
321 IF( nb.LE.1 ) nb = 0
322
323* Find global spectral radius
324 gl = d(1)
325 gu = d(1)
326 DO 5 i = 1,n
327 gl = min( gl, gers( 2*i - 1))
328 gu = max( gu, gers(2*i) )
329 5 CONTINUE
330* Compute global Gerschgorin bounds and spectral diameter
331 tnorm = max( abs( gl ), abs( gu ) )
332 gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
333 gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
334 spdiam = gu - gl
335* Input arguments for SLAEBZ:
336* The relative tolerance. An interval (a,b] lies within
337* "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
338 rtoli = reltol
339 atoli = fudge*two*uflow + fudge*two*pivmin
340
341 IF( irange.EQ.3 ) THEN
342
343* RANGE='I': Compute an interval containing eigenvalues
344* IL through IU. The initial interval [GL,GU] from the global
345* Gerschgorin bounds GL and GU is refined by SLAEBZ.
346 itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
347 $ log( two ) ) + 2
348 work( n+1 ) = gl
349 work( n+2 ) = gl
350 work( n+3 ) = gu
351 work( n+4 ) = gu
352 work( n+5 ) = gl
353 work( n+6 ) = gu
354 iwork( 1 ) = -1
355 iwork( 2 ) = -1
356 iwork( 3 ) = n + 1
357 iwork( 4 ) = n + 1
358 iwork( 5 ) = il - 1
359 iwork( 6 ) = iu
360*
361 CALL slaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin,
362 $ d, e, e2, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
363 $ iwork, w, iblock, iinfo )
364 IF( iinfo .NE. 0 ) THEN
365 info = iinfo
366 RETURN
367 END IF
368* On exit, output intervals may not be ordered by ascending negcount
369 IF( iwork( 6 ).EQ.iu ) THEN
370 wl = work( n+1 )
371 wlu = work( n+3 )
372 nwl = iwork( 1 )
373 wu = work( n+4 )
374 wul = work( n+2 )
375 nwu = iwork( 4 )
376 ELSE
377 wl = work( n+2 )
378 wlu = work( n+4 )
379 nwl = iwork( 2 )
380 wu = work( n+3 )
381 wul = work( n+1 )
382 nwu = iwork( 3 )
383 END IF
384* On exit, the interval [WL, WLU] contains a value with negcount NWL,
385* and [WUL, WU] contains a value with negcount NWU.
386 IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
387 info = 4
388 RETURN
389 END IF
390
391 ELSEIF( irange.EQ.2 ) THEN
392 wl = vl
393 wu = vu
394
395 ELSEIF( irange.EQ.1 ) THEN
396 wl = gl
397 wu = gu
398 ENDIF
399
400
401
402* Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
403* NWL accumulates the number of eigenvalues .le. WL,
404* NWU accumulates the number of eigenvalues .le. WU
405 m = 0
406 iend = 0
407 info = 0
408 nwl = 0
409 nwu = 0
410*
411 DO 70 jblk = 1, nsplit
412 ioff = iend
413 ibegin = ioff + 1
414 iend = isplit( jblk )
415 in = iend - ioff
416*
417 IF( irange.EQ.1 ) THEN
418 IF( (iend.LT.dol).OR.(ibegin.GT.dou) ) THEN
419* the local block contains none of eigenvalues that matter
420* to this processor
421 nwu = nwu + in
422 DO 30 j = 1, in
423 m = m + 1
424 iblock( m ) = jblk
425 30 CONTINUE
426 GO TO 70
427 END IF
428 END IF
429
430 IF( in.EQ.1 ) THEN
431* 1x1 block
432 IF( wl.GE.d( ibegin )-pivmin )
433 $ nwl = nwl + 1
434 IF( wu.GE.d( ibegin )-pivmin )
435 $ nwu = nwu + 1
436 IF( irange.EQ.1 .OR. ( wl.LT.d( ibegin )-pivmin .AND. wu.GE.
437 $ d( ibegin )-pivmin ) ) THEN
438 m = m + 1
439 w( m ) = d( ibegin )
440 werr(m) = zero
441* The gap for a single block doesn't matter for the later
442* algorithm and is assigned an arbitrary large value
443 iblock( m ) = jblk
444 indexw( m ) = 1
445 END IF
446 ELSE
447* General Case - block of size IN > 2
448* Compute local Gerschgorin interval and use it as the initial
449* interval for SLAEBZ
450 gu = d( ibegin )
451 gl = d( ibegin )
452 tmp1 = zero
453
454 DO 40 j = ibegin, iend
455 gl = min( gl, gers( 2*j - 1))
456 gu = max( gu, gers(2*j) )
457 40 CONTINUE
458 spdiam = gu - gl
459 gl = gl - fudge*tnorm*eps*in - fudge*pivmin
460 gu = gu + fudge*tnorm*eps*in + fudge*pivmin
461*
462 IF( irange.GT.1 ) THEN
463 IF( gu.LT.wl ) THEN
464* the local block contains none of the wanted eigenvalues
465 nwl = nwl + in
466 nwu = nwu + in
467 GO TO 70
468 END IF
469* refine search interval if possible, only range (WL,WU] matters
470 gl = max( gl, wl )
471 gu = min( gu, wu )
472 IF( gl.GE.gu )
473 $ GO TO 70
474 END IF
475
476* Find negcount of initial interval boundaries GL and GU
477 work( n+1 ) = gl
478 work( n+in+1 ) = gu
479 CALL slaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
480 $ d( ibegin ), e( ibegin ), e2( ibegin ),
481 $ idumma, work( n+1 ), work( n+2*in+1 ), im,
482 $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
483 IF( iinfo .NE. 0 ) THEN
484 info = iinfo
485 RETURN
486 END IF
487*
488 nwl = nwl + iwork( 1 )
489 nwu = nwu + iwork( in+1 )
490 iwoff = m - iwork( 1 )
491
492* Compute Eigenvalues
493 itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /
494 $ log( two ) ) + 2
495 CALL slaebz( 2, itmax, in, in, 1, nb, atoli, rtoli, pivmin,
496 $ d( ibegin ), e( ibegin ), e2( ibegin ),
497 $ idumma, work( n+1 ), work( n+2*in+1 ), iout,
498 $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
499 IF( iinfo .NE. 0 ) THEN
500 info = iinfo
501 RETURN
502 END IF
503*
504* Copy eigenvalues into W and IBLOCK
505* Use -JBLK for block number for unconverged eigenvalues.
506* Loop over the number of output intervals from SLAEBZ
507 DO 60 j = 1, iout
508* eigenvalue approximation is middle point of interval
509 tmp1 = half*( work( j+n )+work( j+in+n ) )
510* semi length of error interval
511 tmp2 = half*abs( work( j+n )-work( j+in+n ) )
512 IF( j.GT.iout-iinfo ) THEN
513* Flag non-convergence.
514 ncnvrg = .true.
515 ib = -jblk
516 ELSE
517 ib = jblk
518 END IF
519 DO 50 je = iwork( j ) + 1 + iwoff,
520 $ iwork( j+in ) + iwoff
521 w( je ) = tmp1
522 werr( je ) = tmp2
523 indexw( je ) = je - iwoff
524 iblock( je ) = ib
525 50 CONTINUE
526 60 CONTINUE
527*
528 m = m + im
529 END IF
530 70 CONTINUE
531
532* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
533* If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
534 IF( irange.EQ.3 ) THEN
535 idiscl = il - 1 - nwl
536 idiscu = nwu - iu
537*
538 IF( idiscl.GT.0 ) THEN
539 im = 0
540 DO 80 je = 1, m
541* Remove some of the smallest eigenvalues from the left so that
542* at the end IDISCL =0. Move all eigenvalues up to the left.
543 IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
544 idiscl = idiscl - 1
545 ELSE
546 im = im + 1
547 w( im ) = w( je )
548 werr( im ) = werr( je )
549 indexw( im ) = indexw( je )
550 iblock( im ) = iblock( je )
551 END IF
552 80 CONTINUE
553 m = im
554 END IF
555 IF( idiscu.GT.0 ) THEN
556* Remove some of the largest eigenvalues from the right so that
557* at the end IDISCU =0. Move all eigenvalues up to the left.
558 im=m+1
559 DO 81 je = m, 1, -1
560 IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
561 idiscu = idiscu - 1
562 ELSE
563 im = im - 1
564 w( im ) = w( je )
565 werr( im ) = werr( je )
566 indexw( im ) = indexw( je )
567 iblock( im ) = iblock( je )
568 END IF
569 81 CONTINUE
570 jee = 0
571 DO 82 je = im, m
572 jee = jee + 1
573 w( jee ) = w( je )
574 werr( jee ) = werr( je )
575 indexw( jee ) = indexw( je )
576 iblock( jee ) = iblock( je )
577 82 CONTINUE
578 m = m-im+1
579 END IF
580
581 IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
582* Code to deal with effects of bad arithmetic. (If N(w) is
583* monotone non-decreasing, this should never happen.)
584* Some low eigenvalues to be discarded are not in (WL,WLU],
585* or high eigenvalues to be discarded are not in (WUL,WU]
586* so just kill off the smallest IDISCL/largest IDISCU
587* eigenvalues, by marking the corresponding IBLOCK = 0
588 IF( idiscl.GT.0 ) THEN
589 wkill = wu
590 DO 100 jdisc = 1, idiscl
591 iw = 0
592 DO 90 je = 1, m
593 IF( iblock( je ).NE.0 .AND.
594 $ ( w( je ).LT.wkill .OR. iw.EQ.0 ) ) THEN
595 iw = je
596 wkill = w( je )
597 END IF
598 90 CONTINUE
599 iblock( iw ) = 0
600 100 CONTINUE
601 END IF
602 IF( idiscu.GT.0 ) THEN
603 wkill = wl
604 DO 120 jdisc = 1, idiscu
605 iw = 0
606 DO 110 je = 1, m
607 IF( iblock( je ).NE.0 .AND.
608 $ ( w( je ).GE.wkill .OR. iw.EQ.0 ) ) THEN
609 iw = je
610 wkill = w( je )
611 END IF
612 110 CONTINUE
613 iblock( iw ) = 0
614 120 CONTINUE
615 END IF
616* Now erase all eigenvalues with IBLOCK set to zero
617 im = 0
618 DO 130 je = 1, m
619 IF( iblock( je ).NE.0 ) THEN
620 im = im + 1
621 w( im ) = w( je )
622 werr( im ) = werr( je )
623 indexw( im ) = indexw( je )
624 iblock( im ) = iblock( je )
625 END IF
626 130 CONTINUE
627 m = im
628 END IF
629 IF( idiscl.LT.0 .OR. idiscu.LT.0 ) THEN
630 toofew = .true.
631 END IF
632 END IF
633*
634 IF(( irange.EQ.1 .AND. m.NE.n ).OR.
635 $ ( irange.EQ.3 .AND. m.NE.iu-il+1 ) ) THEN
636 toofew = .true.
637 END IF
638
639* If ORDER='B',(IBLOCK = 2), do nothing the eigenvalues are already sorted
640* by block.
641* If ORDER='E',(IBLOCK = 1), sort the eigenvalues from smallest to largest
642
643 IF( iorder.EQ.1 .AND. nsplit.GT.1 ) THEN
644 DO 150 je = 1, m - 1
645 ie = 0
646 tmp1 = w( je )
647 DO 140 j = je + 1, m
648 IF( w( j ).LT.tmp1 ) THEN
649 ie = j
650 tmp1 = w( j )
651 END IF
652 140 CONTINUE
653 IF( ie.NE.0 ) THEN
654 tmp2 = werr( ie )
655 itmp1 = iblock( ie )
656 itmp2 = indexw( ie )
657 w( ie ) = w( je )
658 werr( ie ) = werr( je )
659 iblock( ie ) = iblock( je )
660 indexw( ie ) = indexw( je )
661 w( je ) = tmp1
662 werr( je ) = tmp2
663 iblock( je ) = itmp1
664 indexw( je ) = itmp2
665 END IF
666 150 CONTINUE
667 END IF
668*
669 info = 0
670 IF( ncnvrg )
671 $ info = info + 1
672 IF( toofew )
673 $ info = info + 2
674 RETURN
675*
676* End of SLARRD2
677*
678 END
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
slarrd2
subroutine slarrd2(range, order, n, vl, vu, il, iu, gers, reltol, d, e, e2, pivmin, nsplit, isplit, m, w, werr, wl, wu, iblock, indexw, work, iwork, dol, dou, info)
Definition slarrd2.f:5