Parent Directory | Revision Log

Revision **1** -
(**show annotations**)
(**download**)
(**as text**)

*Fri Oct 29 20:54:12 2004 UTC*
(15 years, 9 months ago)
by *aw0a*

File MIME type: text/x-csrc

File size: 3651 byte(s)

File MIME type: text/x-csrc

File size: 3651 byte(s)

Setting up web subdirectory in repository

1 | /*********************************************************************\ |

2 | d1mach: Ascend Replacement for d1mach.f |

3 | by Ben Allan |

4 | Created: September, 1994 |

5 | Version: $Revision: 1.3 $ |

6 | Date last modified: $Date: 1998/07/06 10:56:12 $ |

7 | |

8 | This file is part of the Ascend fortran subroutine collection. |

9 | |

10 | Copyright (C) Benjamin Andrew Allan |

11 | |

12 | The ascend fortran subroutine collection is free software; you can redistribute |

13 | it and/or modify it under the terms of the GNU General Public License as |

14 | published by the Free Software Foundation; either version 2 of the |

15 | License, or (at your option) any later version. |

16 | Most of the sources in the ascend fortran subroutine collection are public |

17 | domain and available from NETLIB. See newsgroup sci.math.numerical-analysis. |

18 | Sources from netlib are not restricted by the GNU license and are marked as |

19 | such. |

20 | |

21 | The Ascend fortran subroutine collection is distributed in hope that it will be |

22 | useful, but WITHOUT ANY WARRANTY; without even the implied warranty of |

23 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

24 | General Public License for more details. |

25 | |

26 | You should have received a copy of the GNU General Public License along with |

27 | the program; if not, write to the Free Software Foundation, Inc., 675 |

28 | Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING. |

29 | COPYING is found in ../compiler. |

30 | \*********************************************************************/ |

31 | |

32 | /* d1mach.c. Ben Allan |

33 | C replacement for d1mach.f in terms of ANSI constants. |

34 | The LINPACK d1mach.f is not such that f77 compilers pick |

35 | the right set of constants automatically. |

36 | |

37 | Observed equivalences: |

38 | F D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. |

39 | C DBL_MIN |

40 | F D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. |

41 | C DBL_MAX |

42 | F D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. |

43 | C DBL_EPSILON/2 |

44 | F D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. |

45 | C DBL_EPSILON |

46 | F D1MACH( 5) = LOG10(B) |

47 | C NONE |

48 | B: FLT_RADIX |

49 | EMIN: DBL_MIN_EXP |

50 | EMAX: DBL_MAX_EXP |

51 | T: DBL_MANT_DIG |

52 | |

53 | On alphas d1mach(3)=DBL_EPSILON/2 for some reason. Returning |

54 | DBL_EPSILON may result in 1 bit of conservatism in some codes, but |

55 | this is the price of portability. |

56 | |

57 | */ |

58 | #include <stdlib.h> |

59 | #include <stdio.h> |

60 | #include <float.h> |

61 | #include <math.h> |

62 | |

63 | |

64 | /* Commentary and Apology: |

65 | * We used to have a bunch of #ifdef's here trying to figure out which |

66 | * platform we were on and then blessing our d1mach function with the |

67 | * proper number of underbars (d1mach vs d1mach_) so that the linker |

68 | * would not whine about missing symbols due to the insanity of |

69 | * whether or not the f77 compiler puts an underbar on the symbols it |

70 | * generates. Of course, just to make life fun, it's not strictly |

71 | * platform dependent, since some f77 compilers accept flags that turn |

72 | * the underbars on or off. Given this lunacy and the wasted time of |

73 | * trying to tack down this bug every time it occurs, we've decided to |

74 | * just duplicate the function and be finished with the underbar |

75 | * madness. We realize this sucks and we apologize for it, but at |

76 | * this point,``Frankly, my dears, we don't give a damn.'' |

77 | */ |

78 | |

79 | double d1mach(int *i) { |

80 | switch (*i) { |

81 | case 1: |

82 | return DBL_MIN; |

83 | case 2: |

84 | return DBL_MAX; |

85 | case 3: |

86 | return DBL_EPSILON; |

87 | case 4: |

88 | return DBL_EPSILON; |

89 | case 5: |

90 | return log10((double)FLT_RADIX); |

91 | default: |

92 | fprintf(stderr," D1MACH - I OUT OF BOUNDS %d",*i); |

93 | abort(); |

94 | } |

95 | } |

96 | |

97 | double d1mach_(int *i) { |

98 | return d1mach(i); |

99 | } |

100 | |

101 | double D1MACH(int *i) { |

102 | return d1mach(i); |

103 | } |

104 | |

105 | |

106 |

john.pye@anu.edu.au | ViewVC Help |

Powered by ViewVC 1.1.22 |