Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module utils_nan_dp
21 :
22 : use const_def, only : dp
23 :
24 : use ISO_FORTRAN_ENV
25 : use ISO_C_BINDING
26 :
27 : implicit none
28 :
29 : integer, parameter :: FRAC_BITS_64 = 52
30 : integer, parameter :: EXPN_BITS_64 = 11
31 :
32 : integer(INT64), parameter :: QNAN_64 = INT(z'7ff8000000000001', INT64)
33 : integer(INT64), parameter :: SNAN_64 = INT(z'7ff0000000000001', INT64)
34 :
35 : interface is_nan
36 : module procedure is_nan_dp
37 : end interface is_nan
38 :
39 : interface is_inf
40 : module procedure is_inf_dp
41 : end interface is_inf
42 :
43 : interface is_bad
44 : module procedure is_bad_dp
45 : end interface is_bad
46 :
47 : interface set_nan
48 : module procedure set_nan_dp_0d
49 : module procedure set_nan_dp_1d
50 : module procedure set_nan_dp_2d
51 : module procedure set_nan_dp_3d
52 : module procedure set_nan_dp_4d
53 : end interface set_nan
54 :
55 : private
56 :
57 : public :: is_nan
58 : public :: is_inf
59 : public :: is_bad
60 : public :: set_nan
61 :
62 : contains
63 :
64 2044783171 : elemental function is_nan_dp (x, signal) result (is_nan)
65 :
66 : real(dp), target, intent(in) :: x
67 : logical, optional, intent(in) :: signal
68 : logical :: is_nan
69 :
70 : integer(INT64) :: ix
71 : integer(INT64) :: frac
72 : integer(INT64) :: expn
73 : integer(INT64) :: sign
74 :
75 : ! Convert x to integer
76 :
77 2044783171 : ix = TRANSFER(x, ix)
78 :
79 : ! Split out IEEE fields
80 :
81 2044783171 : frac = IBITS(ix, 0, FRAC_BITS_64)
82 2044783171 : expn = IBITS(ix, FRAC_BITS_64, EXPN_BITS_64)
83 2044783171 : sign = IBITS(ix, FRAC_BITS_64+EXPN_BITS_64, 1)
84 :
85 : ! Check for NaN
86 :
87 2044783171 : is_nan = expn == MASKR(EXPN_BITS_64, INT64) .AND. frac /= 0_INT64
88 :
89 0 : if (PRESENT(signal)) then
90 0 : is_nan = is_nan .AND. (BTEST(frac, FRAC_BITS_64-1) .EQV. .NOT. signal)
91 : end if
92 :
93 : return
94 :
95 : end function is_nan_dp
96 :
97 :
98 2044783171 : elemental function is_inf_dp (x) result (is_inf)
99 :
100 : real(dp), target, intent(in) :: x
101 : logical :: is_inf
102 :
103 : integer(INT64) :: ix
104 : integer(INT64) :: frac
105 : integer(INT64) :: expn
106 : integer(INT64) :: sign
107 :
108 : ! Convert x to integer
109 :
110 2044783171 : ix = TRANSFER(x, ix)
111 :
112 : ! Split out IEEE fields
113 :
114 2044783171 : frac = IBITS(ix, 0, FRAC_BITS_64)
115 2044783171 : expn = IBITS(ix, FRAC_BITS_64, EXPN_BITS_64)
116 2044783171 : sign = IBITS(ix, FRAC_BITS_64+EXPN_BITS_64, 1)
117 :
118 : ! Check for infinity
119 :
120 2044783171 : is_inf = expn == MASKR(EXPN_BITS_64, INT64) .AND. frac == 0_INT64
121 :
122 : return
123 :
124 : end function is_inf_dp
125 :
126 :
127 2044783171 : elemental function is_bad_dp (x) result (is_bad)
128 :
129 : real(dp), intent(in) :: x
130 : logical :: is_bad
131 :
132 : ! Check for NaN or infinity
133 :
134 2044783171 : is_bad = is_nan(x) .OR. is_inf(x)
135 :
136 : return
137 :
138 : end function is_bad_dp
139 :
140 :
141 0 : subroutine set_nan_dp_0d (x, signal)
142 :
143 : real(dp), target, intent(out) :: x
144 : logical, optional, intent(in) :: signal
145 :
146 : integer(INT64), pointer :: ix
147 :
148 : ! Convert x to a signaling or quiet NaN
149 :
150 0 : call C_F_POINTER(C_LOC(x), ix)
151 :
152 0 : if (PRESENT(signal)) then
153 0 : if (signal) then
154 0 : ix = SNAN_64
155 : else
156 0 : ix = QNAN_64
157 : end if
158 : else
159 0 : ix = SNAN_64
160 : end if
161 :
162 0 : return
163 :
164 : end subroutine set_nan_dp_0d
165 :
166 :
167 0 : subroutine set_nan_dp_1d (x, signal)
168 :
169 : real(dp), target, intent(out) :: x(:)
170 : logical, optional, intent(in) :: signal
171 :
172 : integer(INT64), pointer :: ix(:)
173 :
174 : ! Convert x to a signaling or quiet NaN
175 :
176 0 : call C_F_POINTER(C_LOC(x), ix, SHAPE(x))
177 :
178 0 : if (PRESENT(signal)) then
179 0 : if (signal) then
180 0 : ix = SNAN_64
181 : else
182 0 : ix = QNAN_64
183 : end if
184 : else
185 0 : ix = SNAN_64
186 : end if
187 :
188 0 : return
189 :
190 : end subroutine set_nan_dp_1d
191 :
192 :
193 19 : subroutine set_nan_dp_2d (x, signal)
194 :
195 : real(dp), target, intent(out) :: x(:,:)
196 : logical, optional, intent(in) :: signal
197 :
198 : integer(INT64), pointer :: ix(:,:)
199 :
200 : ! Convert x to a signaling or quiet NaN
201 :
202 57 : call C_F_POINTER(C_LOC(x), ix, SHAPE(x))
203 :
204 19 : if (PRESENT(signal)) then
205 0 : if (signal) then
206 0 : ix = SNAN_64
207 : else
208 0 : ix = QNAN_64
209 : end if
210 : else
211 10641159 : ix = SNAN_64
212 : end if
213 :
214 19 : return
215 :
216 : end subroutine set_nan_dp_2d
217 :
218 :
219 0 : subroutine set_nan_dp_3d (x, signal)
220 :
221 : real(dp), target, intent(out) :: x(:,:,:)
222 : logical, optional, intent(in) :: signal
223 :
224 : integer(INT64), pointer :: ix(:,:,:)
225 :
226 : ! Convert x to a signaling or quiet NaN
227 :
228 0 : call C_F_POINTER(C_LOC(x), ix, SHAPE(x))
229 :
230 0 : if (PRESENT(signal)) then
231 0 : if (signal) then
232 0 : ix = SNAN_64
233 : else
234 0 : ix = QNAN_64
235 : end if
236 : else
237 0 : ix = SNAN_64
238 : end if
239 :
240 0 : return
241 :
242 : end subroutine set_nan_dp_3d
243 :
244 :
245 138 : subroutine set_nan_dp_4d (x, signal)
246 :
247 : real(dp), target, intent(out) :: x(:,:,:,:)
248 : logical, optional, intent(in) :: signal
249 :
250 : integer(INT64), pointer :: ix(:,:,:,:)
251 :
252 : ! Convert x to a signaling or quiet NaN
253 :
254 690 : call C_F_POINTER(C_LOC(x), ix, SHAPE(x))
255 :
256 138 : if (PRESENT(signal)) then
257 0 : if (signal) then
258 0 : ix = SNAN_64
259 : else
260 0 : ix = QNAN_64
261 : end if
262 : else
263 19887042 : ix = SNAN_64
264 : end if
265 :
266 138 : return
267 :
268 : end subroutine set_nan_dp_4d
269 :
270 0 : end module utils_nan_dp
|
