OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dlarrb2.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dlarrb2 (n, d, lld, ifirst, ilast, rtol1, rtol2, offset, w, wgap, werr, work, iwork, pivmin, lgpvmn, lgspdm, twist, info)
integer function dlaneg2 (n, d, lld, sigma, pivmin, r)
integer function dlaneg2a (n, dlld, sigma, pivmin, r)

Function/Subroutine Documentation

◆ dlaneg2()

integer function dlaneg2 ( integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) lld,
double precision sigma,
double precision pivmin,
integer r )

Definition at line 336 of file dlarrb2.f.

337*
338 IMPLICIT NONE
339*
340 INTEGER DLANEG2
341*
342* .. Scalar Arguments ..
343 INTEGER N, R
344 DOUBLE PRECISION PIVMIN, SIGMA
345* ..
346* .. Array Arguments ..
347 DOUBLE PRECISION D( * ), LLD( * )
348*
349 DOUBLE PRECISION ZERO
350 parameter( zero = 0.0d0 )
351
352 INTEGER BLKLEN
353 parameter( blklen = 2048 )
354* ..
355* .. Local Scalars ..
356 INTEGER BJ, J, NEG1, NEG2, NEGCNT, TO
357 DOUBLE PRECISION DMINUS, DPLUS, GAMMA, P, S, T, TMP, XSAV
358 LOGICAL SAWNAN
359* ..
360* .. External Functions ..
361 LOGICAL DISNAN
362 EXTERNAL disnan
363
364 negcnt = 0
365*
366* I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
367* run dstqds block-wise to avoid excessive work when NaNs occur
368*
369 s = zero
370 DO 210 bj = 1, r-1, blklen
371 neg1 = 0
372 xsav = s
373 to = bj+blklen-1
374 IF ( to.LE.r-1 ) THEN
375 DO 21 j = bj, to
376 t = s - sigma
377 dplus = d( j ) + t
378 IF( dplus.LT.zero ) neg1=neg1 + 1
379 s = t*lld( j ) / dplus
380 21 CONTINUE
381 ELSE
382 DO 22 j = bj, r-1
383 t = s - sigma
384 dplus = d( j ) + t
385 IF( dplus.LT.zero ) neg1=neg1 + 1
386 s = t*lld( j ) / dplus
387 22 CONTINUE
388 ENDIF
389 sawnan = disnan( s )
390*
391 IF( sawnan ) THEN
392 neg1 = 0
393 s = xsav
394 to = bj+blklen-1
395 IF ( to.LE.r-1 ) THEN
396 DO 23 j = bj, to
397 t = s - sigma
398 dplus = d( j ) + t
399 IF(abs(dplus).LT.pivmin)
400 $ dplus = -pivmin
401 tmp = lld( j ) / dplus
402 IF( dplus.LT.zero )
403 $ neg1 = neg1 + 1
404 s = t*tmp
405 IF( tmp.EQ.zero ) s = lld( j )
406 23 CONTINUE
407 ELSE
408 DO 24 j = bj, r-1
409 t = s - sigma
410 dplus = d( j ) + t
411 IF(abs(dplus).LT.pivmin)
412 $ dplus = -pivmin
413 tmp = lld( j ) / dplus
414 IF( dplus.LT.zero ) neg1=neg1+1
415 s = t*tmp
416 IF( tmp.EQ.zero ) s = lld( j )
417 24 CONTINUE
418 ENDIF
419 END IF
420 negcnt = negcnt + neg1
421 210 CONTINUE
422*
423* II) lower part: L D L^T - SIGMA I = U- D- U-^T
424*
425 p = d( n ) - sigma
426 DO 230 bj = n-1, r, -blklen
427 neg2 = 0
428 xsav = p
429 to = bj-blklen+1
430 IF ( to.GE.r ) THEN
431 DO 25 j = bj, to, -1
432 dminus = lld( j ) + p
433 IF( dminus.LT.zero ) neg2=neg2+1
434 tmp = p / dminus
435 p = tmp * d( j ) - sigma
436 25 CONTINUE
437 ELSE
438 DO 26 j = bj, r, -1
439 dminus = lld( j ) + p
440 IF( dminus.LT.zero ) neg2=neg2+1
441 tmp = p / dminus
442 p = tmp * d( j ) - sigma
443 26 CONTINUE
444 ENDIF
445 sawnan = disnan( p )
446*
447 IF( sawnan ) THEN
448 neg2 = 0
449 p = xsav
450 to = bj-blklen+1
451 IF ( to.GE.r ) THEN
452 DO 27 j = bj, to, -1
453 dminus = lld( j ) + p
454 IF(abs(dminus).LT.pivmin)
455 $ dminus = -pivmin
456 tmp = d( j ) / dminus
457 IF( dminus.LT.zero )
458 $ neg2 = neg2 + 1
459 p = p*tmp - sigma
460 IF( tmp.EQ.zero )
461 $ p = d( j ) - sigma
462 27 CONTINUE
463 ELSE
464 DO 28 j = bj, r, -1
465 dminus = lld( j ) + p
466 IF(abs(dminus).LT.pivmin)
467 $ dminus = -pivmin
468 tmp = d( j ) / dminus
469 IF( dminus.LT.zero )
470 $ neg2 = neg2 + 1
471 p = p*tmp - sigma
472 IF( tmp.EQ.zero )
473 $ p = d( j ) - sigma
474 28 CONTINUE
475 ENDIF
476 END IF
477 negcnt = negcnt + neg2
478 230 CONTINUE
479*
480* III) Twist index
481*
482 gamma = s + p
483 IF( gamma.LT.zero ) negcnt = negcnt+1
484
485 dlaneg2 = negcnt
dlaneg2
integer function dlaneg2(n, d, lld, sigma, pivmin, r)
Definition dlarrb2.f:337
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:59

