ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
dlarrd2.f
Go to the documentation of this file.
1  SUBROUTINE dlarrd2( 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  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 *
678  END
max
#define max(A, B)
Definition: pcgemr.c:180
dlarrd2
subroutine dlarrd2(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: dlarrd2.f:5
min
#define min(A, B)
Definition: pcgemr.c:181