SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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
min
#define min(A, B)
Definition pcgemr.c:181
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