◆ dlaneg2a()

integer function dlaneg2a ( integer n,
double precision, dimension( * ) dlld,
double precision sigma,
double precision pivmin,
integer r )

Definition at line 490 of file dlarrb2.f.

491*
492 IMPLICIT NONE
493*
494 INTEGER DLANEG2A
495*
496* .. Scalar Arguments ..
497 INTEGER N, R
498 DOUBLE PRECISION PIVMIN, SIGMA
499* ..
500* .. Array Arguments ..
501 DOUBLE PRECISION DLLD( * )
502*
503 DOUBLE PRECISION ZERO
504 parameter( zero = 0.0d0 )
505
506 INTEGER BLKLEN
507 parameter( blklen = 512 )
508*
509* ..
510* .. Intrinsic Functions ..
511 INTRINSIC int
512* ..
513* .. Local Scalars ..
514 INTEGER BJ, I, J, NB, NEG1, NEG2, NEGCNT, NX
515 DOUBLE PRECISION DMINUS, DPLUS, GAMMA, P, S, T, TMP, XSAV
516 LOGICAL SAWNAN
517* ..
518* .. External Functions ..
519 LOGICAL DISNAN
520 EXTERNAL disnan
521
522 negcnt = 0
523*
524* I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
525* run dstqds block-wise to avoid excessive work when NaNs occur,
526* first in chunks of size BLKLEN and then the remainder
527*
528 nb = int((r-1)/blklen)
529 nx = nb*blklen
530 s = zero
531 DO 210 bj = 1, nx, blklen
532 neg1 = 0
533 xsav = s
534 DO 21 j = bj, bj+blklen-1
535 i = 2*j
536 t = s - sigma
537 dplus = dlld( i-1 ) + t
538 IF( dplus.LT.zero ) neg1=neg1 + 1
539 s = t*dlld( i ) / dplus
540 21 CONTINUE
541 sawnan = disnan( s )
542*
543 IF( sawnan ) THEN
544 neg1 = 0
545 s = xsav
546 DO 23 j = bj, bj+blklen-1
547 i = 2*j
548 t = s - sigma
549 dplus = dlld( i-1 ) + t
550 IF(abs(dplus).LT.pivmin)
551 $ dplus = -pivmin
552 tmp = dlld( i ) / dplus
553 IF( dplus.LT.zero )
554 $ neg1 = neg1 + 1
555 s = t*tmp
556 IF( tmp.EQ.zero ) s = dlld( i )
557 23 CONTINUE
558 END IF
559 negcnt = negcnt + neg1
560 210 CONTINUE
561*
562 neg1 = 0
563 xsav = s
564 DO 22 j = nx+1, r-1
565 i = 2*j
566 t = s - sigma
567 dplus = dlld( i-1 ) + t
568 IF( dplus.LT.zero ) neg1=neg1 + 1
569 s = t*dlld( i ) / dplus
570 22 CONTINUE
571 sawnan = disnan( s )
572*
573 IF( sawnan ) THEN
574 neg1 = 0
575 s = xsav
576 DO 24 j = nx+1, r-1
577 i = 2*j
578 t = s - sigma
579 dplus = dlld( i-1 ) + t
580 IF(abs(dplus).LT.pivmin)
581 $ dplus = -pivmin
582 tmp = dlld( i ) / dplus
583 IF( dplus.LT.zero ) neg1=neg1+1
584 s = t*tmp
585 IF( tmp.EQ.zero ) s = dlld( i )
586 24 CONTINUE
587 ENDIF
588 negcnt = negcnt + neg1
589*
590* II) lower part: L D L^T - SIGMA I = U- D- U-^T
591*
592 nb = int((n-r)/blklen)
593 nx = n-nb*blklen
594 p = dlld( 2*n-1 ) - sigma
595 DO 230 bj = n-1, nx, -blklen
596 neg2 = 0
597 xsav = p
598 DO 25 j = bj, bj-blklen+1, -1
599 i = 2*j
600 dminus = dlld( i ) + p
601 IF( dminus.LT.zero ) neg2=neg2+1
602 tmp = p / dminus
603 p = tmp * dlld( i-1 ) - sigma
604 25 CONTINUE
605 sawnan = disnan( p )
606*
607 IF( sawnan ) THEN
608 neg2 = 0
609 p = xsav
610 DO 27 j = bj, bj-blklen+1, -1
611 i = 2*j
612 dminus = dlld( i ) + p
613 IF(abs(dminus).LT.pivmin)
614 $ dminus = -pivmin
615 tmp = dlld( i-1 ) / dminus
616 IF( dminus.LT.zero )
617 $ neg2 = neg2 + 1
618 p = p*tmp - sigma
619 IF( tmp.EQ.zero )
620 $ p = dlld( i-1 ) - sigma
621 27 CONTINUE
622 END IF
623 negcnt = negcnt + neg2
624 230 CONTINUE
625
626 neg2 = 0
627 xsav = p
628 DO 26 j = nx-1, r, -1
629 i = 2*j
630 dminus = dlld( i ) + p
631 IF( dminus.LT.zero ) neg2=neg2+1
632 tmp = p / dminus
633 p = tmp * dlld( i-1 ) - sigma
634 26 CONTINUE
635 sawnan = disnan( p )
636*
637 IF( sawnan ) THEN
638 neg2 = 0
639 p = xsav
640 DO 28 j = nx-1, r, -1
641 i = 2*j
642 dminus = dlld( i ) + p
643 IF(abs(dminus).LT.pivmin)
644 $ dminus = -pivmin
645 tmp = dlld( i-1 ) / dminus
646 IF( dminus.LT.zero )
647 $ neg2 = neg2 + 1
648 p = p*tmp - sigma
649 IF( tmp.EQ.zero )
650 $ p = dlld( i-1 ) - sigma
651 28 CONTINUE
652 END IF
653 negcnt = negcnt + neg2
654*
655* III) Twist index
656*
657 gamma = s + p
658 IF( gamma.LT.zero ) negcnt = negcnt+1
659
660 dlaneg2a = negcnt
dlaneg2a
integer function dlaneg2a(n, dlld, sigma, pivmin, r)
Definition dlarrb2.f:491

◆ dlarrb2()

subroutine dlarrb2 ( integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) lld,
integer ifirst,
integer ilast,
double precision rtol1,
double precision rtol2,
integer offset,
double precision, dimension( * ) w,
double precision, dimension( * ) wgap,
double precision, dimension( * ) werr,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
double precision pivmin,
double precision lgpvmn,
double precision lgspdm,
integer twist,
integer info )

Definition at line 1 of file dlarrb2.f.

4*
5* -- ScaLAPACK auxiliary routine (version 2.0) --
6* Univ. of Tennessee, Univ. of California Berkeley, Univ of Colorado Denver
7* July 4, 2010
8*
9 IMPLICIT NONE
10*
11* .. Scalar Arguments ..
12 INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
13 DOUBLE PRECISION LGPVMN, LGSPDM, PIVMIN,
14 $ RTOL1, RTOL2
15* ..
16* .. Array Arguments ..
17 INTEGER IWORK( * )
18 DOUBLE PRECISION D( * ), LLD( * ), W( * ),
19 $ WERR( * ), WGAP( * ), WORK( * )
20* ..
21*
22* Purpose
23* =======
24*
25* Given the relatively robust representation(RRR) L D L^T, DLARRB2
26* does "limited" bisection to refine the eigenvalues of L D L^T,
27* W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
28* guesses for these eigenvalues are input in W, the corresponding estimate
29* of the error in these guesses and their gaps are input in WERR
30* and WGAP, respectively. During bisection, intervals
31* [left, right] are maintained by storing their mid-points and
32* semi-widths in the arrays W and WERR respectively.
33*
34* NOTE:
35* There are very few minor differences between DLARRB from LAPACK
36* and this current subroutine DLARRB2.
37* The most important reason for creating this nearly identical copy
38* is profiling: in the ScaLAPACK MRRR algorithm, eigenvalue computation
39* using DLARRB2 is used for refinement in the construction of
40* the representation tree, as opposed to the initial computation of the
41* eigenvalues for the root RRR which uses DLARRB. When profiling,
42* this allows an easy quantification of refinement work vs. computing
43* eigenvalues of the root.
44*
45* Arguments
46* =========
47*
48* N (input) INTEGER
49* The order of the matrix.
50*
51* D (input) DOUBLE PRECISION array, dimension (N)
52* The N diagonal elements of the diagonal matrix D.
53*
54* LLD (input) DOUBLE PRECISION array, dimension (N-1)
55* The (N-1) elements L(i)*L(i)*D(i).
56*
57* IFIRST (input) INTEGER
58* The index of the first eigenvalue to be computed.
59*
60* ILAST (input) INTEGER
61* The index of the last eigenvalue to be computed.
62*
63* RTOL1 (input) DOUBLE PRECISION
64* RTOL2 (input) DOUBLE PRECISION
65* Tolerance for the convergence of the bisection intervals.
66* An interval [LEFT,RIGHT] has converged if
67* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
68* where GAP is the (estimated) distance to the nearest
69* eigenvalue.
70*
71* OFFSET (input) INTEGER
72* Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
73* through ILAST-OFFSET elements of these arrays are to be used.
74*
75* W (input/output) DOUBLE PRECISION array, dimension (N)
76* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
77* estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST.
78* On output, these estimates are refined.
79*
80* WGAP (input/output) DOUBLE PRECISION array, dimension (N-1)
81* On input, the (estimated) gaps between consecutive
82* eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
83* eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
84* then WGAP(IFIRST-OFFSET) must be set to ZERO.
85* On output, these gaps are refined.
86*
87* WERR (input/output) DOUBLE PRECISION array, dimension (N)
88* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
89* the errors in the estimates of the corresponding elements in W.
90* On output, these errors are refined.
91*
92* WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
93* Workspace.
94*
95* IWORK (workspace) INTEGER array, dimension (2*N)
96* Workspace.
97*
98* PIVMIN (input) DOUBLE PRECISION
99* The minimum pivot in the sturm sequence.
100*
101* LGPVMN (input) DOUBLE PRECISION
102* Logarithm of PIVMIN, precomputed.
103*
104* LGSPDM (input) DOUBLE PRECISION
105* Logarithm of the spectral diameter, precomputed.
106*
107* TWIST (input) INTEGER
108* The twist index for the twisted factorization that is used
109* for the negcount.
110* TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
111* TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
112* TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
113*
114* INFO (output) INTEGER
115* Error flag.
116*
117* .. Parameters ..
118 DOUBLE PRECISION ZERO, TWO, HALF
119 parameter( zero = 0.0d0, two = 2.0d0,
120 $ half = 0.5d0 )
121 INTEGER MAXITR
122* ..
123* .. Local Scalars ..
124 INTEGER I, I1, II, INDLLD, IP, ITER, J, K, NEGCNT,
125 $ NEXT, NINT, OLNINT, PREV, R
126 DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
127 $ RGAP, RIGHT, SAVGAP, TMP, WIDTH
128 LOGICAL PARANOID
129* ..
130* .. External Functions ..
131 LOGICAL DISNAN
132 DOUBLE PRECISION DLAMCH
133 INTEGER DLANEG2A
134 EXTERNAL disnan, dlamch,
135 $ dlaneg2a
136*
137* ..
138* .. Intrinsic Functions ..
139 INTRINSIC abs, max, min
140* ..
141* .. Executable Statements ..
142*
143 info = 0
144*
145* Turn on paranoid check for rounding errors
146* invalidating uncertainty intervals of eigenvalues
147*
148 paranoid = .true.
149*
150 maxitr = int( ( lgspdm - lgpvmn ) / log( two ) ) + 2
151 mnwdth = two * pivmin
152*
153 r = twist
154*
155 indlld = 2*n
156 DO 5 j = 1, n-1
157 i=2*j
158 work(indlld+i-1) = d(j)
159 work(indlld+i) = lld(j)
160 5 CONTINUE
161 work(indlld+2*n-1) = d(n)
162*
163 IF((r.LT.1).OR.(r.GT.n)) r = n
164*
165* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
166* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
167* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
168* for an unconverged interval is set to the index of the next unconverged
169* interval, and is -1 or 0 for a converged interval. Thus a linked
170* list of unconverged intervals is set up.
171*
172 i1 = ifirst
173* The number of unconverged intervals
174 nint = 0
175* The last unconverged interval found
176 prev = 0
177
178 rgap = wgap( i1-offset )
179 DO 75 i = i1, ilast
180 k = 2*i
181 ii = i - offset
182 left = w( ii ) - werr( ii )
183 right = w( ii ) + werr( ii )
184 lgap = rgap
185 rgap = wgap( ii )
186 gap = min( lgap, rgap )
187
188 IF((abs(left).LE.16*pivmin).OR.(abs(right).LE.16*pivmin))
189 $ THEN
190 info = -1
191 RETURN
192 ENDIF
193
194 IF( paranoid ) THEN
195* Make sure that [LEFT,RIGHT] contains the desired eigenvalue
196* Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
197*
198* Do while( NEGCNT(LEFT).GT.I-1 )
199*
200 back = werr( ii )
201 20 CONTINUE
202 negcnt = dlaneg2a( n, work(indlld+1), left, pivmin, r )
203 IF( negcnt.GT.i-1 ) THEN
204 left = left - back
205 back = two*back
206 GO TO 20
207 END IF
208*
209* Do while( NEGCNT(RIGHT).LT.I )
210* Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
211*
212 back = werr( ii )
213 50 CONTINUE
214 negcnt = dlaneg2a( n, work(indlld+1),right, pivmin, r )
215
216 IF( negcnt.LT.i ) THEN
217 right = right + back
218 back = two*back
219 GO TO 50
220 END IF
221 ENDIF
222
223 width = half*abs( left - right )
224 tmp = max( abs( left ), abs( right ) )
225 cvrgd = max(rtol1*gap,rtol2*tmp)
226 IF( width.LE.cvrgd .OR. width.LE.mnwdth ) THEN
227* This interval has already converged and does not need refinement.
228* (Note that the gaps might change through refining the
229* eigenvalues, however, they can only get bigger.)
230* Remove it from the list.
231 iwork( k-1 ) = -1
232* Make sure that I1 always points to the first unconverged interval
233 IF((i.EQ.i1).AND.(i.LT.ilast)) i1 = i + 1
234 IF((prev.GE.i1).AND.(i.LE.ilast)) iwork( 2*prev-1 ) = i + 1
235 ELSE
236* unconverged interval found
237 prev = i
238 nint = nint + 1
239 iwork( k-1 ) = i + 1
240 iwork( k ) = negcnt
241 END IF
242 work( k-1 ) = left
243 work( k ) = right
244 75 CONTINUE
245
246*
247* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
248* and while (ITER.LT.MAXITR)
249*
250 iter = 0
251 80 CONTINUE
252 prev = i1 - 1
253 i = i1
254 olnint = nint
255
256 DO 100 ip = 1, olnint
257 k = 2*i
258 ii = i - offset
259 rgap = wgap( ii )
260 lgap = rgap
261 IF(ii.GT.1) lgap = wgap( ii-1 )
262 gap = min( lgap, rgap )
263 next = iwork( k-1 )
264 left = work( k-1 )
265 right = work( k )
266 mid = half*( left + right )
267* semiwidth of interval
268 width = right - mid
269 tmp = max( abs( left ), abs( right ) )
270 cvrgd = max(rtol1*gap,rtol2*tmp)
271 IF( ( width.LE.cvrgd ) .OR. ( width.LE.mnwdth ).OR.
272 $ ( iter.EQ.maxitr ) )THEN
273* reduce number of unconverged intervals
274 nint = nint - 1
275* Mark interval as converged.
276 iwork( k-1 ) = 0
277 IF( i1.EQ.i ) THEN
278 i1 = next
279 ELSE
280* Prev holds the last unconverged interval previously examined
281 IF(prev.GE.i1) iwork( 2*prev-1 ) = next
282 END IF
283 i = next
284 GO TO 100
285 END IF
286 prev = i
287*
288* Perform one bisection step
289*
290 negcnt = dlaneg2a( n, work(indlld+1), mid, pivmin, r )
291 IF( negcnt.LE.i-1 ) THEN
292 work( k-1 ) = mid
293 ELSE
294 work( k ) = mid
295 END IF
296 i = next
297 100 CONTINUE
298 iter = iter + 1
299* do another loop if there are still unconverged intervals
300* However, in the last iteration, all intervals are accepted
301* since this is the best we can do.
302 IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
303*
304*
305* At this point, all the intervals have converged
306*
307* save this gap to restore it after the loop
308 savgap = wgap( ilast-offset )
309*
310 left = work( 2*ifirst-1 )
311 DO 110 i = ifirst, ilast
312 k = 2*i
313 ii = i - offset
314* RIGHT is the right boundary of this current interval
315 right = work( k )
316* All intervals marked by '0' have been refined.
317 IF( iwork( k-1 ).EQ.0 ) THEN
318 w( ii ) = half*( left+right )
319 werr( ii ) = right - w( ii )
320 END IF
321* Left is the boundary of the next interval
322 left = work( k +1 )
323 wgap( ii ) = max( zero, left - right )
324 110 CONTINUE
325* restore the last gap which was overwritten by garbage
326 wgap( ilast-offset ) = savgap
327
328 RETURN
329*
330* End of DLARRB2
331*
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21