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_qp
21 :
22 : use const_def, only : qp
23 :
24 : use ISO_FORTRAN_ENV
25 : use ISO_C_BINDING
26 :
27 : implicit none
28 :
29 : integer, parameter :: FRAC_BITS_128_H = 48
30 : integer, parameter :: FRAC_BITS_128_L = 64
31 : integer, parameter :: EXPN_BITS_128_H = 15
32 :
33 : integer(INT64), parameter :: QNAN_128_H = INT(z'7fff100000000000', INT64)
34 : integer(INT64), parameter :: QNAN_128_L = INT(z'0000000000000001', INT64)
35 : integer(INT64), parameter :: SNAN_128_H = INT(z'7fff000000000000', INT64)
36 : integer(INT64), parameter :: SNAN_128_L = INT(z'0000000000000001', INT64)
37 :
38 : interface is_nan
39 : module procedure is_nan_qp
40 : end interface is_nan
41 :
42 : interface is_inf
43 : module procedure is_inf_qp
44 : end interface is_inf
45 :
46 : interface is_bad
47 : module procedure is_bad_qp
48 : end interface is_bad
49 :
50 : interface set_nan
51 : module procedure set_nan_qp_0d
52 : module procedure set_nan_qp_1d
53 : module procedure set_nan_qp_2d
54 : module procedure set_nan_qp_3d
55 : module procedure set_nan_qp_4d
56 : end interface set_nan
57 :
58 : private
59 :
60 : public :: is_nan
61 : public :: is_inf
62 : public :: is_bad
63 : public :: set_nan
64 :
65 : contains
66 :
67 120738 : elemental function is_nan_qp (x, signal) result (is_nan)
68 :
69 : real(qp), target, intent(in) :: x
70 : logical, optional, intent(in) :: signal
71 : logical :: is_nan
72 :
73 : integer(INT64) :: ix(2)
74 : integer(INT64) :: frac_l
75 : integer(INT64) :: frac_h
76 : integer(INT64) :: expn
77 : integer(INT64) :: sign
78 :
79 : ! Convert x to integer
80 :
81 362214 : ix = TRANSFER(x, ix)
82 :
83 : ! Split out IEEE fields
84 :
85 120738 : frac_l = IBITS(ix(1), 0, FRAC_BITS_128_L)
86 120738 : frac_h = IBITS(ix(2), 0, FRAC_BITS_128_H)
87 120738 : expn = IBITS(ix(2), FRAC_BITS_128_H, EXPN_BITS_128_H)
88 120738 : sign = IBITS(ix(2), FRAC_BITS_128_H+EXPN_BITS_128_H, 1)
89 :
90 : ! Check for NaN
91 :
92 : is_nan = expn == MASKR(EXPN_BITS_128_H, INT64) .AND. &
93 120738 : (frac_h /= 0_INT64 .OR. frac_l /= 0_INT64)
94 :
95 120738 : if (PRESENT(signal)) then
96 0 : is_nan = is_nan .AND. (BTEST(frac_h, FRAC_BITS_128_H-1) .EQV. .NOT. signal)
97 : end if
98 :
99 : return
100 :
101 : end function is_nan_qp
102 :
103 :
104 120738 : elemental function is_inf_qp (x) result (is_inf)
105 :
106 : real(qp), target, intent(in) :: x
107 : logical :: is_inf
108 :
109 : integer(INT64) :: ix(2)
110 : integer(INT64) :: frac_l
111 : integer(INT64) :: frac_h
112 : integer(INT64) :: expn
113 : integer(INT64) :: sign
114 :
115 : ! Convert x to integer
116 :
117 362214 : ix = TRANSFER(x, ix)
118 :
119 120738 : frac_l = IBITS(ix(1), 0, FRAC_BITS_128_L)
120 120738 : frac_h = IBITS(ix(2), 0, FRAC_BITS_128_H)
121 120738 : expn = IBITS(ix(2), FRAC_BITS_128_H, EXPN_BITS_128_H)
122 120738 : sign = IBITS(ix(2), FRAC_BITS_128_H+EXPN_BITS_128_H, 1)
123 :
124 : ! Check for infinity
125 :
126 : is_inf = expn == MASKR(EXPN_BITS_128_H, INT64) .AND. &
127 120738 : (frac_l == 0_INT64 .AND. frac_h == 0_INT64)
128 :
129 : return
130 :
131 : end function is_inf_qp
132 :
133 :
134 120738 : elemental function is_bad_qp (x) result (is_bad)
135 :
136 : real(qp), intent(in) :: x
137 : logical :: is_bad
138 :
139 : ! Check for NaN or infinity
140 :
141 120738 : is_bad = is_nan(x) .OR. is_inf(x)
142 :
143 : return
144 :
145 : end function is_bad_qp
146 :
147 :
148 0 : subroutine set_nan_qp_0d (x, signal)
149 :
150 : real(qp), target, intent(out) :: x
151 : logical, optional, intent(in) :: signal
152 :
153 : integer(INT64), pointer :: ix(:)
154 :
155 : ! Convert x to a signaling or quiet NaN
156 :
157 0 : call C_F_POINTER(C_LOC(x), ix, [2])
158 :
159 0 : if (PRESENT(signal)) then
160 0 : if (signal) then
161 0 : ix(1) = SNAN_128_L
162 0 : ix(2) = SNAN_128_H
163 : else
164 0 : ix(1) = QNAN_128_L
165 0 : ix(2) = QNAN_128_H
166 : end if
167 : else
168 0 : ix(1) = SNAN_128_L
169 0 : ix(2) = SNAN_128_H
170 : end if
171 :
172 0 : return
173 :
174 : end subroutine set_nan_qp_0d
175 :
176 :
177 0 : subroutine set_nan_qp_1d (x, signal)
178 :
179 : real(qp), target, intent(out) :: x(:)
180 : logical, optional, intent(in) :: signal
181 :
182 : integer(INT64), pointer :: ix(:,:)
183 :
184 : ! Convert x to a signaling or quiet NaN
185 :
186 0 : call C_F_POINTER(C_LOC(x), ix, [2,SHAPE(x)])
187 :
188 0 : if (PRESENT(signal)) then
189 0 : if (signal) then
190 0 : ix(1,:) = SNAN_128_L
191 0 : ix(2,:) = SNAN_128_H
192 : else
193 0 : ix(1,:) = QNAN_128_L
194 0 : ix(2,:) = QNAN_128_H
195 : end if
196 : else
197 0 : ix(1,:) = SNAN_128_L
198 0 : ix(2,:) = SNAN_128_H
199 : end if
200 :
201 0 : return
202 :
203 : end subroutine set_nan_qp_1d
204 :
205 :
206 0 : subroutine set_nan_qp_2d (x, signal)
207 :
208 : real(qp), target, intent(out) :: x(:,:)
209 : logical, optional, intent(in) :: signal
210 :
211 : integer(INT64), pointer :: ix(:,:,:)
212 :
213 : ! Convert x to a signaling or quiet NaN
214 :
215 0 : call C_F_POINTER(C_LOC(x), ix, [2,SHAPE(x)])
216 :
217 0 : if (PRESENT(signal)) then
218 0 : if (signal) then
219 0 : ix(1,:,:) = SNAN_128_L
220 0 : ix(2,:,:) = SNAN_128_H
221 : else
222 0 : ix(1,:,:) = QNAN_128_L
223 0 : ix(2,:,:) = QNAN_128_H
224 : end if
225 : else
226 0 : ix(1,:,:) = SNAN_128_L
227 0 : ix(2,:,:) = SNAN_128_H
228 : end if
229 :
230 0 : return
231 :
232 : end subroutine set_nan_qp_2d
233 :
234 :
235 0 : subroutine set_nan_qp_3d (x, signal)
236 :
237 : real(qp), target, intent(out) :: x(:,:,:)
238 : logical, optional, intent(in) :: signal
239 :
240 : integer(INT64), pointer :: ix(:,:,:,:)
241 :
242 : ! Convert x to a signaling or quiet NaN
243 :
244 0 : call C_F_POINTER(C_LOC(x), ix, [2,SHAPE(x)])
245 :
246 0 : if (PRESENT(signal)) then
247 0 : if (signal) then
248 0 : ix(1,:,:,:) = SNAN_128_L
249 0 : ix(2,:,:,:) = SNAN_128_H
250 : else
251 0 : ix(1,:,:,:) = QNAN_128_L
252 0 : ix(2,:,:,:) = QNAN_128_H
253 : end if
254 : else
255 0 : ix(1,:,:,:) = SNAN_128_L
256 0 : ix(2,:,:,:) = SNAN_128_H
257 : end if
258 :
259 0 : return
260 :
261 : end subroutine set_nan_qp_3d
262 :
263 :
264 0 : subroutine set_nan_qp_4d (x, signal)
265 :
266 : real(qp), target, intent(out) :: x(:,:,:,:)
267 : logical, optional, intent(in) :: signal
268 :
269 : integer(INT64), pointer :: ix(:,:,:,:,:)
270 :
271 : ! Convert x to a signaling or quiet NaN
272 :
273 0 : call C_F_POINTER(C_LOC(x), ix, [2,SHAPE(x)])
274 :
275 0 : if (PRESENT(signal)) then
276 0 : if (signal) then
277 0 : ix(1,:,:,:,:) = SNAN_128_L
278 0 : ix(2,:,:,:,:) = SNAN_128_H
279 : else
280 0 : ix(1,:,:,:,:) = QNAN_128_L
281 0 : ix(2,:,:,:,:) = QNAN_128_H
282 : end if
283 : else
284 0 : ix(1,:,:,:,:) = SNAN_128_L
285 0 : ix(2,:,:,:,:) = SNAN_128_H
286 : end if
287 :
288 0 : return
289 :
290 : end subroutine set_nan_qp_4d
291 :
292 0 : end module utils_nan_qp
|
