libflame revision_anchor
Functions
fla_slamch.c File Reference

(r)

Functions

double fla_pow_ri (real *ap, integer *bp)
 
real fla_slamch (char *cmach, ftnlen cmach_len)
 
int fla_slamc1 (integer *beta, integer *t, logical *rnd, logical *ieee1)
 
int fla_slamc2 (integer *beta, integer *t, logical *rnd, real *eps, integer *emin, real *rmin, integer *emax, real *rmax)
 
real fla_slamc3 (real *a, real *b)
 
int fla_slamc4 (integer *emin, real *start, integer *base)
 
int fla_slamc5 (integer *beta, integer *p, integer *emin, logical *ieee, integer *emax, real *rmax)
 

Function Documentation

◆ fla_pow_ri()

double fla_pow_ri ( real *  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}
integer
int integer
Definition FLA_f2c.h:25
i
int i
Definition bl1_axmyv2.c:145

References i.

Referenced by fla_slamc2(), and fla_slamch().

◆ fla_slamc1()

int fla_slamc1 ( integer *  beta,
integer *  t,
logical *  rnd,
logical *  ieee1 
)
203{
204 /* Initialized data */
205
206 static logical first = TRUE_;
207
208 /* System generated locals */
209 real r__1, r__2;
210
211 /* Local variables */
212 static logical lrnd;
213 static real a, b, c__, f;
214 static integer lbeta;
215 static real savec;
216 static logical lieee1;
217 static real t1, t2;
218 extern real fla_slamc3(real *, real *);
219 static integer lt;
220 static real one, qtr;
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/* SLAMC1 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 = (float)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 SLAMC3 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 = (float)1.;
296 c__ = (float)1.;
297
298/* + WHILE( C.EQ.ONE )LOOP */
299L10:
300 if (c__ == one) {
301 a *= 2;
302 c__ = fla_slamc3(&a, &one);
303 r__1 = -a;
304 c__ = fla_slamc3(&c__, &r__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 = (float)1.;
315 c__ = fla_slamc3(&a, &b);
316
317/* + WHILE( C.EQ.A )LOOP */
318L20:
319 if (c__ == a) {
320 b *= 2;
321 c__ = fla_slamc3(&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 r__1 = -a;
334 c__ = fla_slamc3(&c__, &r__1);
335 lbeta = c__ + qtr;
336
337/* Now determine whether rounding or chopping occurs, by adding a */
338/* bit less than beta/2 and a bit more than beta/2 to a. */
339
340 b = (real) lbeta;
341 r__1 = b / 2;
342 r__2 = -b / 100;
343 f = fla_slamc3(&r__1, &r__2);
344 c__ = fla_slamc3(&f, &a);
345 if (c__ == a) {
346 lrnd = TRUE_;
347 } else {
348 lrnd = FALSE_;
349 }
350 r__1 = b / 2;
351 r__2 = b / 100;
352 f = fla_slamc3(&r__1, &r__2);
353 c__ = fla_slamc3(&f, &a);
354 if (lrnd && c__ == a) {
355 lrnd = FALSE_;
356 }
357
358/* Try and decide whether rounding is done in the IEEE 'round to */
359/* nearest' style. B/2 is half a unit in the last place of the two */
360/* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */
361/* zero, and SAVEC is odd. Thus adding B/2 to A should not change */
362/* A, but adding B/2 to SAVEC should change SAVEC. */
363
364 r__1 = b / 2;
365 t1 = fla_slamc3(&r__1, &a);
366 r__1 = b / 2;
367 t2 = fla_slamc3(&r__1, &savec);
368 lieee1 = t1 == a && t2 > savec && lrnd;
369
370/* Now find the mantissa, t. It should be the integer part of */
371/* log to the base beta of a, however it is safer to determine t */
372/* by powering. So we find t as the smallest positive integer for */
373/* which */
374
375/* fl( beta**t + 1.0 ) = 1.0. */
376
377 lt = 0;
378 a = (float)1.;
379 c__ = (float)1.;
380
381/* + WHILE( C.EQ.ONE )LOOP */
382L30:
383 if (c__ == one) {
384 ++lt;
385 a *= lbeta;
386 c__ = fla_slamc3(&a, &one);
387 r__1 = -a;
388 c__ = fla_slamc3(&c__, &r__1);
389 goto L30;
390 }
391/* + END WHILE */
392
393 }
394
395 *beta = lbeta;
396 *t = lt;
397 *rnd = lrnd;
398 *ieee1 = lieee1;
399 first = FALSE_;
400 return 0;
401
402/* End of SLAMC1 */
403
404} /* fla_slamc1_ */
logical
int logical
Definition FLA_f2c.h:36
real
float real
Definition FLA_f2c.h:30
fla_slamc3
real fla_slamc3(real *a, real *b)
Definition fla_slamch.c:721

References fla_slamc3(), and i.

Referenced by fla_slamc2().

◆ fla_slamc2()

int fla_slamc2 ( integer *  beta,
integer *  t,
logical *  rnd,
real *  eps,
integer *  emin,
real *  rmin,
integer *  emax,
real *  rmax 
)
411{
412 /* Initialized data */
413
414 static logical first = TRUE_;
415 static logical iwarn = FALSE_;
416
417 /* Format strings */
418 static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre\
419ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the value EMIN loo\
420ks\002,\002 acceptable please comment out \002,/\002 the IF block as marked \
421within the code of routine\002,\002 SLAMC2,\002,/\002 otherwise supply EMIN \
422explicitly.\002,/)";
423
424 /* System generated locals */
425 integer i__1;
426 real r__1, r__2, r__3, r__4, r__5;
427
428 /* Builtin functions */
429 double fla_pow_ri(real *, integer *);
430 //integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe();
431
432 /* Local variables */
433 static logical ieee;
434 static real half;
435 static logical lrnd;
436 static real leps, zero, a, b, c__;
437 static integer i__, lbeta;
438 static real rbase;
439 static integer lemin, lemax, gnmin;
440 static real small;
441 static integer gpmin;
442 static real third, lrmin, lrmax, sixth;
443 static logical lieee1;
444 extern /* Subroutine */ int fla_slamc1(integer *, integer *, logical *,
445 logical *);
446 extern real fla_slamc3(real *, real *);
447 extern /* Subroutine */ int fla_slamc4(integer *, real *, integer *),
448 fla_slamc5(integer *, integer *, integer *, logical *, integer *,
449 real *);
450 static integer lt, ngnmin, ngpmin;
451 static real one, two;
452
453 /* Fortran I/O blocks */
454 //static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
455
456
457
458/* -- LAPACK auxiliary routine (version 3.2) -- */
459/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
460/* November 2006 */
461
462/* .. Scalar Arguments .. */
463/* .. */
464
465/* Purpose */
466/* ======= */
467
468/* SLAMC2 determines the machine parameters specified in its argument */
469/* list. */
470
471/* Arguments */
472/* ========= */
473
474/* BETA (output) INTEGER */
475/* The base of the machine. */
476
477/* T (output) INTEGER */
478/* The number of ( BETA ) digits in the mantissa. */
479
480/* RND (output) LOGICAL */
481/* Specifies whether proper rounding ( RND = .TRUE. ) or */
482/* chopping ( RND = .FALSE. ) occurs in addition. This may not */
483/* be a reliable guide to the way in which the machine performs */
484/* its arithmetic. */
485
486/* EPS (output) REAL */
487/* The smallest positive number such that */
488
489/* fl( 1.0 - EPS ) .LT. 1.0, */
490
491/* where fl denotes the computed value. */
492
493/* EMIN (output) INTEGER */
494/* The minimum exponent before (gradual) underflow occurs. */
495
496/* RMIN (output) REAL */
497/* The smallest normalized number for the machine, given by */
498/* BASE**( EMIN - 1 ), where BASE is the floating point value */
499/* of BETA. */
500
501/* EMAX (output) INTEGER */
502/* The maximum exponent before overflow occurs. */
503
504/* RMAX (output) REAL */
505/* The largest positive number for the machine, given by */
506/* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
507/* value of BETA. */
508
509/* Further Details */
510/* =============== */
511
512/* The computation of EPS is based on a routine PARANOIA by */
513/* W. Kahan of the University of California at Berkeley. */
514
515/* ===================================================================== */
516
517/* .. Local Scalars .. */
518/* .. */
519/* .. External Functions .. */
520/* .. */
521/* .. External Subroutines .. */
522/* .. */
523/* .. Intrinsic Functions .. */
524/* .. */
525/* .. Save statement .. */
526/* .. */
527/* .. Data statements .. */
528/* .. */
529/* .. Executable Statements .. */
530
531 if (first) {
532 zero = (float)0.;
533 one = (float)1.;
534 two = (float)2.;
535
536/* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */
537/* BETA, T, RND, EPS, EMIN and RMIN. */
538
539/* Throughout this routine we use the function SLAMC3 to ensure */
540/* that relevant values are stored and not held in registers, or */
541/* are not affected by optimizers. */
542
543/* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */
544
545 fla_slamc1(&lbeta, &lt, &lrnd, &lieee1);
546
547/* Start to find EPS. */
548
549 b = (real) lbeta;
550 i__1 = -lt;
551 a = fla_pow_ri(&b, &i__1);
552 leps = a;
553
554/* Try some tricks to see whether or not this is the correct EPS. */
555
556 b = two / 3;
557 half = one / 2;
558 r__1 = -half;
559 sixth = fla_slamc3(&b, &r__1);
560 third = fla_slamc3(&sixth, &sixth);
561 r__1 = -half;
562 b = fla_slamc3(&third, &r__1);
563 b = fla_slamc3(&b, &sixth);
564 b = abs(b);
565 if (b < leps) {
566 b = leps;
567 }
568
569 leps = (float)1.;
570
571/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
572L10:
573 if (leps > b && b > zero) {
574 leps = b;
575 r__1 = half * leps;
576/* Computing 5th power */
577 r__3 = two, r__4 = r__3, r__3 *= r__3;
578/* Computing 2nd power */
579 r__5 = leps;
580 r__2 = r__4 * (r__3 * r__3) * (r__5 * r__5);
581 c__ = fla_slamc3(&r__1, &r__2);
582 r__1 = -c__;
583 c__ = fla_slamc3(&half, &r__1);
584 b = fla_slamc3(&half, &c__);
585 r__1 = -b;
586 c__ = fla_slamc3(&half, &r__1);
587 b = fla_slamc3(&half, &c__);
588 goto L10;
589 }
590/* + END WHILE */
591
592 if (a < leps) {
593 leps = a;
594 }
595
596/* Computation of EPS complete. */
597
598/* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */
599/* Keep dividing A by BETA until (gradual) underflow occurs. This */
600/* is detected when we cannot recover the previous A. */
601
602 rbase = one / lbeta;
603 small = one;
604 for (i__ = 1; i__ <= 3; ++i__) {
605 r__1 = small * rbase;
606 small = fla_slamc3(&r__1, &zero);
607/* L20: */
608 }
609 a = fla_slamc3(&one, &small);
610 fla_slamc4(&ngpmin, &one, &lbeta);
611 r__1 = -one;
612 fla_slamc4(&ngnmin, &r__1, &lbeta);
613 fla_slamc4(&gpmin, &a, &lbeta);
614 r__1 = -a;
615 fla_slamc4(&gnmin, &r__1, &lbeta);
616 ieee = FALSE_;
617
618 if (ngpmin == ngnmin && gpmin == gnmin) {
619 if (ngpmin == gpmin) {
620 lemin = ngpmin;
621/* ( Non twos-complement machines, no gradual underflow; */
622/* e.g., VAX ) */
623 } else if (gpmin - ngpmin == 3) {
624 lemin = ngpmin - 1 + lt;
625 ieee = TRUE_;
626/* ( Non twos-complement machines, with gradual underflow; */
627/* e.g., IEEE standard followers ) */
628 } else {
629 lemin = min(ngpmin,gpmin);
630/* ( A guess; no known machine ) */
631 iwarn = TRUE_;
632 }
633
634 } else if (ngpmin == gpmin && ngnmin == gnmin) {
635 if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
636 lemin = max(ngpmin,ngnmin);
637/* ( Twos-complement machines, no gradual underflow; */
638/* e.g., CYBER 205 ) */
639 } else {
640 lemin = min(ngpmin,ngnmin);
641/* ( A guess; no known machine ) */
642 iwarn = TRUE_;
643 }
644
645 } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
646 {
647 if (gpmin - min(ngpmin,ngnmin) == 3) {
648 lemin = max(ngpmin,ngnmin) - 1 + lt;
649/* ( Twos-complement machines with gradual underflow; */
650/* no known machine ) */
651 } else {
652 lemin = min(ngpmin,ngnmin);
653/* ( A guess; no known machine ) */
654 iwarn = TRUE_;
655 }
656
657 } else {
658/* Computing MIN */
659 i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin);
660 lemin = min(i__1,gnmin);
661/* ( A guess; no known machine ) */
662 iwarn = TRUE_;
663 }
664 first = FALSE_;
665/* ** */
666/* Comment out this if block if EMIN is ok */
667 if (iwarn) {
668 first = TRUE_;
669/*
670 s_wsfe(&io___58);
671 do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer));
672 e_wsfe();
673*/
674 printf( "%s", fmt_9999 );
675 }
676/* ** */
677
678/* Assume IEEE arithmetic if we found denormalised numbers above, */
679/* or if arithmetic seems to round in the IEEE style, determined */
680/* in routine SLAMC1. A true IEEE machine should have both things */
681/* true; however, faulty machines may have one or the other. */
682
683 ieee = ieee || lieee1;
684
685/* Compute RMIN by successive division by BETA. We could compute */
686/* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */
687/* this computation. */
688
689 lrmin = (float)1.;
690 i__1 = 1 - lemin;
691 for (i__ = 1; i__ <= i__1; ++i__) {
692 r__1 = lrmin * rbase;
693 lrmin = fla_slamc3(&r__1, &zero);
694/* L30: */
695 }
696
697/* Finally, call SLAMC5 to compute EMAX and RMAX. */
698
699 fla_slamc5(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);
700 }
701
702 *beta = lbeta;
703 *t = lt;
704 *rnd = lrnd;
705 *eps = leps;
706 *emin = lemin;
707 *rmin = lrmin;
708 *emax = lemax;
709 *rmax = lrmax;
710
711 return 0;
712
713
714/* End of SLAMC2 */
715
716} /* fla_slamc2_ */
fla_slamc5
int fla_slamc5(integer *beta, integer *p, integer *emin, logical *ieee, integer *emax, real *rmax)
Definition fla_slamch.c:861
fla_slamc4
int fla_slamc4(integer *emin, real *start, integer *base)
Definition fla_slamch.c:763
fla_pow_ri
double fla_pow_ri(real *ap, integer *bp)
Definition fla_slamch.c:26
fla_slamc1
int fla_slamc1(integer *beta, integer *t, logical *rnd, logical *ieee1)
Definition fla_slamch.c:201

