Mercurial > hg > qm-dsp

comparison ext/cblas/src/dlamch.c @ 430:335af74a25b6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge from branch clapack-included
author Chris Cannam <c.cannam@qmul.ac.uk>
date Fri, 30 Sep 2016 16:24:24 +0100
parents 905e45637745
children
comparison
equal deleted inserted replaced
426:a23b9f8b4a59 430:335af74a25b6
1 /* dlamch.f -- translated by f2c (version 20061008).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #include "f2c.h"
14 #include "blaswrap.h"
15
16 /* Table of constant values */
17
18 static integer c__1 = 1;
19 static doublereal c_b32 = 0.;
20
21 doublereal dlamch_(char *cmach)
22 {
23 /* Initialized data */
24
25 static logical first = TRUE_;
26
27 /* System generated locals */
28 integer i__1;
29 doublereal ret_val;
30
31 /* Builtin functions */
32 double pow_di(doublereal *, integer *);
33
34 /* Local variables */
35 static doublereal t;
36 integer it;
37 static doublereal rnd, eps, base;
38 integer beta;
39 static doublereal emin, prec, emax;
40 integer imin, imax;
41 logical lrnd;
42 static doublereal rmin, rmax;
43 doublereal rmach;
44 extern logical lsame_(char *, char *);
45 doublereal small;
46 static doublereal sfmin;
47 extern /* Subroutine */ int dlamc2_(integer *, integer *, logical *,
48 doublereal *, integer *, doublereal *, integer *, doublereal *);
49
50
51 /* -- LAPACK auxiliary routine (version 3.2) -- */
52 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
53 /* November 2006 */
54
55 /* .. Scalar Arguments .. */
56 /* .. */
57
58 /* Purpose */
59 /* ======= */
60
61 /* DLAMCH determines double precision machine parameters. */
62
63 /* Arguments */
64 /* ========= */
65
66 /* CMACH (input) CHARACTER*1 */
67 /* Specifies the value to be returned by DLAMCH: */
68 /* = 'E' or 'e', DLAMCH := eps */
69 /* = 'S' or 's , DLAMCH := sfmin */
70 /* = 'B' or 'b', DLAMCH := base */
71 /* = 'P' or 'p', DLAMCH := eps*base */
72 /* = 'N' or 'n', DLAMCH := t */
73 /* = 'R' or 'r', DLAMCH := rnd */
74 /* = 'M' or 'm', DLAMCH := emin */
75 /* = 'U' or 'u', DLAMCH := rmin */
76 /* = 'L' or 'l', DLAMCH := emax */
77 /* = 'O' or 'o', DLAMCH := rmax */
78
79 /* where */
80
81 /* eps = relative machine precision */
82 /* sfmin = safe minimum, such that 1/sfmin does not overflow */
83 /* base = base of the machine */
84 /* prec = eps*base */
85 /* t = number of (base) digits in the mantissa */
86 /* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */
87 /* emin = minimum exponent before (gradual) underflow */
88 /* rmin = underflow threshold - base**(emin-1) */
89 /* emax = largest exponent before overflow */
90 /* rmax = overflow threshold - (base**emax)*(1-eps) */
91
92 /* ===================================================================== */
93
94 /* .. Parameters .. */
95 /* .. */
96 /* .. Local Scalars .. */
97 /* .. */
98 /* .. External Functions .. */
99 /* .. */
100 /* .. External Subroutines .. */
101 /* .. */
102 /* .. Save statement .. */
103 /* .. */
104 /* .. Data statements .. */
105 /* .. */
106 /* .. Executable Statements .. */
107
108 if (first) {
109 dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
110 base = (doublereal) beta;
111 t = (doublereal) it;
112 if (lrnd) {
113 rnd = 1.;
114 i__1 = 1 - it;
115 eps = pow_di(&base, &i__1) / 2;
116 } else {
117 rnd = 0.;
118 i__1 = 1 - it;
119 eps = pow_di(&base, &i__1);
120 }
121 prec = eps * base;
122 emin = (doublereal) imin;
123 emax = (doublereal) imax;
124 sfmin = rmin;
125 small = 1. / rmax;
126 if (small >= sfmin) {
127
128 /* Use SMALL plus a bit, to avoid the possibility of rounding */
129 /* causing overflow when computing 1/sfmin. */
130
131 sfmin = small * (eps + 1.);
132 }
133 }
134
135 if (lsame_(cmach, "E")) {
136 rmach = eps;
137 } else if (lsame_(cmach, "S")) {
138 rmach = sfmin;
139 } else if (lsame_(cmach, "B")) {
140 rmach = base;
141 } else if (lsame_(cmach, "P")) {
142 rmach = prec;
143 } else if (lsame_(cmach, "N")) {
144 rmach = t;
145 } else if (lsame_(cmach, "R")) {
146 rmach = rnd;
147 } else if (lsame_(cmach, "M")) {
148 rmach = emin;
149 } else if (lsame_(cmach, "U")) {
150 rmach = rmin;
151 } else if (lsame_(cmach, "L")) {
152 rmach = emax;
153 } else if (lsame_(cmach, "O")) {
154 rmach = rmax;
155 }
156
157 ret_val = rmach;
158 first = FALSE_;
159 return ret_val;
160
161 /* End of DLAMCH */
162
163 } /* dlamch_ */
164
165
166 /* *********************************************************************** */
167
168 /* Subroutine */ int dlamc1_(integer *beta, integer *t, logical *rnd, logical
169 *ieee1)
170 {
171 /* Initialized data */
172
173 static logical first = TRUE_;
174
175 /* System generated locals */
176 doublereal d__1, d__2;
177
178 /* Local variables */
179 doublereal a, b, c__, f, t1, t2;
180 static integer lt;
181 doublereal one, qtr;
182 static logical lrnd;
183 static integer lbeta;
184 doublereal savec;
185 extern doublereal dlamc3_(doublereal *, doublereal *);
186 static logical lieee1;
187
188
189 /* -- LAPACK auxiliary routine (version 3.2) -- */
190 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
191 /* November 2006 */
192
193 /* .. Scalar Arguments .. */
194 /* .. */
195
196 /* Purpose */
197 /* ======= */
198
199 /* DLAMC1 determines the machine parameters given by BETA, T, RND, and */
200 /* IEEE1. */
201
202 /* Arguments */
203 /* ========= */
204
205 /* BETA (output) INTEGER */
206 /* The base of the machine. */
207
208 /* T (output) INTEGER */
209 /* The number of ( BETA ) digits in the mantissa. */
210
211 /* RND (output) LOGICAL */
212 /* Specifies whether proper rounding ( RND = .TRUE. ) or */
213 /* chopping ( RND = .FALSE. ) occurs in addition. This may not */
214 /* be a reliable guide to the way in which the machine performs */
215 /* its arithmetic. */
216
217 /* IEEE1 (output) LOGICAL */
218 /* Specifies whether rounding appears to be done in the IEEE */
219 /* 'round to nearest' style. */
220
221 /* Further Details */
222 /* =============== */
223
224 /* The routine is based on the routine ENVRON by Malcolm and */
225 /* incorporates suggestions by Gentleman and Marovich. See */
226
227 /* Malcolm M. A. (1972) Algorithms to reveal properties of */
228 /* floating-point arithmetic. Comms. of the ACM, 15, 949-951. */
229
230 /* Gentleman W. M. and Marovich S. B. (1974) More on algorithms */
231 /* that reveal properties of floating point arithmetic units. */
232 /* Comms. of the ACM, 17, 276-277. */
233
234 /* ===================================================================== */
235
236 /* .. Local Scalars .. */
237 /* .. */
238 /* .. External Functions .. */
239 /* .. */
240 /* .. Save statement .. */
241 /* .. */
242 /* .. Data statements .. */
243 /* .. */
244 /* .. Executable Statements .. */
245
246 if (first) {
247 one = 1.;
248
249 /* LBETA, LIEEE1, LT and LRND are the local values of BETA, */
250 /* IEEE1, T and RND. */
251
252 /* Throughout this routine we use the function DLAMC3 to ensure */
253 /* that relevant values are stored and not held in registers, or */
254 /* are not affected by optimizers. */
255
256 /* Compute a = 2.0**m with the smallest positive integer m such */
257 /* that */
258
259 /* fl( a + 1.0 ) = a. */
260
261 a = 1.;
262 c__ = 1.;
263
264 /* + WHILE( C.EQ.ONE )LOOP */
265 L10:
266 if (c__ == one) {
267 a *= 2;
268 c__ = dlamc3_(&a, &one);
269 d__1 = -a;
270 c__ = dlamc3_(&c__, &d__1);
271 goto L10;
272 }
273 /* + END WHILE */
274
275 /* Now compute b = 2.0**m with the smallest positive integer m */
276 /* such that */
277
278 /* fl( a + b ) .gt. a. */
279
280 b = 1.;
281 c__ = dlamc3_(&a, &b);
282
283 /* + WHILE( C.EQ.A )LOOP */
284 L20:
285 if (c__ == a) {
286 b *= 2;
287 c__ = dlamc3_(&a, &b);
288 goto L20;
289 }
290 /* + END WHILE */
291
292 /* Now compute the base. a and c are neighbouring floating point */
293 /* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so */
294 /* their difference is beta. Adding 0.25 to c is to ensure that it */
295 /* is truncated to beta and not ( beta - 1 ). */
296
297 qtr = one / 4;
298 savec = c__;
299 d__1 = -a;
300 c__ = dlamc3_(&c__, &d__1);
301 lbeta = (integer) (c__ + qtr);
302
303 /* Now determine whether rounding or chopping occurs, by adding a */
304 /* bit less than beta/2 and a bit more than beta/2 to a. */
305
306 b = (doublereal) lbeta;
307 d__1 = b / 2;
308 d__2 = -b / 100;
309 f = dlamc3_(&d__1, &d__2);
310 c__ = dlamc3_(&f, &a);
311 if (c__ == a) {
312 lrnd = TRUE_;
313 } else {
314 lrnd = FALSE_;
315 }
316 d__1 = b / 2;
317 d__2 = b / 100;
318 f = dlamc3_(&d__1, &d__2);
319 c__ = dlamc3_(&f, &a);
320 if (lrnd && c__ == a) {
321 lrnd = FALSE_;
322 }
323
324 /* Try and decide whether rounding is done in the IEEE 'round to */
325 /* nearest' style. B/2 is half a unit in the last place of the two */
326 /* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */
327 /* zero, and SAVEC is odd. Thus adding B/2 to A should not change */
328 /* A, but adding B/2 to SAVEC should change SAVEC. */
329
330 d__1 = b / 2;
331 t1 = dlamc3_(&d__1, &a);
332 d__1 = b / 2;
333 t2 = dlamc3_(&d__1, &savec);
334 lieee1 = t1 == a && t2 > savec && lrnd;
335
336 /* Now find the mantissa, t. It should be the integer part of */
337 /* log to the base beta of a, however it is safer to determine t */
338 /* by powering. So we find t as the smallest positive integer for */
339 /* which */
340
341 /* fl( beta**t + 1.0 ) = 1.0. */
342
343 lt = 0;
344 a = 1.;
345 c__ = 1.;
346
347 /* + WHILE( C.EQ.ONE )LOOP */
348 L30:
349 if (c__ == one) {
350 ++lt;
351 a *= lbeta;
352 c__ = dlamc3_(&a, &one);
353 d__1 = -a;
354 c__ = dlamc3_(&c__, &d__1);
355 goto L30;
356 }
357 /* + END WHILE */
358
359 }
360
361 *beta = lbeta;
362 *t = lt;
363 *rnd = lrnd;
364 *ieee1 = lieee1;
365 first = FALSE_;
366 return 0;
367
368 /* End of DLAMC1 */
369
370 } /* dlamc1_ */
371
372
373 /* *********************************************************************** */
374
375 /* Subroutine */ int dlamc2_(integer *beta, integer *t, logical *rnd,
376 doublereal *eps, integer *emin, doublereal *rmin, integer *emax,
377 doublereal *rmax)
378 {
379 /* Initialized data */
380
381 static logical first = TRUE_;
382 static logical iwarn = FALSE_;
383
384 /* Format strings */
385 static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre"
386 "ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va"
387 "lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
388 " the IF block as marked within the code of routine\002,\002 DLAM"
389 "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
390
391 /* System generated locals */
392 integer i__1;
393 doublereal d__1, d__2, d__3, d__4, d__5;
394
395 /* Builtin functions */
396 double pow_di(doublereal *, integer *);
397 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
398
399 /* Local variables */
400 doublereal a, b, c__;
401 integer i__;
402 static integer lt;
403 doublereal one, two;
404 logical ieee;
405 doublereal half;
406 logical lrnd;
407 static doublereal leps;
408 doublereal zero;
409 static integer lbeta;
410 doublereal rbase;
411 static integer lemin, lemax;
412 integer gnmin;
413 doublereal small;
414 integer gpmin;
415 doublereal third;
416 static doublereal lrmin, lrmax;
417 doublereal sixth;
418 extern /* Subroutine */ int dlamc1_(integer *, integer *, logical *,
419 logical *);
420 extern doublereal dlamc3_(doublereal *, doublereal *);
421 logical lieee1;
422 extern /* Subroutine */ int dlamc4_(integer *, doublereal *, integer *),
423 dlamc5_(integer *, integer *, integer *, logical *, integer *,
424 doublereal *);
425 integer ngnmin, ngpmin;
426
427 /* Fortran I/O blocks */
428 static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
429
430
431
432 /* -- LAPACK auxiliary routine (version 3.2) -- */
433 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
434 /* November 2006 */
435
436 /* .. Scalar Arguments .. */
437 /* .. */
438
439 /* Purpose */
440 /* ======= */
441
442 /* DLAMC2 determines the machine parameters specified in its argument */
443 /* list. */
444
445 /* Arguments */
446 /* ========= */
447
448 /* BETA (output) INTEGER */
449 /* The base of the machine. */
450
451 /* T (output) INTEGER */
452 /* The number of ( BETA ) digits in the mantissa. */
453
454 /* RND (output) LOGICAL */
455 /* Specifies whether proper rounding ( RND = .TRUE. ) or */
456 /* chopping ( RND = .FALSE. ) occurs in addition. This may not */
457 /* be a reliable guide to the way in which the machine performs */
458 /* its arithmetic. */
459
460 /* EPS (output) DOUBLE PRECISION */
461 /* The smallest positive number such that */
462
463 /* fl( 1.0 - EPS ) .LT. 1.0, */
464
465 /* where fl denotes the computed value. */
466
467 /* EMIN (output) INTEGER */
468 /* The minimum exponent before (gradual) underflow occurs. */
469
470 /* RMIN (output) DOUBLE PRECISION */
471 /* The smallest normalized number for the machine, given by */
472 /* BASE**( EMIN - 1 ), where BASE is the floating point value */
473 /* of BETA. */
474
475 /* EMAX (output) INTEGER */
476 /* The maximum exponent before overflow occurs. */
477
478 /* RMAX (output) DOUBLE PRECISION */
479 /* The largest positive number for the machine, given by */
480 /* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
481 /* value of BETA. */
482
483 /* Further Details */
484 /* =============== */
485
486 /* The computation of EPS is based on a routine PARANOIA by */
487 /* W. Kahan of the University of California at Berkeley. */
488
489 /* ===================================================================== */
490
491 /* .. Local Scalars .. */
492 /* .. */
493 /* .. External Functions .. */
494 /* .. */
495 /* .. External Subroutines .. */
496 /* .. */
497 /* .. Intrinsic Functions .. */
498 /* .. */
499 /* .. Save statement .. */
500 /* .. */
501 /* .. Data statements .. */
502 /* .. */
503 /* .. Executable Statements .. */
504
505 if (first) {
506 zero = 0.;
507 one = 1.;
508 two = 2.;
509
510 /* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */
511 /* BETA, T, RND, EPS, EMIN and RMIN. */
512
513 /* Throughout this routine we use the function DLAMC3 to ensure */
514 /* that relevant values are stored and not held in registers, or */
515 /* are not affected by optimizers. */
516
517 /* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */
518
519 dlamc1_(&lbeta, &lt, &lrnd, &lieee1);
520
521 /* Start to find EPS. */
522
523 b = (doublereal) lbeta;
524 i__1 = -lt;
525 a = pow_di(&b, &i__1);
526 leps = a;
527
528 /* Try some tricks to see whether or not this is the correct EPS. */
529
530 b = two / 3;
531 half = one / 2;
532 d__1 = -half;
533 sixth = dlamc3_(&b, &d__1);
534 third = dlamc3_(&sixth, &sixth);
535 d__1 = -half;
536 b = dlamc3_(&third, &d__1);
537 b = dlamc3_(&b, &sixth);
538 b = abs(b);
539 if (b < leps) {
540 b = leps;
541 }
542
543 leps = 1.;
544
545 /* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
546 L10:
547 if (leps > b && b > zero) {
548 leps = b;
549 d__1 = half * leps;
550 /* Computing 5th power */
551 d__3 = two, d__4 = d__3, d__3 *= d__3;
552 /* Computing 2nd power */
553 d__5 = leps;
554 d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
555 c__ = dlamc3_(&d__1, &d__2);
556 d__1 = -c__;
557 c__ = dlamc3_(&half, &d__1);
558 b = dlamc3_(&half, &c__);
559 d__1 = -b;
560 c__ = dlamc3_(&half, &d__1);
561 b = dlamc3_(&half, &c__);
562 goto L10;
563 }
564 /* + END WHILE */
565
566 if (a < leps) {
567 leps = a;
568 }
569
570 /* Computation of EPS complete. */
571
572 /* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */
573 /* Keep dividing A by BETA until (gradual) underflow occurs. This */
574 /* is detected when we cannot recover the previous A. */
575
576 rbase = one / lbeta;
577 small = one;
578 for (i__ = 1; i__ <= 3; ++i__) {
579 d__1 = small * rbase;
580 small = dlamc3_(&d__1, &zero);
581 /* L20: */
582 }
583 a = dlamc3_(&one, &small);
584 dlamc4_(&ngpmin, &one, &lbeta);
585 d__1 = -one;
586 dlamc4_(&ngnmin, &d__1, &lbeta);
587 dlamc4_(&gpmin, &a, &lbeta);
588 d__1 = -a;
589 dlamc4_(&gnmin, &d__1, &lbeta);
590 ieee = FALSE_;
591
592 if (ngpmin == ngnmin && gpmin == gnmin) {
593 if (ngpmin == gpmin) {
594 lemin = ngpmin;
595 /* ( Non twos-complement machines, no gradual underflow; */
596 /* e.g., VAX ) */
597 } else if (gpmin - ngpmin == 3) {
598 lemin = ngpmin - 1 + lt;
599 ieee = TRUE_;
600 /* ( Non twos-complement machines, with gradual underflow; */
601 /* e.g., IEEE standard followers ) */
602 } else {
603 lemin = min(ngpmin,gpmin);
604 /* ( A guess; no known machine ) */
605 iwarn = TRUE_;
606 }
607
608 } else if (ngpmin == gpmin && ngnmin == gnmin) {
609 if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
610 lemin = max(ngpmin,ngnmin);
611 /* ( Twos-complement machines, no gradual underflow; */
612 /* e.g., CYBER 205 ) */
613 } else {
614 lemin = min(ngpmin,ngnmin);
615 /* ( A guess; no known machine ) */
616 iwarn = TRUE_;
617 }
618
619 } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
620 {
621 if (gpmin - min(ngpmin,ngnmin) == 3) {
622 lemin = max(ngpmin,ngnmin) - 1 + lt;
623 /* ( Twos-complement machines with gradual underflow; */
624 /* no known machine ) */
625 } else {
626 lemin = min(ngpmin,ngnmin);
627 /* ( A guess; no known machine ) */
628 iwarn = TRUE_;
629 }
630
631 } else {
632 /* Computing MIN */
633 i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin);
634 lemin = min(i__1,gnmin);
635 /* ( A guess; no known machine ) */
636 iwarn = TRUE_;
637 }
638 first = FALSE_;
639 /* ** */
640 /* Comment out this if block if EMIN is ok */
641 if (iwarn) {
642 first = TRUE_;
643 s_wsfe(&io___58);
644 do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer));
645 e_wsfe();
646 }
647 /* ** */
648
649 /* Assume IEEE arithmetic if we found denormalised numbers above, */
650 /* or if arithmetic seems to round in the IEEE style, determined */
651 /* in routine DLAMC1. A true IEEE machine should have both things */
652 /* true; however, faulty machines may have one or the other. */
653
654 ieee = ieee || lieee1;
655
656 /* Compute RMIN by successive division by BETA. We could compute */
657 /* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */
658 /* this computation. */
659
660 lrmin = 1.;
661 i__1 = 1 - lemin;
662 for (i__ = 1; i__ <= i__1; ++i__) {
663 d__1 = lrmin * rbase;
664 lrmin = dlamc3_(&d__1, &zero);
665 /* L30: */
666 }
667
668 /* Finally, call DLAMC5 to compute EMAX and RMAX. */
669
670 dlamc5_(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);
671 }
672
673 *beta = lbeta;
674 *t = lt;
675 *rnd = lrnd;
676 *eps = leps;
677 *emin = lemin;
678 *rmin = lrmin;
679 *emax = lemax;
680 *rmax = lrmax;
681
682 return 0;
683
684
685 /* End of DLAMC2 */
686
687 } /* dlamc2_ */
688
689
690 /* *********************************************************************** */
691
692 doublereal dlamc3_(doublereal *a, doublereal *b)
693 {
694 /* System generated locals */
695 doublereal ret_val;
696
697
698 /* -- LAPACK auxiliary routine (version 3.2) -- */
699 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
700 /* November 2006 */
701
702 /* .. Scalar Arguments .. */
703 /* .. */
704
705 /* Purpose */
706 /* ======= */
707
708 /* DLAMC3 is intended to force A and B to be stored prior to doing */
709 /* the addition of A and B , for use in situations where optimizers */
710 /* might hold one of these in a register. */
711
712 /* Arguments */
713 /* ========= */
714
715 /* A (input) DOUBLE PRECISION */
716 /* B (input) DOUBLE PRECISION */
717 /* The values A and B. */
718
719 /* ===================================================================== */
720
721 /* .. Executable Statements .. */
722
723 ret_val = *a + *b;
724
725 return ret_val;
726
727 /* End of DLAMC3 */
728
729 } /* dlamc3_ */
730
731
732 /* *********************************************************************** */
733
734 /* Subroutine */ int dlamc4_(integer *emin, doublereal *start, integer *base)
735 {
736 /* System generated locals */
737 integer i__1;
738 doublereal d__1;
739
740 /* Local variables */
741 doublereal a;
742 integer i__;
743 doublereal b1, b2, c1, c2, d1, d2, one, zero, rbase;
744 extern doublereal dlamc3_(doublereal *, doublereal *);
745
746
747 /* -- LAPACK auxiliary routine (version 3.2) -- */
748 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
749 /* November 2006 */
750
751 /* .. Scalar Arguments .. */
752 /* .. */
753
754 /* Purpose */
755 /* ======= */
756
757 /* DLAMC4 is a service routine for DLAMC2. */
758
759 /* Arguments */
760 /* ========= */
761
762 /* EMIN (output) INTEGER */
763 /* The minimum exponent before (gradual) underflow, computed by */
764 /* setting A = START and dividing by BASE until the previous A */
765 /* can not be recovered. */
766
767 /* START (input) DOUBLE PRECISION */
768 /* The starting point for determining EMIN. */
769
770 /* BASE (input) INTEGER */
771 /* The base of the machine. */
772
773 /* ===================================================================== */
774
775 /* .. Local Scalars .. */
776 /* .. */
777 /* .. External Functions .. */
778 /* .. */
779 /* .. Executable Statements .. */
780
781 a = *start;
782 one = 1.;
783 rbase = one / *base;
784 zero = 0.;
785 *emin = 1;
786 d__1 = a * rbase;
787 b1 = dlamc3_(&d__1, &zero);
788 c1 = a;
789 c2 = a;
790 d1 = a;
791 d2 = a;
792 /* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
793 /* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
794 L10:
795 if (c1 == a && c2 == a && d1 == a && d2 == a) {
796 --(*emin);
797 a = b1;
798 d__1 = a / *base;
799 b1 = dlamc3_(&d__1, &zero);
800 d__1 = b1 * *base;
801 c1 = dlamc3_(&d__1, &zero);
802 d1 = zero;
803 i__1 = *base;
804 for (i__ = 1; i__ <= i__1; ++i__) {
805 d1 += b1;
806 /* L20: */
807 }
808 d__1 = a * rbase;
809 b2 = dlamc3_(&d__1, &zero);
810 d__1 = b2 / rbase;
811 c2 = dlamc3_(&d__1, &zero);
812 d2 = zero;
813 i__1 = *base;
814 for (i__ = 1; i__ <= i__1; ++i__) {
815 d2 += b2;
816 /* L30: */
817 }
818 goto L10;
819 }
820 /* + END WHILE */
821
822 return 0;
823
824 /* End of DLAMC4 */
825
826 } /* dlamc4_ */
827
828
829 /* *********************************************************************** */
830
831 /* Subroutine */ int dlamc5_(integer *beta, integer *p, integer *emin,
832 logical *ieee, integer *emax, doublereal *rmax)
833 {
834 /* System generated locals */
835 integer i__1;
836 doublereal d__1;
837
838 /* Local variables */
839 integer i__;
840 doublereal y, z__;
841 integer try__, lexp;
842 doublereal oldy;
843 integer uexp, nbits;
844 extern doublereal dlamc3_(doublereal *, doublereal *);
845 doublereal recbas;
846 integer exbits, expsum;
847
848
849 /* -- LAPACK auxiliary routine (version 3.2) -- */
850 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
851 /* November 2006 */
852
853 /* .. Scalar Arguments .. */
854 /* .. */
855
856 /* Purpose */
857 /* ======= */
858
859 /* DLAMC5 attempts to compute RMAX, the largest machine floating-point */
860 /* number, without overflow. It assumes that EMAX + abs(EMIN) sum */
861 /* approximately to a power of 2. It will fail on machines where this */
862 /* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
863 /* EMAX = 28718). It will also fail if the value supplied for EMIN is */
864 /* too large (i.e. too close to zero), probably with overflow. */
865
866 /* Arguments */
867 /* ========= */
868
869 /* BETA (input) INTEGER */
870 /* The base of floating-point arithmetic. */
871
872 /* P (input) INTEGER */
873 /* The number of base BETA digits in the mantissa of a */
874 /* floating-point value. */
875
876 /* EMIN (input) INTEGER */
877 /* The minimum exponent before (gradual) underflow. */
878
879 /* IEEE (input) LOGICAL */
880 /* A logical flag specifying whether or not the arithmetic */
881 /* system is thought to comply with the IEEE standard. */
882
883 /* EMAX (output) INTEGER */
884 /* The largest exponent before overflow */
885
886 /* RMAX (output) DOUBLE PRECISION */
887 /* The largest machine floating-point number. */
888
889 /* ===================================================================== */
890
891 /* .. Parameters .. */
892 /* .. */
893 /* .. Local Scalars .. */
894 /* .. */
895 /* .. External Functions .. */
896 /* .. */
897 /* .. Intrinsic Functions .. */
898 /* .. */
899 /* .. Executable Statements .. */
900
901 /* First compute LEXP and UEXP, two powers of 2 that bound */
902 /* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
903 /* approximately to the bound that is closest to abs(EMIN). */
904 /* (EMAX is the exponent of the required number RMAX). */
905
906 lexp = 1;
907 exbits = 1;
908 L10:
909 try__ = lexp << 1;
910 if (try__ <= -(*emin)) {
911 lexp = try__;
912 ++exbits;
913 goto L10;
914 }
915 if (lexp == -(*emin)) {
916 uexp = lexp;
917 } else {
918 uexp = try__;
919 ++exbits;
920 }
921
922 /* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
923 /* than or equal to EMIN. EXBITS is the number of bits needed to */
924 /* store the exponent. */
925
926 if (uexp + *emin > -lexp - *emin) {
927 expsum = lexp << 1;
928 } else {
929 expsum = uexp << 1;
930 }
931
932 /* EXPSUM is the exponent range, approximately equal to */
933 /* EMAX - EMIN + 1 . */
934
935 *emax = expsum + *emin - 1;
936 nbits = exbits + 1 + *p;
937
938 /* NBITS is the total number of bits needed to store a */
939 /* floating-point number. */
940
941 if (nbits % 2 == 1 && *beta == 2) {
942
943 /* Either there are an odd number of bits used to store a */
944 /* floating-point number, which is unlikely, or some bits are */
945 /* not used in the representation of numbers, which is possible, */
946 /* (e.g. Cray machines) or the mantissa has an implicit bit, */
947 /* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
948 /* most likely. We have to assume the last alternative. */
949 /* If this is true, then we need to reduce EMAX by one because */
950 /* there must be some way of representing zero in an implicit-bit */
951 /* system. On machines like Cray, we are reducing EMAX by one */
952 /* unnecessarily. */
953
954 --(*emax);
955 }
956
957 if (*ieee) {
958
959 /* Assume we are on an IEEE machine which reserves one exponent */
960 /* for infinity and NaN. */
961
962 --(*emax);
963 }
964
965 /* Now create RMAX, the largest machine number, which should */
966 /* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
967
968 /* First compute 1.0 - BETA**(-P), being careful that the */
969 /* result is less than 1.0 . */
970
971 recbas = 1. / *beta;
972 z__ = *beta - 1.;
973 y = 0.;
974 i__1 = *p;
975 for (i__ = 1; i__ <= i__1; ++i__) {
976 z__ *= recbas;
977 if (y < 1.) {
978 oldy = y;
979 }
980 y = dlamc3_(&y, &z__);
981 /* L20: */
982 }
983 if (y >= 1.) {
984 y = oldy;
985 }
986
987 /* Now multiply by BETA**EMAX to get RMAX. */
988
989 i__1 = *emax;
990 for (i__ = 1; i__ <= i__1; ++i__) {
991 d__1 = y * *beta;
992 y = dlamc3_(&d__1, &c_b32);
993 /* L30: */
994 }
995
996 *rmax = y;
997 return 0;
998
999 /* End of DLAMC5 */
1000
1001 } /* dlamc5_ */