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_sp
21 :
22 : use const_def, only : sp
23 :
24 : use ISO_FORTRAN_ENV
25 : use ISO_C_BINDING
26 :
27 : implicit none
28 :
29 : integer, parameter :: FRAC_BITS_32 = 23
30 : integer, parameter :: EXPN_BITS_32 = 8
31 :
32 : integer(INT32), parameter :: QNAN_32 = INT(z'7fc00001', INT32)
33 : integer(INT32), parameter :: SNAN_32 = INT(z'7f800001', INT32)
34 :
35 : interface is_nan
36 : module procedure is_nan_sp
37 : end interface is_nan
38 :
39 : interface is_inf
40 : module procedure is_inf_sp
41 : end interface is_inf
42 :
43 : interface is_bad
44 : module procedure is_bad_sp
45 : end interface is_bad
46 :
47 : interface set_nan
48 : module procedure set_nan_sp_0d
49 : module procedure set_nan_sp_1d
50 : module procedure set_nan_sp_2d
51 : module procedure set_nan_sp_3d
52 : module procedure set_nan_sp_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 0 : elemental function is_nan_sp (x, signal) result (is_nan)
65 :
66 : real(sp), target, intent(in) :: x
67 : logical, optional, intent(in) :: signal
68 : logical :: is_nan
69 :
70 : integer(INT32) :: ix
71 : integer(INT32) :: frac
72 : integer(INT32) :: expn
73 : integer(INT32) :: sign
74 :
75 : ! Convert x to integer
76 :
77 0 : ix = TRANSFER(x, ix)
78 :
79 : ! Split out IEEE fields
80 :
81 0 : frac = IBITS(ix, 0, FRAC_BITS_32)
82 0 : expn = IBITS(ix, FRAC_BITS_32, EXPN_BITS_32)
83 0 : sign = IBITS(ix, FRAC_BITS_32+EXPN_BITS_32, 1)
84 :
85 : ! Check for NaN
86 :
87 0 : is_nan = expn == MASKR(EXPN_BITS_32, INT32) .AND. frac /= 0_INT32
88 :
89 0 : if (PRESENT(signal)) then
90 0 : is_nan = is_nan .AND. (BTEST(frac, FRAC_BITS_32-1) .EQV. .NOT. signal)
91 : end if
92 :
93 :
94 : return
95 :
96 : end function is_nan_sp
97 :
98 :
99 0 : elemental function is_inf_sp (x) result (is_inf)
100 :
101 : real(sp), target, intent(in) :: x
102 : logical :: is_inf
103 :
104 : integer(INT32) :: ix
105 : integer(INT32) :: frac
106 : integer(INT32) :: expn
107 : integer(INT32) :: sign
108 :
109 : ! Convert x to integer
110 :
111 0 : ix = TRANSFER(x, ix)
112 :
113 : ! Split out IEEE fields
114 :
115 0 : frac = IBITS(ix, 0, FRAC_BITS_32)
116 0 : expn = IBITS(ix, FRAC_BITS_32, EXPN_BITS_32)
117 0 : sign = IBITS(ix, FRAC_BITS_32+EXPN_BITS_32, 1)
118 :
119 : ! Check for infinity
120 :
121 0 : is_inf = expn == MASKR(EXPN_BITS_32, INT32) .AND. frac == 0_INT32
122 :
123 :
124 : return
125 :
126 : end function is_inf_sp
127 :
128 :
129 0 : elemental function is_bad_sp (x) result (is_bad)
130 :
131 : real(sp), intent(in) :: x
132 : logical :: is_bad
133 :
134 : ! Check for NaN or infinity
135 :
136 0 : is_bad = is_nan(x) .OR. is_inf(x)
137 :
138 :
139 : return
140 :
141 : end function is_bad_sp
142 :
143 :
144 0 : subroutine set_nan_sp_0d (x, signal)
145 :
146 : real(sp), target, intent(out) :: x
147 : logical, optional, intent(in) :: signal
148 :
149 : integer(INT32), pointer :: ix
150 :
151 : ! Convert x to a signaling or quiet NaN
152 :
153 0 : call C_F_POINTER(C_LOC(x), ix)
154 :
155 0 : if (PRESENT(signal)) then
156 0 : if (signal) then
157 0 : ix = SNAN_32
158 : else
159 0 : ix = QNAN_32
160 : end if
161 : else
162 0 : ix = SNAN_32
163 : end if
164 :
165 0 : end subroutine set_nan_sp_0d
166 :
167 :
168 0 : subroutine set_nan_sp_1d (x, signal)
169 :
170 : real(sp), target, intent(out) :: x(:)
171 : logical, optional, intent(in) :: signal
172 :
173 : integer(INT32), pointer :: ix(:)
174 :
175 : ! Convert x to a signaling or quiet NaN
176 :
177 0 : call C_F_POINTER(C_LOC(x), ix, SHAPE(x))
178 :
179 0 : if (PRESENT(signal)) then
180 0 : if (signal) then
181 0 : ix = SNAN_32
182 : else
183 0 : ix = QNAN_32
184 : end if
185 : else
186 0 : ix = SNAN_32
187 : end if
188 :
189 :
190 0 : return
191 :
192 : end subroutine set_nan_sp_1d
193 :
194 :
195 0 : subroutine set_nan_sp_2d (x, signal)
196 :
197 : real(sp), target, intent(out) :: x(:,:)
198 : logical, optional, intent(in) :: signal
199 :
200 : integer(INT32), pointer :: ix(:,:)
201 :
202 : ! Convert x to a signaling or quiet NaN
203 :
204 0 : call C_F_POINTER(C_LOC(x), ix, SHAPE(x))
205 :
206 0 : if (PRESENT(signal)) then
207 0 : if (signal) then
208 0 : ix = SNAN_32
209 : else
210 0 : ix = QNAN_32
211 : end if
212 : else
213 0 : ix = SNAN_32
214 : end if
215 :
216 :
217 0 : return
218 :
219 : end subroutine set_nan_sp_2d
220 :
221 :
222 0 : subroutine set_nan_sp_3d (x, signal)
223 :
224 : real(sp), target, intent(out) :: x(:,:,:)
225 : logical, optional, intent(in) :: signal
226 :
227 : integer(INT32), pointer :: ix(:,:,:)
228 :
229 : ! Convert x to a signaling or quiet NaN
230 :
231 0 : call C_F_POINTER(C_LOC(x), ix, SHAPE(x))
232 :
233 0 : if (PRESENT(signal)) then
234 0 : if (signal) then
235 0 : ix = SNAN_32
236 : else
237 0 : ix = QNAN_32
238 : end if
239 : else
240 0 : ix = SNAN_32
241 : end if
242 :
243 :
244 0 : return
245 :
246 : end subroutine set_nan_sp_3d
247 :
248 :
249 0 : subroutine set_nan_sp_4d (x, signal)
250 :
251 : real(sp), target, intent(out) :: x(:,:,:,:)
252 : logical, optional, intent(in) :: signal
253 :
254 : integer(INT32), pointer :: ix(:,:,:,:)
255 :
256 : ! Convert x to a signaling or quiet NaN
257 :
258 0 : call C_F_POINTER(C_LOC(x), ix, SHAPE(x))
259 :
260 0 : if (PRESENT(signal)) then
261 0 : if (signal) then
262 0 : ix = SNAN_32
263 : else
264 0 : ix = QNAN_32
265 : end if
266 : else
267 0 : ix = SNAN_32
268 : end if
269 :
270 :
271 0 : return
272 :
273 : end subroutine set_nan_sp_4d
274 :
275 0 : end module utils_nan_sp
|
