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