portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> A module with intrinsic functions that are used by MOM but are not supported | |
| 6 | !! by some compilers. | |
| 7 | module MOM_intrinsic_functions | |
| 8 | ||
| 9 | use iso_fortran_env, only : stdout => output_unit, stderr => error_unit | |
| 10 | use iso_fortran_env, only : int64, real64 | |
| 11 | ||
| 12 | implicit none ; private | |
| 13 | ||
| 14 | public :: invcosh, cuberoot, nth_root | |
| 15 | public :: intrinsic_functions_unit_tests | |
| 16 | ||
| 17 | ! Floating point model, if bit layout from high to low is (sign, exp, frac) | |
| 18 | ||
| 19 | integer, parameter :: bias = maxexponent(1.) - 1 | |
| 20 | !< The double precision exponent offset | |
| 21 | integer, parameter :: signbit = storage_size(1.) - 1 | |
| 22 | !< Position of sign bit | |
| 23 | integer, parameter :: explen = 1 + ceiling(log(real(bias))/log(2.)) | |
| 24 | !< Bit size of exponent | |
| 25 | integer, parameter :: expbit = signbit - explen | |
| 26 | !< Position of lowest exponent bit | |
| 27 | integer, parameter :: fraclen = expbit | |
| 28 | !< Length of fractional part | |
| 29 | ||
| 30 | contains | |
| 31 | ||
| 32 | !> Evaluate the inverse cosh, either using a math library or an | |
| 33 | !! equivalent expression | |
| 34 | 87120 | function invcosh(x) |
| 35 | real, intent(in) :: x !< The argument of the inverse of cosh [nondim]. NaNs will | |
| 36 | !! occur if x<1, but there is no error checking | |
| 37 | real :: invcosh ! The inverse of cosh of x [nondim] | |
| 38 | ||
| 39 | #ifdef __INTEL_COMPILER | |
| 40 | invcosh = acosh(x) | |
| 41 | #else | |
| 42 | 87120 | invcosh = log(x+sqrt(x*x-1)) |
| 43 | #endif | |
| 44 | ||
| 45 | 87120 | end function invcosh |
| 46 | ||
| 47 | ||
| 48 | !> Returns the cube root of a real argument at roundoff accuracy, in a form that works properly with | |
| 49 | !! rescaling of the argument by integer powers of 8. If the argument is a NaN, a NaN is returned. | |
| 50 | 9674344 | elemental function cuberoot(x) result(root) |
| 51 | !$omp declare target | |
| 52 | real, intent(in) :: x !< The argument of cuberoot in arbitrary units cubed [A3] | |
| 53 | real :: root !< The real cube root of x in arbitrary units [A] | |
| 54 | ||
| 55 | real :: asx ! The absolute value of x rescaled by an integer power of 8 to put it into | |
| 56 | ! the range from 0.125 < asx <= 1.0, in ambiguous units cubed [B3] | |
| 57 | real :: root_asx ! The cube root of asx [B] | |
| 58 | real :: ra_3 ! root_asx cubed [B3] | |
| 59 | real :: num ! The numerator of an expression for the evolving estimate of the cube root of asx | |
| 60 | ! in arbitrary units that can grow or shrink with each iteration [B C] | |
| 61 | real :: den ! The denominator of an expression for the evolving estimate of the cube root of asx | |
| 62 | ! in arbitrary units that can grow or shrink with each iteration [C] | |
| 63 | real :: num_prev ! The numerator of an expression for the previous iteration of the evolving estimate | |
| 64 | ! of the cube root of asx in arbitrary units that can grow or shrink with each iteration [B D] | |
| 65 | real :: np_3 ! num_prev cubed [B3 D3] | |
| 66 | real :: den_prev ! The denominator of an expression for the previous iteration of the evolving estimate of | |
| 67 | ! the cube root of asx in arbitrary units that can grow or shrink with each iteration [D] | |
| 68 | real :: dp_3 ! den_prev cubed [C3] | |
| 69 | real :: r0 ! Initial value of the iterative solver. [B C] | |
| 70 | real :: r0_3 ! r0 cubed [B3 C3] | |
| 71 | integer :: itt | |
| 72 | ||
| 73 | integer(kind=int64) :: e_x, s_x | |
| 74 | ||
| 75 | 9674344 | if ((x >= 0.0) .eqv. (x <= 0.0)) then |
| 76 | ! Return 0 for an input of 0, or NaN for a NaN input. | |
| 77 | 5856 | root = x |
| 78 | else | |
| 79 | 9668488 | call rescale_cbrt(x, asx, e_x, s_x) |
| 80 | ||
| 81 | ! Iteratively determine root_asx = asx**1/3 using Halley's method and then Newton's method, | |
| 82 | ! noting that Halley's method onverges monotonically and needs no bounding. Halley's method is | |
| 83 | ! slightly more complicated that Newton's method, but converges in a third fewer iterations. | |
| 84 | ! Keeping the estimates in a fractional form Root = num / den allows this calculation with | |
| 85 | ! no real divisions during the iterations before doing a single real division at the end, | |
| 86 | ! and it is therefore more computationally efficient. | |
| 87 | ||
| 88 | ! This first estimate gives the same magnitude of errors for 0.125 and 1.0 after two iterations. | |
| 89 | ! The first iteration is applied explicitly. | |
| 90 | 9668488 | r0 = 0.707106 |
| 91 | 9668488 | r0_3 = r0 * r0 * r0 |
| 92 | 9668488 | num = r0 * (r0_3 + 2.0 * asx) |
| 93 | 9668488 | den = 2.0 * r0_3 + asx |
| 94 | ||
| 95 | 29005464 | do itt=1,2 |
| 96 | ! Halley's method iterates estimates as Root = Root * (Root**3 + 2.*asx) / (2.*Root**3 + asx). | |
| 97 | 19336976 | num_prev = num ; den_prev = den |
| 98 | ||
| 99 | ! Pre-compute these as integer powers, to avoid `pow()`-like intrinsics. | |
| 100 | 19336976 | np_3 = num_prev * num_prev * num_prev |
| 101 | 19336976 | dp_3 = den_prev * den_prev * den_prev |
| 102 | ||
| 103 | 19336976 | num = num_prev * (np_3 + 2.0 * asx * dp_3) |
| 104 | 29005464 | den = den_prev * (2.0 * np_3 + asx * dp_3) |
| 105 | ! Equivalent to: root_asx = root_asx * (root_asx**3 + 2.*asx) / (2.*root_asx**3 + asx) | |
| 106 | enddo | |
| 107 | ! At this point the error in root_asx is better than 1 part in 3e14. | |
| 108 | 9668488 | root_asx = num / den |
| 109 | ||
| 110 | ! One final iteration with Newton's method polishes up the root and gives a solution | |
| 111 | ! that is within the last bit of the true solution. | |
| 112 | 9668488 | ra_3 = root_asx * root_asx * root_asx |
| 113 | 9668488 | root_asx = root_asx - (ra_3 - asx) / (3.0 * (root_asx * root_asx)) |
| 114 | ||
| 115 | 9668488 | root = descale(root_asx, e_x, s_x) |
| 116 | endif | |
| 117 | 9674344 | end function cuberoot |
| 118 | ||
| 119 | ||
| 120 | !> Bit-stable n-th root of x for x in (0, +inf) and integer n >= 1, suitable | |
| 121 | !! for evaluation inside `!$omp target` / `do concurrent` offloaded regions. | |
| 122 | !! | |
| 123 | !! Lowering `x**(1.0/n)` via the compiler produces `exp((1.0/n)*log(x))` — two | |
| 124 | !! transcendentals whose last-bit rounding differs between host libm and CUDA | |
| 125 | !! libdevice. This routine avoids that path entirely: it uses fixed-iteration | |
| 126 | !! Newton on y^n - x = 0, with y^(n-1) evaluated as repeated multiplication, | |
| 127 | !! and one bit-precision-polishing iteration at the end. | |
| 128 | !! | |
| 129 | !! For x in [0, 1] (the case in `MOM_barotropic.F90`'s `bt_rem = av_rem**Instep`) | |
| 130 | !! convergence is rapid because the linear initial guess y0 = 1 - (1-x)/n is | |
| 131 | !! already within a few percent of the true root. | |
| 132 | 0 | elemental function nth_root(x, n) result(root) |
| 133 | !$omp declare target | |
| 134 | real, intent(in) :: x !< Argument, x >= 0 [arbitrary] | |
| 135 | integer, intent(in) :: n !< Root index, n >= 1 | |
| 136 | real :: root !< x**(1/n) in the same units as x | |
| 137 | ||
| 138 | integer, parameter :: maxitt = 20 ! Fixed (deterministic) iteration count | |
| 139 | real :: y, ypow_nm1 | |
| 140 | real :: x_n_r, x_nm1_r | |
| 141 | integer :: itt, k | |
| 142 | ||
| 143 | ! Trivial cases — return early to keep loop tight and avoid 0/0 below. | |
| 144 | 0 | if (n <= 1) then |
| 145 | 0 | root = x |
| 146 | 0 | return |
| 147 | endif | |
| 148 | 0 | if (x == 0.0) then |
| 149 | 0 | root = 0.0 |
| 150 | 0 | return |
| 151 | endif | |
| 152 | ||
| 153 | 0 | x_n_r = real(n) |
| 154 | 0 | x_nm1_r = real(n - 1) |
| 155 | ||
| 156 | ! Linear initial guess valid for x in [0, 1] and decent for moderate x > 1. | |
| 157 | ! For our caller (av_rem in [0, 1], typically near 1), this is within ~1%. | |
| 158 | 0 | y = 1.0 - (1.0 - x) / x_n_r |
| 159 | ||
| 160 | ! Newton iteration: y_{k+1} = ((n-1)*y_k + x / y_k^{n-1}) / n | |
| 161 | ! All ops are *, +, /. Integer power y^{n-1} is repeated multiplication, | |
| 162 | ! so there is no `pow`/`exp/log` lowering anywhere in the iteration. | |
| 163 | 0 | do itt = 1, maxitt |
| 164 | 0 | ypow_nm1 = 1.0 |
| 165 | 0 | do k = 1, n - 1 |
| 166 | 0 | ypow_nm1 = ypow_nm1 * y |
| 167 | enddo | |
| 168 | 0 | y = (x_nm1_r * y + x / ypow_nm1) / x_n_r |
| 169 | enddo | |
| 170 | ||
| 171 | 0 | root = y |
| 172 | 0 | end function nth_root |
| 173 | ||
| 174 | ||
| 175 | !> Rescale `a` to the range [0.125, 1) and compute its cube-root exponent. | |
| 176 | !! | |
| 177 | !! This function decomposes `a` into the form `s * x * 2**e` so that `x` is | |
| 178 | !! in the desired range. This is accomplished by computing the integral cube | |
| 179 | !! root of `e` (as a division) and applying the residual to `x`. | |
| 180 | 9668488 | pure subroutine rescale_cbrt(a, x, e_r, s_a) |
| 181 | !$omp declare target | |
| 182 | real, intent(in) :: a | |
| 183 | !< The number to be rescaled for cube-root computation [A3] | |
| 184 | real, intent(out) :: x | |
| 185 | !< The rescaled value of `a` in the range [0.125, 1) [B3] | |
| 186 | integer(kind=int64), intent(out) :: e_r | |
| 187 | !< The integral component of the cube-root exponent of `a`. | |
| 188 | integer(kind=int64), intent(out) :: s_a | |
| 189 | !< Sign bit of `a`. A nonzero value indicates negative sign. | |
| 190 | ||
| 191 | integer(kind=int64) :: xb | |
| 192 | ! Floating point integer representation of `a` | |
| 193 | integer(kind=int64) :: e_a | |
| 194 | ! Exponent of `a` | |
| 195 | integer(kind=int64) :: e_x | |
| 196 | ! Exponent of `x` | |
| 197 | ||
| 198 | ! Pack bits of a into xb and extract its exponent and sign. | |
| 199 | 9668488 | xb = transfer(a, 1_int64) |
| 200 | 9668488 | s_a = ibits(xb, signbit, 1) |
| 201 | 9668488 | e_a = ibits(xb, expbit, explen) - bias |
| 202 | ||
| 203 | ! The floating-point form of `a` with exponent `e` is | |
| 204 | ! | |
| 205 | ! a = s * (1 + m) * 2**e | |
| 206 | ! | |
| 207 | ! where (1+m) ∈ [1,2). We want to split 2**e so that (1+m) is rescaled to | |
| 208 | ! the range [0.125, 1); that is, [2**-3, 2**0). | |
| 209 | ! | |
| 210 | ! First decompose the exponent `e` into quotient-remainder form: | |
| 211 | ! | |
| 212 | ! e = 3⌊e/3⌋ + modulo(e,3) | |
| 213 | ! | |
| 214 | ! Since modulo(e,3) ∈ {0,1,2}, the second term of the following expression is | |
| 215 | ! in {-3,-2,-1}. | |
| 216 | ! | |
| 217 | ! e = 3 * (⌊e/3⌋ + 1) + (modulo(e,3) - 3). | |
| 218 | ! | |
| 219 | ! Here, (modulo(e,3) - 3) is in the range [2**-3, 1) and holds the | |
| 220 | ! floating-point exponent of `x`. | |
| 221 | ! | |
| 222 | ! Fortran integer division is round-to-zero. To convert to floor division, | |
| 223 | ! we use the sign() intrinsic to shift negative values downward. | |
| 224 | ! | |
| 225 | ! ⌊e/3⌋ = (e + sign(1,e) - 1) / 3 | |
| 226 | ! | |
| 227 | ! ⌊e/3⌋ + 1 reduces to the form below. This is what we call the integral | |
| 228 | ! cube-root of `a` in the description above. | |
| 229 | ||
| 230 | 9668488 | e_r = (e_a + sign(1_int64, e_a) + 2) / 3 |
| 231 | ||
| 232 | ! modulo() is not implemented on all systems, so compute the remainder as | |
| 233 | ! r = n - 3*q. | |
| 234 | ||
| 235 | 9668488 | e_x = e_a - e_r * 3 |
| 236 | ||
| 237 | ! Insert the new 11-bit exponent into xb and write to x and extend the | |
| 238 | ! bitcount to 12, so that the sign bit is zero and x is always positive. | |
| 239 | 9668488 | call mvbits(e_x + bias, 0, explen + 1, xb, fraclen) |
| 240 | 9668488 | x = transfer(xb, 1.) |
| 241 | 9668488 | end subroutine rescale_cbrt |
| 242 | ||
| 243 | ||
| 244 | !> Undo the rescaling of a real number back to its original base. | |
| 245 | 9668488 | pure function descale(x, e_a, s_a) result(a) |
| 246 | !$omp declare target | |
| 247 | real, intent(in) :: x | |
| 248 | !< The rescaled value which is to be restored in ambiguous units [B] | |
| 249 | integer(kind=int64), intent(in) :: e_a | |
| 250 | !< Exponent of the unscaled value | |
| 251 | integer(kind=int64), intent(in) :: s_a | |
| 252 | !< Sign bit of the unscaled value | |
| 253 | real :: a | |
| 254 | !< Restored value with the corrected exponent and sign in arbitrary units [A] | |
| 255 | ||
| 256 | integer(kind=int64) :: xb | |
| 257 | ! Bit-packed real number into integer form | |
| 258 | integer(kind=int64) :: e_x | |
| 259 | ! Biased exponent of x | |
| 260 | ||
| 261 | ! Apply the corrected exponent and sign to x. | |
| 262 | 9668488 | xb = transfer(x, 1_int64) |
| 263 | 9668488 | e_x = ibits(xb, expbit, explen) |
| 264 | 9668488 | call mvbits(e_a + e_x, 0, explen, xb, expbit) |
| 265 | 9668488 | call mvbits(s_a, 0, 1, xb, signbit) |
| 266 | 9668488 | a = transfer(xb, 1.) |
| 267 | 9668488 | end function descale |
| 268 | ||
| 269 | ||
| 270 | !> Returns true if any unit test of intrinsic_functions fails, or false if they all pass. | |
| 271 | 0 | function intrinsic_functions_unit_tests(verbose) result(fail) |
| 272 | logical, intent(in) :: verbose !< If true, write results to stdout | |
| 273 | logical :: fail !< True if any of the unit tests fail | |
| 274 | ||
| 275 | ! Local variables | |
| 276 | real :: testval ! A test value for self-consistency testing [nondim] | |
| 277 | logical :: v | |
| 278 | integer :: n | |
| 279 | ||
| 280 | 0 | fail = .false. |
| 281 | 0 | v = verbose |
| 282 | 0 | write(stdout,*) '==== MOM_intrinsic_functions: intrinsic_functions_unit_tests ===' |
| 283 | ||
| 284 | 0 | fail = fail .or. Test_cuberoot(v, 1.2345678901234e9) |
| 285 | 0 | fail = fail .or. Test_cuberoot(v, -9.8765432109876e-21) |
| 286 | 0 | fail = fail .or. Test_cuberoot(v, 64.0) |
| 287 | 0 | fail = fail .or. Test_cuberoot(v, -0.5000000000001) |
| 288 | 0 | fail = fail .or. Test_cuberoot(v, 0.0) |
| 289 | 0 | fail = fail .or. Test_cuberoot(v, 1.0) |
| 290 | 0 | fail = fail .or. Test_cuberoot(v, 0.125) |
| 291 | 0 | fail = fail .or. Test_cuberoot(v, 0.965) |
| 292 | 0 | fail = fail .or. Test_cuberoot(v, 1.0 - epsilon(1.0)) |
| 293 | 0 | fail = fail .or. Test_cuberoot(v, 1.0 - 0.5*epsilon(1.0)) |
| 294 | ||
| 295 | 0 | testval = 1.0e-99 |
| 296 | 0 | v = .false. |
| 297 | 0 | do n=-160,160 |
| 298 | 0 | fail = fail .or. Test_cuberoot(v, testval) |
| 299 | 0 | testval = (-2.908 * (1.414213562373 + 1.2345678901234e-5*n)) * testval |
| 300 | enddo | |
| 301 | 0 | end function intrinsic_functions_unit_tests |
| 302 | ||
| 303 | !> True if the cube of cuberoot(val) does not closely match val. False otherwise. | |
| 304 | 0 | logical function Test_cuberoot(verbose, val) |
| 305 | logical, intent(in) :: verbose !< If true, write results to stdout | |
| 306 | real, intent(in) :: val !< The real value to test, in arbitrary units [A] | |
| 307 | ! Local variables | |
| 308 | real :: diff ! The difference between val and the cube root of its cube [A]. | |
| 309 | ||
| 310 | 0 | diff = val - cuberoot(val)**3 |
| 311 | 0 | Test_cuberoot = (abs(diff) > 2.0e-15*abs(val)) |
| 312 | ||
| 313 | 0 | if (Test_cuberoot) then |
| 314 | 0 | write(stdout, '("For val = ",ES22.15,", (val - cuberoot(val**3))) = ",ES9.2," <-- FAIL")') val, diff |
| 315 | 0 | elseif (verbose) then |
| 316 | 0 | write(stdout, '("For val = ",ES22.15,", (val - cuberoot(val**3))) = ",ES9.2)') val, diff |
| 317 | ||
| 318 | endif | |
| 319 | 0 | end function Test_cuberoot |
| 320 | ||
| 321 | end module MOM_intrinsic_functions |