ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
dlarre2.f
Go to the documentation of this file.
1  SUBROUTINE dlarre2( RANGE, N, VL, VU, IL, IU, D, E, E2,
2  $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT,
3  $ M, DOL, DOU,
4  $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
5  $ WORK, IWORK, INFO )
6 *
7 * -- ScaLAPACK auxiliary routine (version 2.0) --
8 * Univ. of Tennessee, Univ. of California Berkeley, Univ of Colorado Denver
9 * July 4, 2010
10 *
11 * .. Scalar Arguments ..
12  CHARACTER RANGE
13  INTEGER DOL, DOU, IL, INFO, IU, M, N, NSPLIT
14  DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
15 * ..
16 * .. Array Arguments ..
17  INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
18  $ indexw( * )
19  DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
20  $ W( * ),WERR( * ), WGAP( * ), WORK( * )
21 *
22 * Purpose
23 * =======
24 *
25 * To find the desired eigenvalues of a given real symmetric
26 * tridiagonal matrix T, DLARRE2 sets, via DLARRA,
27 * "small" off-diagonal elements to zero. For each block T_i, it finds
28 * (a) a suitable shift at one end of the block's spectrum,
29 * (b) the root RRR, T_i - sigma_i I = L_i D_i L_i^T, and
30 * (c) eigenvalues of each L_i D_i L_i^T.
31 * The representations and eigenvalues found are then returned to
32 * DSTEGR2 to compute the eigenvectors T.
33 *
34 * DLARRE2 is more suitable for parallel computation than the
35 * original LAPACK code for computing the root RRR and its eigenvalues.
36 * When computing eigenvalues in parallel and the input tridiagonal
37 * matrix splits into blocks, DLARRE2
38 * can skip over blocks which contain none of the eigenvalues from
39 * DOL to DOU for which the processor responsible. In extreme cases (such
40 * as large matrices consisting of many blocks of small size, e.g. 2x2,
41 * the gain can be substantial.
42 *
43 * Arguments
44 * =========
45 *
46 * RANGE (input) CHARACTER
47 * = 'A': ("All") all eigenvalues will be found.
48 * = 'V': ("Value") all eigenvalues in the half-open interval
49 * (VL, VU] will be found.
50 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
51 * entire matrix) will be found.
52 *
53 * N (input) INTEGER
54 * The order of the matrix. N > 0.
55 *
56 * VL (input/output) DOUBLE PRECISION
57 * VU (input/output) DOUBLE PRECISION
58 * If RANGE='V', the lower and upper bounds for the eigenvalues.
59 * Eigenvalues less than or equal to VL, or greater than VU,
60 * will not be returned. VL < VU.
61 * If RANGE='I' or ='A', DLARRE2 computes bounds on the desired
62 * part of the spectrum.
63 *
64 * IL (input) INTEGER
65 * IU (input) INTEGER
66 * If RANGE='I', the indices (in ascending order) of the
67 * smallest and largest eigenvalues to be returned.
68 * 1 <= IL <= IU <= N.
69 *
70 * D (input/output) DOUBLE PRECISION array, dimension (N)
71 * On entry, the N diagonal elements of the tridiagonal
72 * matrix T.
73 * On exit, the N diagonal elements of the diagonal
74 * matrices D_i.
75 *
76 * E (input/output) DOUBLE PRECISION array, dimension (N)
77 * On entry, the first (N-1) entries contain the subdiagonal
78 * elements of the tridiagonal matrix T; E(N) need not be set.
79 * On exit, E contains the subdiagonal elements of the unit
80 * bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
81 * 1 <= I <= NSPLIT, contain the base points sigma_i on output.
82 *
83 * E2 (input/output) DOUBLE PRECISION array, dimension (N)
84 * On entry, the first (N-1) entries contain the SQUARES of the
85 * subdiagonal elements of the tridiagonal matrix T;
86 * E2(N) need not be set.
87 * On exit, the entries E2( ISPLIT( I ) ),
88 * 1 <= I <= NSPLIT, have been set to zero
89 *
90 * RTOL1 (input) DOUBLE PRECISION
91 * RTOL2 (input) DOUBLE PRECISION
92 * Parameters for bisection.
93 * An interval [LEFT,RIGHT] has converged if
94 * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
95 *
96 * SPLTOL (input) DOUBLE PRECISION
97 * The threshold for splitting.
98 *
99 * NSPLIT (output) INTEGER
100 * The number of blocks T splits into. 1 <= NSPLIT <= N.
101 *
102 * ISPLIT (output) INTEGER array, dimension (N)
103 * The splitting points, at which T breaks up into blocks.
104 * The first block 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 *
109 * M (output) INTEGER
110 * The total number of eigenvalues (of all L_i D_i L_i^T)
111 * found.
112 *
113 * DOL (input) INTEGER
114 * DOU (input) INTEGER
115 * If the user wants to work on only a selected part of the
116 * representation tree, he can specify an index range DOL:DOU.
117 * Otherwise, the setting DOL=1, DOU=N should be applied.
118 * Note that DOL and DOU refer to the order in which the eigenvalues
119 * are stored in W.
120 *
121 * W (output) DOUBLE PRECISION array, dimension (N)
122 * The first M elements contain the eigenvalues. The
123 * eigenvalues of each of the blocks, L_i D_i L_i^T, are
124 * sorted in ascending order ( DLARRE2 may use the
125 * remaining N-M elements as workspace).
126 * Note that immediately after exiting this routine, only
127 * the eigenvalues from position DOL:DOU in W might be
128 * reliable on this processor
129 * when the eigenvalue computation is done in parallel.
130 *
131 * WERR (output) DOUBLE PRECISION array, dimension (N)
132 * The error bound on the corresponding eigenvalue in W.
133 * Note that immediately after exiting this routine, only
134 * the uncertainties from position DOL:DOU in WERR might be
135 * reliable on this processor
136 * when the eigenvalue computation is done in parallel.
137 *
138 * WGAP (output) DOUBLE PRECISION array, dimension (N)
139 * The separation from the right neighbor eigenvalue in W.
140 * The gap is only with respect to the eigenvalues of the same block
141 * as each block has its own representation tree.
142 * Exception: at the right end of a block we store the left gap
143 * Note that immediately after exiting this routine, only
144 * the gaps from position DOL:DOU in WGAP might be
145 * reliable on this processor
146 * when the eigenvalue computation is done in parallel.
147 *
148 * IBLOCK (output) INTEGER array, dimension (N)
149 * The indices of the blocks (submatrices) associated with the
150 * corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
151 * W(i) belongs to the first block from the top, =2 if W(i)
152 * belongs to the second block, etc.
153 *
154 * INDEXW (output) INTEGER array, dimension (N)
155 * The indices of the eigenvalues within each block (submatrix);
156 * for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
157 * i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
158 *
159 * GERS (output) DOUBLE PRECISION array, dimension (2*N)
160 * The N Gerschgorin intervals (the i-th Gerschgorin interval
161 * is (GERS(2*i-1), GERS(2*i)).
162 *
163 * PIVMIN (output) DOUBLE PRECISION
164 * The minimum pivot in the sturm sequence for T.
165 *
166 * WORK (workspace) DOUBLE PRECISION array, dimension (6*N)
167 * Workspace.
168 *
169 * IWORK (workspace) INTEGER array, dimension (5*N)
170 * Workspace.
171 *
172 * INFO (output) INTEGER
173 * = 0: successful exit
174 * > 0: A problem occured in DLARRE2.
175 * < 0: One of the called subroutines signaled an internal problem.
176 * Needs inspection of the corresponding parameter IINFO
177 * for further information.
178 *
179 * =-1: Problem in DLARRD.
180 * = 2: No base representation could be found in MAXTRY iterations.
181 * Increasing MAXTRY and recompilation might be a remedy.
182 * =-3: Problem in DLARRB when computing the refined root
183 * representation for DLASQ2.
184 * =-4: Problem in DLARRB when preforming bisection on the
185 * desired part of the spectrum.
186 * =-5: Problem in DLASQ2.
187 * =-6: Problem in DLASQ2.
188 *
189 * =====================================================================
190 *
191 * .. Parameters ..
192  DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
193  $ MAXGROWTH, ONE, PERT, TWO, ZERO
194  PARAMETER ( ZERO = 0.0d0, one = 1.0d0,
195  $ two = 2.0d0, four=4.0d0,
196  $ hndrd = 100.0d0,
197  $ pert = 8.0d0,
198  $ half = one/two, fourth = one/four, fac= half,
199  $ maxgrowth = 64.0d0, fudge = 2.0d0 )
200  INTEGER MAXTRY
201  PARAMETER ( MAXTRY = 6 )
202 * ..
203 * .. Local Scalars ..
204  LOGICAL FORCEB, NOREP, RNDPRT, USEDQD
205  INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
206  $ IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
207  $ WBEGIN, WEND
208  DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
209  $ emax, eold, eps, gl, gu, isleft, isrght, rtl,
210  $ rtol, s1, s2, safmin, sgndef, sigma, spdiam,
211  $ tau, tmp, tmp1
212 
213 
214 * ..
215 * .. Local Arrays ..
216  INTEGER ISEED( 4 )
217 * ..
218 * .. External Functions ..
219  LOGICAL LSAME
220  DOUBLE PRECISION DLAMCH
221  EXTERNAL DLAMCH, LSAME
222 
223 * ..
224 * .. External Subroutines ..
225  EXTERNAL dcopy, dlarnv, dlarra, dlarrb, dlarrc,
226  $ dlarrd, dlasq2
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, max, min
230 
231 * ..
232 * .. Executable Statements ..
233 *
234 
235  info = 0
236 
237 * Dis-/Enable a small random perturbation of the root representation
238  rndprt = .true.
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  END IF
249 
250  m = 0
251 
252 * Get machine constants
253  safmin = dlamch( 'S' )
254  eps = dlamch( 'P' )
255 
256 * Set parameters
257  rtl = sqrt(eps)
258  bsrtol = 1.0d-1
259 
260 * Treat case of 1x1 matrix for quick return
261  IF( n.EQ.1 ) THEN
262  IF( (irange.EQ.1).OR.
263  $ ((irange.EQ.2).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
264  $ ((irange.EQ.3).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
265  m = 1
266  w(1) = d(1)
267 * The computation error of the eigenvalue is zero
268  werr(1) = zero
269  wgap(1) = zero
270  iblock( 1 ) = 1
271  indexw( 1 ) = 1
272  gers(1) = d( 1 )
273  gers(2) = d( 1 )
274  ENDIF
275 * store the shift for the initial RRR, which is zero in this case
276  e(1) = zero
277  RETURN
278  END IF
279 
280 * General case: tridiagonal matrix of order > 1
281 *
282 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
283 * Compute maximum off-diagonal entry and pivmin.
284  gl = d(1)
285  gu = d(1)
286  eold = zero
287  emax = zero
288  e(n) = zero
289  DO 5 i = 1,n
290  werr(i) = zero
291  wgap(i) = zero
292  eabs = abs( e(i) )
293  IF( eabs .GE. emax ) THEN
294  emax = eabs
295  END IF
296  tmp1 = eabs + eold
297  gers( 2*i-1) = d(i) - tmp1
298  gl = min( gl, gers( 2*i - 1))
299  gers( 2*i ) = d(i) + tmp1
300  gu = max( gu, gers(2*i) )
301  eold = eabs
302  5 CONTINUE
303 * The minimum pivot allowed in the sturm sequence for T
304  pivmin = safmin * max( one, emax**2 )
305 * Compute spectral diameter. The Gerschgorin bounds give an
306 * estimate that is wrong by at most a factor of SQRT(2)
307  spdiam = gu - gl
308 
309 * Compute splitting points
310  CALL dlarra( n, d, e, e2, spltol, spdiam,
311  $ nsplit, isplit, iinfo )
312 
313 * Can force use of bisection instead of faster DQDS
314  forceb = .false.
315 
316  IF( (irange.EQ.1) .AND. (.NOT. forceb) ) THEN
317 * Set interval [VL,VU] that contains all eigenvalues
318  vl = gl
319  vu = gu
320  ELSE
321 * We call DLARRD to find crude approximations to the eigenvalues
322 * in the desired range. In case IRANGE = 3, we also obtain the
323 * interval (VL,VU] that contains all the wanted eigenvalues.
324 * An interval [LEFT,RIGHT] has converged if
325 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
326 * DLARRD needs a WORK of size 4*N, IWORK of size 3*N
327  CALL dlarrd( range, 'B', n, vl, vu, il, iu, gers,
328  $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
329  $ mm, w, werr, vl, vu, iblock, indexw,
330  $ work, iwork, iinfo )
331  IF( iinfo.NE.0 ) THEN
332  info = -1
333  RETURN
334  ENDIF
335 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
336  DO 14 i = mm+1,n
337  w( i ) = zero
338  werr( i ) = zero
339  iblock( i ) = 0
340  indexw( i ) = 0
341  14 CONTINUE
342  END IF
343 
344 
345 ***
346 * Loop over unreduced blocks
347  ibegin = 1
348  wbegin = 1
349  DO 170 jblk = 1, nsplit
350  iend = isplit( jblk )
351  in = iend - ibegin + 1
352 
353 * 1 X 1 block
354  IF( in.EQ.1 ) THEN
355  IF( (irange.EQ.1).OR.( (irange.EQ.2).AND.
356  $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
357  $ .OR. ( (irange.EQ.3).AND.(iblock(wbegin).EQ.jblk))
358  $ ) THEN
359  m = m + 1
360  w( m ) = d( ibegin )
361  werr(m) = zero
362 * The gap for a single block doesn't matter for the later
363 * algorithm and is assigned an arbitrary large value
364  wgap(m) = zero
365  iblock( m ) = jblk
366  indexw( m ) = 1
367  wbegin = wbegin + 1
368  ENDIF
369 * E( IEND ) holds the shift for the initial RRR
370  e( iend ) = zero
371  ibegin = iend + 1
372  GO TO 170
373  END IF
374 *
375 * Blocks of size larger than 1x1
376 *
377 * E( IEND ) will hold the shift for the initial RRR, for now set it =0
378  e( iend ) = zero
379 *
380 * Find local outer bounds GL,GU for the block
381  gl = d(ibegin)
382  gu = d(ibegin)
383  DO 15 i = ibegin , iend
384  gl = min( gers( 2*i-1 ), gl )
385  gu = max( gers( 2*i ), gu )
386  15 CONTINUE
387  spdiam = gu - gl
388 
389  IF(.NOT. ((irange.EQ.1).AND.(.NOT.forceb)) ) THEN
390 * Count the number of eigenvalues in the current block.
391  mb = 0
392  DO 20 i = wbegin,mm
393  IF( iblock(i).EQ.jblk ) THEN
394  mb = mb+1
395  ELSE
396  GOTO 21
397  ENDIF
398  20 CONTINUE
399  21 CONTINUE
400 
401  IF( mb.EQ.0) THEN
402 * No eigenvalue in the current block lies in the desired range
403 * E( IEND ) holds the shift for the initial RRR
404  e( iend ) = zero
405  ibegin = iend + 1
406  GO TO 170
407  ELSE
408 
409 * Decide whether dqds or bisection is more efficient
410  usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
411  wend = wbegin + mb - 1
412 * Calculate gaps for the current block
413 * In later stages, when representations for individual
414 * eigenvalues are different, we use SIGMA = E( IEND ).
415  sigma = zero
416  DO 30 i = wbegin, wend - 1
417  wgap( i ) = max( zero,
418  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
419  30 CONTINUE
420  wgap( wend ) = max( zero,
421  $ vu - sigma - (w( wend )+werr( wend )))
422 * Find local index of the first and last desired evalue.
423  indl = indexw(wbegin)
424  indu = indexw( wend )
425  ENDIF
426  ELSE
427 * MB = number of eigenvalues to compute
428  mb = in
429  wend = wbegin + mb - 1
430  indl = 1
431  indu = in
432  ENDIF
433 
434  IF( (wend.LT.dol).OR.(wbegin.GT.dou) ) THEN
435 * if this subblock contains no desired eigenvalues,
436 * skip the computation of this representation tree
437  ibegin = iend + 1
438  wbegin = wend + 1
439  m = m + indu - indl + 1
440  GO TO 170
441  END IF
442 
443 * Find approximations to the extremal eigenvalues of the block
444  CALL dlarrk( in, 1, gl, gu, d(ibegin),
445  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
446  IF( iinfo.NE.0 ) THEN
447  info = -1
448  RETURN
449  ENDIF
450  isleft = max(gl, tmp - tmp1
451  $ - hndrd * eps* abs(tmp - tmp1))
452  CALL dlarrk( in, in, gl, gu, d(ibegin),
453  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
454  IF( iinfo.NE.0 ) THEN
455  info = -1
456  RETURN
457  ENDIF
458  isrght = min(gu, tmp + tmp1
459  $ + hndrd * eps * abs(tmp + tmp1))
460  IF(( (irange.EQ.1) .AND. (.NOT. forceb) ).OR.usedqd) THEN
461 * Case of DQDS
462 * Improve the estimate of the spectral diameter
463  spdiam = isrght - isleft
464  ELSE
465 * Case of bisection
466 * Find approximations to the wanted extremal eigenvalues
467  isleft = max(gl, w(wbegin) - werr(wbegin)
468  $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
469  isrght = min(gu,w(wend) + werr(wend)
470  $ + hndrd * eps * abs(w(wend)+ werr(wend)))
471  ENDIF
472 
473 
474 * Decide whether the base representation for the current block
475 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
476 * should be on the left or the right end of the current block.
477 * The strategy is to shift to the end which is "more populated"
478 * Furthermore, decide whether to use DQDS for the computation of
479 * the eigenvalue approximations at the end of DLARRE2 or bisection.
480 * dqds is chosen if all eigenvalues are desired or the number of
481 * eigenvalues to be computed is large compared to the blocksize.
482  IF( ( irange.EQ.1 ) .AND. (.NOT.forceb) ) THEN
483 * If all the eigenvalues have to be computed, we use dqd
484  usedqd = .true.
485 * INDL is the local index of the first eigenvalue to compute
486  indl = 1
487  indu = in
488 * MB = number of eigenvalues to compute
489  mb = in
490  wend = wbegin + mb - 1
491 * Define 1/4 and 3/4 points of the spectrum
492  s1 = isleft + fourth * spdiam
493  s2 = isrght - fourth * spdiam
494  ELSE
495 * DLARRD has computed IBLOCK and INDEXW for each eigenvalue
496 * approximation.
497 * choose sigma
498  IF( usedqd ) THEN
499  s1 = isleft + fourth * spdiam
500  s2 = isrght - fourth * spdiam
501  ELSE
502  tmp = min(isrght,vu) - max(isleft,vl)
503  s1 = max(isleft,vl) + fourth * tmp
504  s2 = min(isrght,vu) - fourth * tmp
505  ENDIF
506  ENDIF
507 
508 * Compute the negcount at the 1/4 and 3/4 points
509  IF(mb.GT.1) THEN
510  CALL dlarrc( 'T', in, s1, s2, d(ibegin),
511  $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
512  ENDIF
513 
514  IF(mb.EQ.1) THEN
515  sigma = gl
516  sgndef = one
517  ELSEIF( cnt1 - indl .GE. indu - cnt2 ) THEN
518  IF( ( irange.EQ.1 ) .AND. (.NOT.forceb) ) THEN
519  sigma = max(isleft,gl)
520  ELSEIF( usedqd ) THEN
521 * use Gerschgorin bound as shift to get pos def matrix
522 * for dqds
523  sigma = isleft
524  ELSE
525 * use approximation of the first desired eigenvalue of the
526 * block as shift
527  sigma = max(isleft,vl)
528  ENDIF
529  sgndef = one
530  ELSE
531  IF( ( irange.EQ.1 ) .AND. (.NOT.forceb) ) THEN
532  sigma = min(isrght,gu)
533  ELSEIF( usedqd ) THEN
534 * use Gerschgorin bound as shift to get neg def matrix
535 * for dqds
536  sigma = isrght
537  ELSE
538 * use approximation of the first desired eigenvalue of the
539 * block as shift
540  sigma = min(isrght,vu)
541  ENDIF
542  sgndef = -one
543  ENDIF
544 
545 
546 * An initial SIGMA has been chosen that will be used for computing
547 * T - SIGMA I = L D L^T
548 * Define the increment TAU of the shift in case the initial shift
549 * needs to be refined to obtain a factorization with not too much
550 * element growth.
551  IF( usedqd ) THEN
552  tau = spdiam*eps*n + two*pivmin
553  tau = max(tau,eps*abs(sigma))
554  ELSE
555  IF(mb.GT.1) THEN
556  clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
557  avgap = abs(clwdth / dble(wend-wbegin))
558  IF( sgndef.EQ.one ) THEN
559  tau = half*max(wgap(wbegin),avgap)
560  tau = max(tau,werr(wbegin))
561  ELSE
562  tau = half*max(wgap(wend-1),avgap)
563  tau = max(tau,werr(wend))
564  ENDIF
565  ELSE
566  tau = werr(wbegin)
567  ENDIF
568  ENDIF
569 *
570  DO 80 idum = 1, maxtry
571 * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
572 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
573 * pivots in WORK(2*IN+1:3*IN)
574  dpivot = d( ibegin ) - sigma
575  work( 1 ) = dpivot
576  dmax = abs( work(1) )
577  j = ibegin
578  DO 70 i = 1, in - 1
579  work( 2*in+i ) = one / work( i )
580  tmp = e( j )*work( 2*in+i )
581  work( in+i ) = tmp
582  dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
583  work( i+1 ) = dpivot
584  dmax = max( dmax, abs(dpivot) )
585  j = j + 1
586  70 CONTINUE
587 * check for element growth
588  IF( dmax .GT. maxgrowth*spdiam ) THEN
589  norep = .true.
590  ELSE
591  norep = .false.
592  ENDIF
593  IF(norep) THEN
594 * Note that in the case of IRANGE=1, we use the Gerschgorin
595 * shift which makes the matrix definite. So we should end up
596 * here really only in the case of IRANGE = 2,3
597  IF( idum.EQ.maxtry-1 ) THEN
598  IF( sgndef.EQ.one ) THEN
599 * The fudged Gerschgorin shift should succeed
600  sigma =
601  $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
602  ELSE
603  sigma =
604  $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
605  END IF
606  ELSE
607  sigma = sigma - sgndef * tau
608  tau = two * tau
609  END IF
610  ELSE
611 * an initial RRR is found
612  GO TO 83
613  END IF
614  80 CONTINUE
615 * if the program reaches this point, no base representation could be
616 * found in MAXTRY iterations.
617  info = 2
618  RETURN
619 
620  83 CONTINUE
621 * At this point, we have found an initial base representation
622 * T - SIGMA I = L D L^T with not too much element growth.
623 * Store the shift.
624  e( iend ) = sigma
625 * Store D and L.
626  CALL dcopy( in, work, 1, d( ibegin ), 1 )
627  CALL dcopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
628 
629 
630  IF(rndprt .AND. mb.GT.1 ) THEN
631 *
632 * Perturb each entry of the base representation by a small
633 * (but random) relative amount to overcome difficulties with
634 * glued matrices.
635 *
636  DO 122 i = 1, 4
637  iseed( i ) = 1
638  122 CONTINUE
639 
640  CALL dlarnv(2, iseed, 2*in-1, work(1))
641  DO 125 i = 1,in-1
642  d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
643  e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
644  125 CONTINUE
645  d(iend) = d(iend)*(one+eps*four*work(in))
646 *
647  ENDIF
648 *
649 * Don't update the Gerschgorin intervals because keeping track
650 * of the updates would be too much work in DLARRV.
651 * We update W instead and use it to locate the proper Gerschgorin
652 * intervals.
653 
654 * Compute the required eigenvalues of L D L' by bisection or dqds
655  IF ( .NOT.usedqd ) THEN
656 * If DLARRD has been used, shift the eigenvalue approximations
657 * according to their representation. This is necessary for
658 * a uniform DLARRV since dqds computes eigenvalues of the
659 * shifted representation. In DLARRV, W will always hold the
660 * UNshifted eigenvalue approximation.
661  DO 134 j=wbegin,wend
662  w(j) = w(j) - sigma
663  werr(j) = werr(j) + abs(w(j)) * eps
664  134 CONTINUE
665 * call DLARRB to reduce eigenvalue error of the approximations
666 * from DLARRD
667  DO 135 i = ibegin, iend-1
668  work( i ) = d( i ) * e( i )**2
669  135 CONTINUE
670 * use bisection to find EV from INDL to INDU
671  CALL dlarrb(in, d(ibegin), work(ibegin),
672  $ indl, indu, rtol1, rtol2, indl-1,
673  $ w(wbegin), wgap(wbegin), werr(wbegin),
674  $ work( 2*n+1 ), iwork, pivmin, spdiam,
675  $ in, iinfo )
676  IF( iinfo .NE. 0 ) THEN
677  info = -4
678  RETURN
679  END IF
680 * DLARRB computes all gaps correctly except for the last one
681 * Record distance to VU/GU
682  wgap( wend ) = max( zero,
683  $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
684  DO 138 i = indl, indu
685  m = m + 1
686  iblock(m) = jblk
687  indexw(m) = i
688  138 CONTINUE
689  ELSE
690 * Call dqds to get all eigs (and then possibly delete unwanted
691 * eigenvalues).
692 * Note that dqds finds the eigenvalues of the L D L^T representation
693 * of T to high relative accuracy. High relative accuracy
694 * might be lost when the shift of the RRR is subtracted to obtain
695 * the eigenvalues of T. However, T is not guaranteed to define its
696 * eigenvalues to high relative accuracy anyway.
697 * Set RTOL to the order of the tolerance used in DLASQ2
698 * This is an ESTIMATED error, the worst case bound is 4*N*EPS
699 * which is usually too large and requires unnecessary work to be
700 * done by bisection when computing the eigenvectors
701  rtol = log(dble(in)) * four * eps
702  j = ibegin
703  DO 140 i = 1, in - 1
704  work( 2*i-1 ) = abs( d( j ) )
705  work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
706  j = j + 1
707  140 CONTINUE
708  work( 2*in-1 ) = abs( d( iend ) )
709  work( 2*in ) = zero
710  CALL dlasq2( in, work, iinfo )
711  IF( iinfo .NE. 0 ) THEN
712 * If IINFO = -5 then an index is part of a tight cluster
713 * and should be changed. The index is in IWORK(1) and the
714 * gap is in WORK(N+1)
715  info = -5
716  RETURN
717  ELSE
718 * Test that all eigenvalues are positive as expected
719  DO 149 i = 1, in
720  IF( work( i ).LT.zero ) THEN
721  info = -6
722  RETURN
723  ENDIF
724  149 CONTINUE
725  END IF
726  IF( sgndef.GT.zero ) THEN
727  DO 150 i = indl, indu
728  m = m + 1
729  w( m ) = work( in-i+1 )
730  iblock( m ) = jblk
731  indexw( m ) = i
732  150 CONTINUE
733  ELSE
734  DO 160 i = indl, indu
735  m = m + 1
736  w( m ) = -work( i )
737  iblock( m ) = jblk
738  indexw( m ) = i
739  160 CONTINUE
740  END IF
741 
742  DO 165 i = m - mb + 1, m
743 * the value of RTOL below should be the tolerance in DLASQ2
744  werr( i ) = rtol * abs( w(i) )
745  165 CONTINUE
746  DO 166 i = m - mb + 1, m - 1
747 * compute the right gap between the intervals
748  wgap( i ) = max( zero,
749  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
750  166 CONTINUE
751  wgap( m ) = max( zero,
752  $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
753  END IF
754 * proceed with next block
755  ibegin = iend + 1
756  wbegin = wend + 1
757  170 CONTINUE
758 *
759 
760  RETURN
761 *
762 * end of DLARRE2
763 *
764  END
max
#define max(A, B)
Definition: pcgemr.c:180
dlarre2
subroutine dlarre2(RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, DOL, DOU, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO)
Definition: dlarre2.f:6
min
#define min(A, B)
Definition: pcgemr.c:181