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