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 
)

Referenced by fla_slamc2(), and fla_slamch().

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 }
int integer
Definition: FLA_f2c.h:25

◆ fla_slamc1()

int fla_slamc1 ( integer beta,
integer t,
logical rnd,
logical ieee1 
)

References fla_slamc3().

Referenced by fla_slamc2().

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 */
299 L10:
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 */
318 L20:
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 */
382 L30:
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_ */
float real
Definition: FLA_f2c.h:30
int logical
Definition: FLA_f2c.h:36
int integer
Definition: FLA_f2c.h:25
real fla_slamc3(real *a, real *b)
Definition: fla_slamch.c:721

◆ fla_slamc2()

int fla_slamc2 ( integer beta,
integer t,
logical rnd,
real eps,
integer emin,
real rmin,
integer emax,
real rmax 
)

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

Referenced by fla_slamch().

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\
419 ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the value EMIN loo\
420 ks\002,\002 acceptable please comment out \002,/\002 the IF block as marked \
421 within the code of routine\002,\002 SLAMC2,\002,/\002 otherwise supply EMIN \
422 explicitly.\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 *),
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 */
572 L10:
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_ */
int fla_slamc1(integer *beta, integer *t, logical *rnd, logical *ieee1)
Definition: fla_slamch.c:201
float real
Definition: FLA_f2c.h:30
double fla_pow_ri(real *ap, integer *bp)
Definition: fla_slamch.c:26
int logical
Definition: FLA_f2c.h:36
int integer
Definition: FLA_f2c.h:25
real fla_slamc3(real *a, real *b)
Definition: fla_slamch.c:721
int fla_slamc4(integer *emin, real *start, integer *base)
Definition: fla_slamch.c:763
int fla_slamc5(integer *beta, integer *p, integer *emin, logical *ieee, integer *emax, real *rmax)
Definition: fla_slamch.c:861

◆ fla_slamc3()

real fla_slamc3 ( real a,
real b 
)

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

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_ */
float real
Definition: FLA_f2c.h:30

◆ fla_slamc4()

int fla_slamc4 ( integer emin,
real start,
integer base 
)

References fla_slamc3().

Referenced by fla_slamc2().

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 */
824 L10:
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_ */
float real
Definition: FLA_f2c.h:30
int integer
Definition: FLA_f2c.h:25
real fla_slamc3(real *a, real *b)
Definition: fla_slamch.c:721

◆ fla_slamc5()

int fla_slamc5 ( integer beta,
integer p,
integer emin,
logical ieee,
integer emax,
real rmax 
)

References fla_slamc3().

Referenced by fla_slamc2().

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;
938 L10:
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_ */
float real
Definition: FLA_f2c.h:30
int integer
Definition: FLA_f2c.h:25
real fla_slamc3(real *a, real *b)
Definition: fla_slamch.c:721

◆ fla_slamch()

real fla_slamch ( char *  cmach,
ftnlen  cmach_len 
)

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

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_ */
short ftnlen
Definition: FLA_f2c.h:61
float real
Definition: FLA_f2c.h:30
int fla_slamc2(integer *beta, integer *t, logical *rnd, real *eps, integer *emin, real *rmin, integer *emax, real *rmax)
Definition: fla_slamch.c:409
double fla_pow_ri(real *ap, integer *bp)
Definition: fla_slamch.c:26
int logical
Definition: FLA_f2c.h:36
int integer
Definition: FLA_f2c.h:25
logical fla_lsame(char *ca, char *cb, ftnlen ca_len, ftnlen cb_len)
Definition: fla_lsame.c:20