References fla_pow_ri(), fla_slamc1(), fla_slamc3(), fla_slamc4(), fla_slamc5(), and i.

Referenced by fla_slamch().

◆ fla_slamc3()

real fla_slamc3 ( real *  a,
real *  b 
)
722{
723 /* System generated locals */
724 real ret_val;
725
726
727/* -- LAPACK auxiliary routine (version 3.2) -- */
728/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
729/* November 2006 */
730
731/* .. Scalar Arguments .. */
732/* .. */
733
734/* Purpose */
735/* ======= */
736
737/* SLAMC3 is intended to force A and B to be stored prior to doing */
738/* the addition of A and B , for use in situations where optimizers */
739/* might hold one of these in a register. */
740
741/* Arguments */
742/* ========= */
743
744/* A (input) REAL */
745/* B (input) REAL */
746/* The values A and B. */
747
748/* ===================================================================== */
749
750/* .. Executable Statements .. */
751
752 ret_val = *a + *b;
753
754 return ret_val;
755
756/* End of SLAMC3 */
757
758} /* fla_slamc3_ */

References i.

Referenced by fla_slamc1(), fla_slamc2(), fla_slamc4(), and fla_slamc5().

◆ fla_slamc4()

int fla_slamc4 ( integer *  emin,
real *  start,
integer *  base 
)
764{
765 /* System generated locals */
766 integer i__1;
767 real r__1;
768
769 /* Local variables */
770 static real zero, a;
771 static integer i__;
772 static real rbase, b1, b2, c1, c2, d1, d2;
773 extern real fla_slamc3(real *, real *);
774 static real one;
775
776
777/* -- LAPACK auxiliary routine (version 3.2) -- */
778/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
779/* November 2006 */
780
781/* .. Scalar Arguments .. */
782/* .. */
783
784/* Purpose */
785/* ======= */
786
787/* SLAMC4 is a service routine for SLAMC2. */
788
789/* Arguments */
790/* ========= */
791
792/* EMIN (output) INTEGER */
793/* The minimum exponent before (gradual) underflow, computed by */
794/* setting A = START and dividing by BASE until the previous A */
795/* can not be recovered. */
796
797/* START (input) REAL */
798/* The starting point for determining EMIN. */
799
800/* BASE (input) INTEGER */
801/* The base of the machine. */
802
803/* ===================================================================== */
804
805/* .. Local Scalars .. */
806/* .. */
807/* .. External Functions .. */
808/* .. */
809/* .. Executable Statements .. */
810
811 a = *start;
812 one = (float)1.;
813 rbase = one / *base;
814 zero = (float)0.;
815 *emin = 1;
816 r__1 = a * rbase;
817 b1 = fla_slamc3(&r__1, &zero);
818 c1 = a;
819 c2 = a;
820 d1 = a;
821 d2 = a;
822/* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
823/* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
824L10:
825 if (c1 == a && c2 == a && d1 == a && d2 == a) {
826 --(*emin);
827 a = b1;
828 r__1 = a / *base;
829 b1 = fla_slamc3(&r__1, &zero);
830 r__1 = b1 * *base;
831 c1 = fla_slamc3(&r__1, &zero);
832 d1 = zero;
833 i__1 = *base;
834 for (i__ = 1; i__ <= i__1; ++i__) {
835 d1 += b1;
836/* L20: */
837 }
838 r__1 = a * rbase;
839 b2 = fla_slamc3(&r__1, &zero);
840 r__1 = b2 / rbase;
841 c2 = fla_slamc3(&r__1, &zero);
842 d2 = zero;
843 i__1 = *base;
844 for (i__ = 1; i__ <= i__1; ++i__) {
845 d2 += b2;
846/* L30: */
847 }
848 goto L10;
849 }
850/* + END WHILE */
851
852 return 0;
853
854/* End of SLAMC4 */
855
856} /* fla_slamc4_ */

References fla_slamc3(), and i.

Referenced by fla_slamc2().

◆ fla_slamc5()

int fla_slamc5 ( integer *  beta,
integer *  p,
integer *  emin,
logical *  ieee,
integer *  emax,
real *  rmax 
)
863{
864 /* System generated locals */
865 integer i__1;
866 real r__1;
867
868 /* Local variables */
869 static integer lexp;
870 static real oldy;
871 static integer uexp, i__;
872 static real y, z__;
873 static integer nbits;
874 extern real fla_slamc3(real *, real *);
875 static real recbas;
876 static integer exbits, expsum, try__;
877
878
879/* -- LAPACK auxiliary routine (version 3.2) -- */
880/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
881/* November 2006 */
882
883/* .. Scalar Arguments .. */
884/* .. */
885
886/* Purpose */
887/* ======= */
888
889/* SLAMC5 attempts to compute RMAX, the largest machine floating-point */
890/* number, without overflow. It assumes that EMAX + abs(EMIN) sum */
891/* approximately to a power of 2. It will fail on machines where this */
892/* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
893/* EMAX = 28718). It will also fail if the value supplied for EMIN is */
894/* too large (i.e. too close to zero), probably with overflow. */
895
896/* Arguments */
897/* ========= */
898
899/* BETA (input) INTEGER */
900/* The base of floating-point arithmetic. */
901
902/* P (input) INTEGER */
903/* The number of base BETA digits in the mantissa of a */
904/* floating-point value. */
905
906/* EMIN (input) INTEGER */
907/* The minimum exponent before (gradual) underflow. */
908
909/* IEEE (input) LOGICAL */
910/* A logical flag specifying whether or not the arithmetic */
911/* system is thought to comply with the IEEE standard. */
912
913/* EMAX (output) INTEGER */
914/* The largest exponent before overflow */
915
916/* RMAX (output) REAL */
917/* The largest machine floating-point number. */
918
919/* ===================================================================== */
920
921/* .. Parameters .. */
922/* .. */
923/* .. Local Scalars .. */
924/* .. */
925/* .. External Functions .. */
926/* .. */
927/* .. Intrinsic Functions .. */
928/* .. */
929/* .. Executable Statements .. */
930
931/* First compute LEXP and UEXP, two powers of 2 that bound */
932/* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
933/* approximately to the bound that is closest to abs(EMIN). */
934/* (EMAX is the exponent of the required number RMAX). */
935
936 lexp = 1;
937 exbits = 1;
938L10:
939 try__ = lexp << 1;
940 if (try__ <= -(*emin)) {
941 lexp = try__;
942 ++exbits;
943 goto L10;
944 }
945 if (lexp == -(*emin)) {
946 uexp = lexp;
947 } else {
948 uexp = try__;
949 ++exbits;
950 }
951
952/* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
953/* than or equal to EMIN. EXBITS is the number of bits needed to */
954/* store the exponent. */
955
956 if (uexp + *emin > -lexp - *emin) {
957 expsum = lexp << 1;
958 } else {
959 expsum = uexp << 1;
960 }
961
962/* EXPSUM is the exponent range, approximately equal to */
963/* EMAX - EMIN + 1 . */
964
965 *emax = expsum + *emin - 1;
966 nbits = exbits + 1 + *p;
967
968/* NBITS is the total number of bits needed to store a */
969/* floating-point number. */
970
971 if (nbits % 2 == 1 && *beta == 2) {
972
973/* Either there are an odd number of bits used to store a */
974/* floating-point number, which is unlikely, or some bits are */
975/* not used in the representation of numbers, which is possible, */
976/* (e.g. Cray machines) or the mantissa has an implicit bit, */
977/* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
978/* most likely. We have to assume the last alternative. */
979/* If this is true, then we need to reduce EMAX by one because */
980/* there must be some way of representing zero in an implicit-bit */
981/* system. On machines like Cray, we are reducing EMAX by one */
982/* unnecessarily. */
983
984 --(*emax);
985 }
986
987 if (*ieee) {
988
989/* Assume we are on an IEEE machine which reserves one exponent */
990/* for infinity and NaN. */
991
992 --(*emax);
993 }
994
995/* Now create RMAX, the largest machine number, which should */
996/* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
997
998/* First compute 1.0 - BETA**(-P), being careful that the */
999/* result is less than 1.0 . */
1000
1001 recbas = (float)1. / *beta;
1002 z__ = *beta - (float)1.;
1003 y = (float)0.;
1004 i__1 = *p;
1005 for (i__ = 1; i__ <= i__1; ++i__) {
1006 z__ *= recbas;
1007 if (y < (float)1.) {
1008 oldy = y;
1009 }
1010 y = fla_slamc3(&y, &z__);
1011/* L20: */
1012 }
1013 if (y >= (float)1.) {
1014 y = oldy;
1015 }
1016
1017/* Now multiply by BETA**EMAX to get RMAX. */
1018
1019 i__1 = *emax;
1020 for (i__ = 1; i__ <= i__1; ++i__) {
1021 r__1 = y * *beta;
1022 y = fla_slamc3(&r__1, &c_b32);
1023/* L30: */
1024 }
1025
1026 *rmax = y;
1027 return 0;
1028
1029/* End of SLAMC5 */
1030
1031} /* fla_slamc5_ */

References fla_slamc3(), and i.

Referenced by fla_slamc2().

◆ fla_slamch()

real fla_slamch ( char *  cmach,
ftnlen  cmach_len 
)
57{
58 /* Initialized data */
59
60 static logical first = TRUE_;
61
62 /* System generated locals */
63 integer i__1;
64 real ret_val;
65
66 /* Builtin functions */
67 double fla_pow_ri(real *, integer *);
68
69 /* Local variables */
70 static real base;
71 static integer beta;
72 static real emin, prec, emax;
73 static integer imin, imax;
74 static logical lrnd;
75 static real rmin, rmax, t, rmach;
76 extern logical fla_lsame(char *, char *, ftnlen, ftnlen);
77 static real small, sfmin;
78 extern /* Subroutine */ int fla_slamc2(integer *, integer *, logical *, real
79 *, integer *, real *, integer *, real *);
80 static integer it;
81 static real rnd, eps;
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/* SLAMCH determines single precision machine parameters. */
95
96/* Arguments */
97/* ========= */
98
99/* CMACH (input) CHARACTER*1 */
100/* Specifies the value to be returned by SLAMCH: */
101/* = 'E' or 'e', SLAMCH := eps */
102/* = 'S' or 's , SLAMCH := sfmin */
103/* = 'B' or 'b', SLAMCH := base */
104/* = 'P' or 'p', SLAMCH := eps*base */
105/* = 'N' or 'n', SLAMCH := t */
106/* = 'R' or 'r', SLAMCH := rnd */
107/* = 'M' or 'm', SLAMCH := emin */
108/* = 'U' or 'u', SLAMCH := rmin */
109/* = 'L' or 'l', SLAMCH := emax */
110/* = 'O' or 'o', SLAMCH := 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_slamc2(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
143 base = (real) beta;
144 t = (real) it;
145 if (lrnd) {
146 rnd = (float)1.;
147 i__1 = 1 - it;
148 eps = fla_pow_ri(&base, &i__1) / 2;
149 } else {
150 rnd = (float)0.;
151 i__1 = 1 - it;
152 eps = fla_pow_ri(&base, &i__1);
153 }
154 prec = eps * base;
155 emin = (real) imin;
156 emax = (real) imax;
157 sfmin = rmin;
158 small = (float)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 + (float)1.);
165 }
166 }
167
168 if (fla_lsame(cmach, "E", (ftnlen)1, (ftnlen)1)) {
169 rmach = eps;
170 } else if (fla_lsame(cmach, "S", (ftnlen)1, (ftnlen)1)) {
171 rmach = sfmin;
172 } else if (fla_lsame(cmach, "B", (ftnlen)1, (ftnlen)1)) {
173 rmach = base;
174 } else if (fla_lsame(cmach, "P", (ftnlen)1, (ftnlen)1)) {
175 rmach = prec;
176 } else if (fla_lsame(cmach, "N", (ftnlen)1, (ftnlen)1)) {
177 rmach = t;
178 } else if (fla_lsame(cmach, "R", (ftnlen)1, (ftnlen)1)) {
179 rmach = rnd;
180 } else if (fla_lsame(cmach, "M", (ftnlen)1, (ftnlen)1)) {
181 rmach = emin;
182 } else if (fla_lsame(cmach, "U", (ftnlen)1, (ftnlen)1)) {
183 rmach = rmin;
184 } else if (fla_lsame(cmach, "L", (ftnlen)1, (ftnlen)1)) {
185 rmach = emax;
186 } else if (fla_lsame(cmach, "O", (ftnlen)1, (ftnlen)1)) {
187 rmach = rmax;
188 }
189
190 ret_val = rmach;
191 first = FALSE_;
192 return ret_val;
193
194/* End of SLAMCH */
195
196} /* fla_slamch_ */
ftnlen
short ftnlen
Definition FLA_f2c.h:61
fla_lsame
logical fla_lsame(char *ca, char *cb, ftnlen ca_len, ftnlen cb_len)
Definition fla_lsame.c:20
fla_slamc2
int fla_slamc2(integer *beta, integer *t, logical *rnd, real *eps, integer *emin, real *rmin, integer *emax, real *rmax)
Definition fla_slamch.c:409

References fla_lsame(), fla_pow_ri(), fla_slamc2(), and i.

Referenced by FLA_Mach_params_ops().

Generated on Thu Sep 23 2021 17:46:17 for libflame by doxygen 1.9.8