SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlarre2()

subroutine dlarre2 ( character  range,
integer  n,
double precision  vl,
double precision  vu,
integer  il,
integer  iu,
double precision, dimension( * )  d,
double precision, dimension( * )  e,
double precision, dimension( * )  e2,
double precision  rtol1,
double precision  rtol2,
double precision  spltol,
integer  nsplit,
integer, dimension( * )  isplit,
integer  m,
integer  dol,
integer  dou,
double precision, dimension( * )  w,
double precision, dimension( * )  werr,
double precision, dimension( * )  wgap,
integer, dimension( * )  iblock,
integer, dimension( * )  indexw,
double precision, dimension( * )  gers,
double precision  pivmin,
double precision, dimension( * )  work,
integer, dimension( * )  iwork,
integer  info 
)

Definition at line 1 of file dlarre2.f.

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*
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
dlamch
double precision function dlamch(cmach)
Definition tools.f:10
Here is the caller graph for this function: