4
#define abs(x) ((x) >= 0 ? (x) : -(x))
5
#define min(a,b) ((a) <= (b) ? (a) : (b))
6
#define max(a,b) ((a) >= (b) ? (a) : (b))
8
double dlamch_(char *cmach)
10
/* -- LAPACK auxiliary routine (version 2.0) --
11
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
12
Courant Institute, Argonne National Lab, and Rice University
18
DLAMCH determines double precision machine parameters.
23
CMACH (input) CHARACTER*1
24
Specifies the value to be returned by DLAMCH:
25
= 'E' or 'e', DLAMCH := eps
26
= 'S' or 's , DLAMCH := sfmin
27
= 'B' or 'b', DLAMCH := base
28
= 'P' or 'p', DLAMCH := eps*base
29
= 'N' or 'n', DLAMCH := t
30
= 'R' or 'r', DLAMCH := rnd
31
= 'M' or 'm', DLAMCH := emin
32
= 'U' or 'u', DLAMCH := rmin
33
= 'L' or 'l', DLAMCH := emax
34
= 'O' or 'o', DLAMCH := rmax
38
eps = relative machine precision
39
sfmin = safe minimum, such that 1/sfmin does not overflow
40
base = base of the machine
42
t = number of (base) digits in the mantissa
43
rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
44
emin = minimum exponent before (gradual) underflow
45
rmin = underflow threshold - base**(emin-1)
46
emax = largest exponent before overflow
47
rmax = overflow threshold - (base**emax)*(1-eps)
49
=====================================================================
52
static int first = TRUE_;
54
/* System generated locals */
57
/* Builtin functions */
58
double pow_di(double *, int *);
62
static double emin, prec, emax;
63
static int imin, imax;
65
static double rmin, rmax, t, rmach;
66
extern int lsame_(char *, char *);
67
static double small, sfmin;
68
extern /* Subroutine */ int dlamc2_(int *, int *, int *,
69
double *, int *, double *, int *, double *);
71
static double rnd, eps;
75
dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
81
eps = pow_di(&base, &i__1) / 2;
85
eps = pow_di(&base, &i__1);
94
/* Use SMALL plus a bit, to avoid the possibility of rounding
95
causing overflow when computing 1/sfmin. */
96
sfmin = small * (eps + 1.);
100
if (lsame_(cmach, "E")) {
102
} else if (lsame_(cmach, "S")) {
104
} else if (lsame_(cmach, "B")) {
106
} else if (lsame_(cmach, "P")) {
108
} else if (lsame_(cmach, "N")) {
110
} else if (lsame_(cmach, "R")) {
112
} else if (lsame_(cmach, "M")) {
114
} else if (lsame_(cmach, "U")) {
116
} else if (lsame_(cmach, "L")) {
118
} else if (lsame_(cmach, "O")) {
130
/* Subroutine */ int dlamc1_(int *beta, int *t, int *rnd, int
133
/* -- LAPACK auxiliary routine (version 2.0) --
134
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
135
Courant Institute, Argonne National Lab, and Rice University
142
DLAMC1 determines the machine parameters given by BETA, T, RND, and
149
The base of the machine.
152
The number of ( BETA ) digits in the mantissa.
155
Specifies whether proper rounding ( RND = .TRUE. ) or
156
chopping ( RND = .FALSE. ) occurs in addition. This may not
158
be a reliable guide to the way in which the machine performs
163
Specifies whether rounding appears to be done in the IEEE
164
'round to nearest' style.
169
The routine is based on the routine ENVRON by Malcolm and
170
incorporates suggestions by Gentleman and Marovich. See
172
Malcolm M. A. (1972) Algorithms to reveal properties of
173
floating-point arithmetic. Comms. of the ACM, 15, 949-951.
175
Gentleman W. M. and Marovich S. B. (1974) More on algorithms
176
that reveal properties of floating point arithmetic units.
177
Comms. of the ACM, 17, 276-277.
179
=====================================================================
181
/* Initialized data */
182
static int first = TRUE_;
183
/* System generated locals */
185
/* Local variables */
187
static double a, b, c, f;
190
extern double dlamc3_(double *, double *);
192
static double t1, t2;
194
static double one, qtr;
200
/* LBETA, LIEEE1, LT and LRND are the local values of BE
204
Throughout this routine we use the function DLAMC3 to ens
206
that relevant values are stored and not held in registers,
208
are not affected by optimizers.
210
Compute a = 2.0**m with the smallest positive integer m s
214
fl( a + 1.0 ) = a. */
219
/* + WHILE( C.EQ.ONE )LOOP */
223
c = dlamc3_(&a, &one);
225
c = dlamc3_(&c, &d__1);
230
Now compute b = 2.0**m with the smallest positive integer
234
fl( a + b ) .gt. a. */
239
/* + WHILE( C.EQ.A )LOOP */
248
Now compute the base. a and c are neighbouring floating po
250
numbers in the interval ( beta**t, beta**( t + 1 ) ) and
252
their difference is beta. Adding 0.25 to c is to ensure that
254
is truncated to beta and not ( beta - 1 ). */
259
c = dlamc3_(&c, &d__1);
260
lbeta = (int) (c + qtr);
262
/* Now determine whether rounding or chopping occurs, by addin
264
bit less than beta/2 and a bit more than beta/2 to
270
f = dlamc3_(&d__1, &d__2);
279
f = dlamc3_(&d__1, &d__2);
281
if (lrnd && c == a) {
285
/* Try and decide whether rounding is done in the IEEE 'round
287
nearest' style. B/2 is half a unit in the last place of the
289
numbers A and SAVEC. Furthermore, A is even, i.e. has last
291
zero, and SAVEC is odd. Thus adding B/2 to A should not cha
293
A, but adding B/2 to SAVEC should change SAVEC. */
296
t1 = dlamc3_(&d__1, &a);
298
t2 = dlamc3_(&d__1, &savec);
299
lieee1 = t1 == a && t2 > savec && lrnd;
301
/* Now find the mantissa, t. It should be the integer part
303
log to the base beta of a, however it is safer to determine
305
by powering. So we find t as the smallest positive integer
309
fl( beta**t + 1.0 ) = 1.0. */
315
/* + WHILE( C.EQ.ONE )LOOP */
320
c = dlamc3_(&a, &one);
322
c = dlamc3_(&c, &d__1);
340
/* Subroutine */ int dlamc2_(int *beta, int *t, int *rnd,
341
double *eps, int *emin, double *rmin, int *emax,
344
/* -- LAPACK auxiliary routine (version 2.0) --
345
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
346
Courant Institute, Argonne National Lab, and Rice University
353
DLAMC2 determines the machine parameters specified in its argument
360
The base of the machine.
363
The number of ( BETA ) digits in the mantissa.
366
Specifies whether proper rounding ( RND = .TRUE. ) or
367
chopping ( RND = .FALSE. ) occurs in addition. This may not
369
be a reliable guide to the way in which the machine performs
373
EPS (output) DOUBLE PRECISION
374
The smallest positive number such that
376
fl( 1.0 - EPS ) .LT. 1.0,
378
where fl denotes the computed value.
381
The minimum exponent before (gradual) underflow occurs.
383
RMIN (output) DOUBLE PRECISION
384
The smallest normalized number for the machine, given by
385
BASE**( EMIN - 1 ), where BASE is the floating point value
390
The maximum exponent before overflow occurs.
392
RMAX (output) DOUBLE PRECISION
393
The largest positive number for the machine, given by
394
BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
401
The computation of EPS is based on a routine PARANOIA by
402
W. Kahan of the University of California at Berkeley.
404
=====================================================================
406
/* Table of constant values */
409
/* Initialized data */
410
static int first = TRUE_;
411
static int iwarn = FALSE_;
412
/* System generated locals */
414
double d__1, d__2, d__3, d__4, d__5;
415
/* Builtin functions */
416
double pow_di(double *, int *);
417
/* Local variables */
421
static double leps, zero, a, b, c;
424
static int lemin, lemax, gnmin;
427
static double third, lrmin, lrmax, sixth;
428
extern /* Subroutine */ int dlamc1_(int *, int *, int *,
430
extern double dlamc3_(double *, double *);
432
extern /* Subroutine */ int dlamc4_(int *, double *, int *),
433
dlamc5_(int *, int *, int *, int *, int *,
435
static int lt, ngnmin, ngpmin;
436
static double one, two;
444
/* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values
446
BETA, T, RND, EPS, EMIN and RMIN.
448
Throughout this routine we use the function DLAMC3 to ens
450
that relevant values are stored and not held in registers,
452
are not affected by optimizers.
454
DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
457
dlamc1_(&lbeta, <, &lrnd, &lieee1);
459
/* Start to find EPS. */
463
a = pow_di(&b, &i__1);
466
/* Try some tricks to see whether or not this is the correct E
472
sixth = dlamc3_(&b, &d__1);
473
third = dlamc3_(&sixth, &sixth);
475
b = dlamc3_(&third, &d__1);
476
b = dlamc3_(&b, &sixth);
484
/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
486
if (leps > b && b > zero) {
489
/* Computing 5th power */
490
d__3 = two, d__4 = d__3, d__3 *= d__3;
491
/* Computing 2nd power */
493
d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
494
c = dlamc3_(&d__1, &d__2);
496
c = dlamc3_(&half, &d__1);
497
b = dlamc3_(&half, &c);
499
c = dlamc3_(&half, &d__1);
500
b = dlamc3_(&half, &c);
509
/* Computation of EPS complete.
511
Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3
513
Keep dividing A by BETA until (gradual) underflow occurs. T
515
is detected when we cannot recover the previous A. */
519
for (i = 1; i <= 3; ++i) {
520
d__1 = small * rbase;
521
small = dlamc3_(&d__1, &zero);
524
a = dlamc3_(&one, &small);
525
dlamc4_(&ngpmin, &one, &lbeta);
527
dlamc4_(&ngnmin, &d__1, &lbeta);
528
dlamc4_(&gpmin, &a, &lbeta);
530
dlamc4_(&gnmin, &d__1, &lbeta);
533
if (ngpmin == ngnmin && gpmin == gnmin) {
534
if (ngpmin == gpmin) {
536
/* ( Non twos-complement machines, no gradual under
539
} else if (gpmin - ngpmin == 3) {
540
lemin = ngpmin - 1 + lt;
542
/* ( Non twos-complement machines, with gradual und
544
e.g., IEEE standard followers ) */
546
lemin = min(ngpmin,gpmin);
547
/* ( A guess; no known machine ) */
551
} else if (ngpmin == gpmin && ngnmin == gnmin) {
552
if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
553
lemin = max(ngpmin,ngnmin);
554
/* ( Twos-complement machines, no gradual underflow
558
lemin = min(ngpmin,ngnmin);
559
/* ( A guess; no known machine ) */
563
} else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
565
if (gpmin - min(ngpmin,ngnmin) == 3) {
566
lemin = max(ngpmin,ngnmin) - 1 + lt;
567
/* ( Twos-complement machines with gradual underflo
569
no known machine ) */
571
lemin = min(ngpmin,ngnmin);
572
/* ( A guess; no known machine ) */
578
i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin);
579
lemin = min(i__1,gnmin);
580
/* ( A guess; no known machine ) */
584
Comment out this if block if EMIN is ok */
587
printf("\n\n WARNING. The value EMIN may be incorrect:- ");
588
printf("EMIN = %8i\n",lemin);
589
printf("If, after inspection, the value EMIN looks acceptable");
590
printf("please comment out \n the IF block as marked within the");
591
printf("code of routine DLAMC2, \n otherwise supply EMIN");
592
printf("explicitly.\n");
596
Assume IEEE arithmetic if we found denormalised numbers abo
598
or if arithmetic seems to round in the IEEE style, determi
600
in routine DLAMC1. A true IEEE machine should have both thi
602
true; however, faulty machines may have one or the other. */
604
ieee = ieee || lieee1;
606
/* Compute RMIN by successive division by BETA. We could comp
608
RMIN as BASE**( EMIN - 1 ), but some machines underflow dur
614
for (i = 1; i <= 1-lemin; ++i) {
615
d__1 = lrmin * rbase;
616
lrmin = dlamc3_(&d__1, &zero);
620
/* Finally, call DLAMC5 to compute EMAX and RMAX. */
622
dlamc5_(&lbeta, <, &lemin, &ieee, &lemax, &lrmax);
642
double dlamc3_(double *a, double *b)
644
/* -- LAPACK auxiliary routine (version 2.0) --
645
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
646
Courant Institute, Argonne National Lab, and Rice University
653
DLAMC3 is intended to force A and B to be stored prior to doing
655
the addition of A and B , for use in situations where optimizers
657
might hold one of these in a register.
662
A, B (input) DOUBLE PRECISION
665
=====================================================================
668
System generated locals */
680
/* Subroutine */ int dlamc4_(int *emin, double *start, int *base)
682
/* -- LAPACK auxiliary routine (version 2.0) --
683
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
684
Courant Institute, Argonne National Lab, and Rice University
691
DLAMC4 is a service routine for DLAMC2.
697
The minimum exponent before (gradual) underflow, computed by
699
setting A = START and dividing by BASE until the previous A
700
can not be recovered.
702
START (input) DOUBLE PRECISION
703
The starting point for determining EMIN.
706
The base of the machine.
708
=====================================================================
710
/* System generated locals */
713
/* Local variables */
714
static double zero, a;
716
static double rbase, b1, b2, c1, c2, d1, d2;
717
extern double dlamc3_(double *, double *);
726
b1 = dlamc3_(&d__1, &zero);
731
/* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
732
$ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
734
if (c1 == a && c2 == a && d1 == a && d2 == a) {
738
b1 = dlamc3_(&d__1, &zero);
740
c1 = dlamc3_(&d__1, &zero);
743
for (i = 1; i <= *base; ++i) {
748
b2 = dlamc3_(&d__1, &zero);
750
c2 = dlamc3_(&d__1, &zero);
753
for (i = 1; i <= *base; ++i) {
768
/* Subroutine */ int dlamc5_(int *beta, int *p, int *emin,
769
int *ieee, int *emax, double *rmax)
771
/* -- LAPACK auxiliary routine (version 2.0) --
772
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
773
Courant Institute, Argonne National Lab, and Rice University
780
DLAMC5 attempts to compute RMAX, the largest machine floating-point
781
number, without overflow. It assumes that EMAX + abs(EMIN) sum
782
approximately to a power of 2. It will fail on machines where this
783
assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
785
EMAX = 28718). It will also fail if the value supplied for EMIN is
786
too large (i.e. too close to zero), probably with overflow.
792
The base of floating-point arithmetic.
795
The number of base BETA digits in the mantissa of a
796
floating-point value.
799
The minimum exponent before (gradual) underflow.
802
A int flag specifying whether or not the arithmetic
803
system is thought to comply with the IEEE standard.
806
The largest exponent before overflow
808
RMAX (output) DOUBLE PRECISION
809
The largest machine floating-point number.
811
=====================================================================
815
First compute LEXP and UEXP, two powers of 2 that bound
816
abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
817
approximately to the bound that is closest to abs(EMIN).
818
(EMAX is the exponent of the required number RMAX). */
819
/* Table of constant values */
820
static double c_b5 = 0.;
822
/* System generated locals */
825
/* Local variables */
831
extern double dlamc3_(double *, double *);
832
static double recbas;
833
static int exbits, expsum, try__;
841
if (try__ <= -(*emin)) {
846
if (lexp == -(*emin)) {
853
/* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
854
than or equal to EMIN. EXBITS is the number of bits needed to
855
store the exponent. */
857
if (uexp + *emin > -lexp - *emin) {
863
/* EXPSUM is the exponent range, approximately equal to
866
*emax = expsum + *emin - 1;
867
nbits = exbits + 1 + *p;
869
/* NBITS is the total number of bits needed to store a
870
floating-point number. */
872
if (nbits % 2 == 1 && *beta == 2) {
874
/* Either there are an odd number of bits used to store a
875
floating-point number, which is unlikely, or some bits are
877
not used in the representation of numbers, which is possible
879
(e.g. Cray machines) or the mantissa has an implicit bit,
880
(e.g. IEEE machines, Dec Vax machines), which is perhaps the
882
most likely. We have to assume the last alternative.
883
If this is true, then we need to reduce EMAX by one because
885
there must be some way of representing zero in an implicit-b
887
system. On machines like Cray, we are reducing EMAX by one
896
/* Assume we are on an IEEE machine which reserves one exponent
898
for infinity and NaN. */
903
/* Now create RMAX, the largest machine number, which should
904
be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
906
First compute 1.0 - BETA**(-P), being careful that the
907
result is less than 1.0 . */
913
for (i = 1; i <= *p; ++i) {
925
/* Now multiply by BETA**EMAX to get RMAX. */
928
for (i = 1; i <= *emax; ++i) {
930
y = dlamc3_(&d__1, &c_b5);
941
double pow_di(double *ap, int *bp)
added
removed
Lib/sparse/SuperLU/SRC/dlamch.c
