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