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

◆ dlarrd2()

subroutine dlarrd2 ( character  range,
character  order,
integer  n,
double precision  vl,
double precision  vu,
integer  il,
integer  iu,
double precision, dimension( * )  gers,
double precision  reltol,
double precision, dimension( * )  d,
double precision, dimension( * )  e,
double precision, dimension( * )  e2,
double precision  pivmin,
integer  nsplit,
integer, dimension( * )  isplit,
integer  m,
double precision, dimension( * )  w,
double precision, dimension( * )  werr,
double precision  wl,
double precision  wu,
integer, dimension( * )  iblock,
integer, dimension( * )  indexw,
double precision, dimension( * )  work,
integer, dimension( * )  iwork,
integer  dol,
integer  dou,
integer  info 
)

Definition at line 1 of file dlarrd2.f.

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 DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
14* ..
15* .. Array Arguments ..
16 INTEGER IBLOCK( * ), INDEXW( * ),
17 $ ISPLIT( * ), IWORK( * )
18 DOUBLE PRECISION D( * ), E( * ), E2( * ),
19 $ GERS( * ), W( * ), WERR( * ), WORK( * )
20* ..
21*
22* Purpose
23* =======
24*
25* DLARRD2 computes the eigenvalues of a symmetric tridiagonal
26* matrix T to limited initial accuracy. This is an auxiliary code to be
27* called from DLARRE2A.
28*
29* DLARRD2 has been created using the LAPACK code DLARRD
30* which itself stems from DSTEBZ. The motivation for creating
31* DLARRD2 is efficiency: When computing eigenvalues in parallel
32* and the input tridiagonal matrix splits into blocks, DLARRD2
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) DOUBLE PRECISION
62* VU (input) DOUBLE PRECISION
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) DOUBLE PRECISION 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) DOUBLE PRECISION
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) DOUBLE PRECISION array, dimension (N)
87* The n diagonal elements of the tridiagonal matrix T.
88*
89* E (input) DOUBLE PRECISION array, dimension (N-1)
90* The (n-1) off-diagonal elements of the tridiagonal matrix T.
91*
92* E2 (input) DOUBLE PRECISION array, dimension (N-1)
93* The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
94*
95* PIVMIN (input) DOUBLE PRECISION
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) DOUBLE PRECISION array, dimension (N)
117* On exit, the first M elements of W will contain the
118* eigenvalue approximations. DLARRD2 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) DOUBLE PRECISION array, dimension (N)
125* The error bound on the corresponding eigenvalue approximation
126* in W.
127*
128* WL (output) DOUBLE PRECISION
129* WU (output) DOUBLE PRECISION
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 DLAEBZ 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. (DLARRD2 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) DOUBLE PRECISION 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 DOUBLE PRECISION, 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 DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE
206 parameter( zero = 0.0d0, one = 1.0d0,
207 $ two = 2.0d0, half = one/two,
208 $ fudge = 10.0d0 )
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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH
228 EXTERNAL lsame, ilaenv, dlamch
229* ..
230* .. External Subroutines ..
231 EXTERNAL dlaebz
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 = dlamch( 'P' )
299 uflow = dlamch( '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, 'DSTEBZ', ' ', 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 DLAEBZ:
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 DLAEBZ.
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 dlaebz( 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 DLAEBZ
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 dlaebz( 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 dlaebz( 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 DLAEBZ
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 DLARRD2
677*
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
dlamch
double precision function dlamch(cmach)
Definition tools.f:10
Here is the caller graph for this function: