libflame  revision_anchor
Functions
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()

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

References fla_dlamc3().

Referenced by fla_dlamc2().

203 {
204  /* Initialized data */
205 
206  static logical first = TRUE_;
207 
208  /* System generated locals */
209  doublereal d__1, d__2;
210 
211  /* Local variables */
212  static logical lrnd;
213  static doublereal a, b, c__, f;
214  static integer lbeta;
215  static doublereal savec;
217  static logical lieee1;
218  static doublereal t1, t2;
219  static integer lt;
220  static doublereal 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 /* 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 */
299 L10:
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 */
318 L20:
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 */
383 L30:
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_ */
double doublereal
Definition: FLA_f2c.h:31
int logical
Definition: FLA_f2c.h:36
int integer
Definition: FLA_f2c.h:25
doublereal fla_dlamc3(doublereal *a, doublereal *b)
Definition: fla_dlamch.c:726

◆ fla_dlamc2()

int fla_dlamc2 ( integer beta,
integer t,
logical rnd,
doublereal eps,
integer emin,
doublereal rmin,
integer emax,
doublereal rmax 
)

References fla_dlamc1(), fla_dlamc3(), fla_dlamc4(), fla_dlamc5(), and fla_pow_di().

Referenced by fla_dlamch().

414 {
415  /* Initialized data */
416 
417  static logical first = TRUE_;
418  static logical iwarn = FALSE_;
419 
420  /* Format strings */
421  static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre\
422 ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the value EMIN loo\
423 ks\002,\002 acceptable please comment out \002,/\002 the IF block as marked \
424 within the code of routine\002,\002 DLAMC2,\002,/\002 otherwise supply EMIN \
425 explicitly.\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 */
432  double fla_pow_di(doublereal *, integer *);
433  //integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe();
434 
435  /* Local variables */
436  static logical ieee;
437  static doublereal half;
438  static logical lrnd;
439  static doublereal leps, zero, a, b, c__;
440  static integer i__, lbeta;
441  static doublereal rbase;
442  static integer lemin, lemax, gnmin;
443  static doublereal small;
444  static integer gpmin;
445  static doublereal third, lrmin, lrmax, sixth;
446  extern /* Subroutine */ int fla_dlamc1(integer *, integer *, logical *,
447  logical *);
449  static logical lieee1;
450  extern /* Subroutine */ int fla_dlamc4(integer *, doublereal *, integer *),
452  doublereal *);
453  static integer lt, ngnmin, ngpmin;
454  static doublereal one, two;
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, &lt, &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 */
575 L10:
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, &lt, &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_dlamc4(integer *emin, doublereal *start, integer *base)
Definition: fla_dlamch.c:768
int fla_dlamc5(integer *beta, integer *p, integer *emin, logical *ieee, integer *emax, doublereal *rmax)
Definition: fla_dlamch.c:866
double doublereal
Definition: FLA_f2c.h:31
int logical
Definition: FLA_f2c.h:36
int integer
Definition: FLA_f2c.h:25
double fla_pow_di(doublereal *ap, integer *bp)
Definition: fla_dlamch.c:26
int fla_dlamc1(integer *beta, integer *t, logical *rnd, logical *ieee1)
Definition: fla_dlamch.c:201
doublereal fla_dlamc3(doublereal *a, doublereal *b)
Definition: fla_dlamch.c:726

◆ fla_dlamc3()

doublereal fla_dlamc3 ( doublereal a,
doublereal b 
)

Referenced by fla_dlamc1(), fla_dlamc2(), fla_dlamc4(), and fla_dlamc5().

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_ */
double doublereal
Definition: FLA_f2c.h:31

◆ fla_dlamc4()

int fla_dlamc4 ( integer emin,
doublereal start,
integer base 
)

References fla_dlamc3().

Referenced by fla_dlamc2().

769 {
770  /* System generated locals */
771  integer i__1;
772  doublereal d__1;
773 
774  /* Local variables */
775  static doublereal zero, a;
776  static integer i__;
777  static doublereal rbase, b1, b2, c1, c2, d1, d2;
779  static doublereal one;
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 */
829 L10:
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_ */
double doublereal
Definition: FLA_f2c.h:31
int integer
Definition: FLA_f2c.h:25
doublereal fla_dlamc3(doublereal *a, doublereal *b)
Definition: fla_dlamch.c:726

◆ fla_dlamc5()

int fla_dlamc5 ( integer beta,
integer p,
integer emin,
logical ieee,
integer emax,
doublereal rmax 
)

References fla_dlamc3().

Referenced by fla_dlamc2().

868 {
869  /* System generated locals */
870  integer i__1;
871  doublereal d__1;
872 
873  /* Local variables */
874  static integer lexp;
875  static doublereal oldy;
876  static integer uexp, i__;
877  static doublereal y, z__;
878  static integer nbits;
880  static doublereal recbas;
881  static integer exbits, expsum, try__;
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;
943 L10:
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_ */
double doublereal
Definition: FLA_f2c.h:31
int integer
Definition: FLA_f2c.h:25
doublereal fla_dlamc3(doublereal *a, doublereal *b)
Definition: fla_dlamch.c:726

◆ fla_dlamch()

doublereal fla_dlamch ( char *  cmach,
ftnlen  cmach_len 
)

References fla_dlamc2(), fla_lsame(), and fla_pow_di().

57 {
58  /* Initialized data */
59 
60  static logical first = TRUE_;
61 
62  /* System generated locals */
63  integer i__1;
64  doublereal ret_val;
65 
66  /* Builtin functions */
67  double fla_pow_di(doublereal *, integer *);
68 
69  /* Local variables */
70  static doublereal base;
71  static integer beta;
72  static doublereal emin, prec, emax;
73  static integer imin, imax;
74  static logical lrnd;
75  static doublereal rmin, rmax, t, rmach;
76  extern logical fla_lsame(char *, char *, ftnlen, ftnlen);
77  static doublereal small, sfmin;
78  extern /* Subroutine */ int fla_dlamc2(integer *, integer *, logical *,
80  static integer it;
81  static doublereal 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 /* 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 
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 DLAMCH */
195 
196 } /* fla_dlamch_ */
short ftnlen
Definition: FLA_f2c.h:61
double doublereal
Definition: FLA_f2c.h:31
int logical
Definition: FLA_f2c.h:36
int integer
Definition: FLA_f2c.h:25
double fla_pow_di(doublereal *ap, integer *bp)
Definition: fla_dlamch.c:26
int fla_dlamc2(integer *beta, integer *t, logical *rnd, doublereal *eps, integer *emin, doublereal *rmin, integer *emax, doublereal *rmax)
Definition: fla_dlamch.c:411
logical fla_lsame(char *ca, char *cb, ftnlen ca_len, ftnlen cb_len)
Definition: fla_lsame.c:20

◆ fla_pow_di()

double fla_pow_di ( doublereal ap,
integer bp 
)

Referenced by fla_dlamc2(), and fla_dlamch().

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