other/m_common.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Common module for the BSpline FEEC package | ||
| 2 | !> | ||
| 3 | !> This module contains common definitions, constants, and interfaces used throughout the package. | ||
| 4 | module m_common | ||
| 5 | |||
| 6 | private | ||
| 7 | public :: wp, pi, factorial, levi_civita, set_union, sort, uppercase, to_lowercase, zero_function, cumsum, prime_factors, balanced_split | ||
| 8 | public :: user_function_1d_interface, user_function_3d_interface, user_function_3d_vec_interface, user_function_3d_mat_interface | ||
| 9 | public :: cartesian_to_linear, brent | ||
| 10 | public :: init_random_seed | ||
| 11 | public :: VERBOSITY_WARN_NEVER, VERBOSITY_WARN_ON_FAILURE, & | ||
| 12 | VERBOSITY_WARN_ALWAYS, VERBOSITY_ERROR_ON_FAILURE | ||
| 13 | |||
| 14 | !> Single precision real numbers | ||
| 15 | integer, parameter :: single_precision = selected_real_kind(6) | ||
| 16 | |||
| 17 | !> Double precision real numbers | ||
| 18 | integer, parameter :: double_precision = selected_real_kind(15) | ||
| 19 | |||
| 20 | !> Default value which unused if CMake is used | ||
| 21 | #ifndef CMAKE_BSPLINE_FEEC_PRECISION | ||
| 22 | #define CMAKE_BSPLINE_FEEC_PRECISION double_precision | ||
| 23 | #endif | ||
| 24 | |||
| 25 | !> The working precision for real-valued numbers | ||
| 26 | integer, parameter :: wp = CMAKE_BSPLINE_FEEC_PRECISION | ||
| 27 | |||
| 28 | !> The value of pi (3.14159...) | ||
| 29 | real(wp), parameter :: pi = 3.1415926535897932384626433832795028841971_wp | ||
| 30 | |||
| 31 | !> Interface for user-defined functions | ||
| 32 | abstract interface | ||
| 33 | pure function user_function_1d_interface(x) result(ans) | ||
| 34 | import wp | ||
| 35 | implicit none | ||
| 36 | real(wp), intent(in) :: x | ||
| 37 | real(wp) :: ans | ||
| 38 | end function | ||
| 39 | |||
| 40 | pure function user_function_3d_interface(x, y, z) result(ans) | ||
| 41 | import wp | ||
| 42 | implicit none | ||
| 43 | real(wp), intent(in) :: x, y, z | ||
| 44 | real(wp) :: ans | ||
| 45 | end function | ||
| 46 | |||
| 47 | pure function user_function_3d_vec_interface(ind, x, y, z) result(ans) | ||
| 48 | import wp | ||
| 49 | implicit none | ||
| 50 | integer, intent(in) :: ind | ||
| 51 | real(wp), intent(in) :: x, y, z | ||
| 52 | real(wp) :: ans | ||
| 53 | end function | ||
| 54 | |||
| 55 | pure function user_function_3d_mat_interface(s1, s2, x, y, z) result(ans) | ||
| 56 | import wp | ||
| 57 | implicit none | ||
| 58 | integer, intent(in) :: s1, s2 | ||
| 59 | real(wp), intent(in) :: x, y, z | ||
| 60 | real(wp) :: ans | ||
| 61 | end function | ||
| 62 | end interface | ||
| 63 | |||
| 64 | ! TODO use enums or some other type for verbosity levels | ||
| 65 | !> Verbosity levels for warnings and errors | ||
| 66 | integer, parameter :: VERBOSITY_WARN_NEVER = 0 | ||
| 67 | integer, parameter :: VERBOSITY_WARN_ON_FAILURE = 1 | ||
| 68 | integer, parameter :: VERBOSITY_WARN_ALWAYS = 2 | ||
| 69 | integer, parameter :: VERBOSITY_ERROR_ON_FAILURE = 3 | ||
| 70 | |||
| 71 | ! Module procedure interfaces for submodules | ||
| 72 | interface | ||
| 73 | !> @brief Sort an array of integers using the quicksort algorithm | ||
| 74 | !> | ||
| 75 | !> @param[inout] array The array of integers to be sorted | ||
| 76 | !> @param[in] inverse _(optional)_ If true, sort in descending order; otherwise, ascending order (default: false) | ||
| 77 | module pure subroutine sort_int(array, inverse) | ||
| 78 | integer, intent(inout) :: array(:) | ||
| 79 | logical, intent(in), optional :: inverse | ||
| 80 | end subroutine | ||
| 81 | |||
| 82 | !> @brief Sort an array of real numbers using the quicksort algorithm | ||
| 83 | !> | ||
| 84 | !> @param[inout] array The array of real numbers to be sorted | ||
| 85 | !> @param[in] inverse _(optional)_ If true, sort in descending order; otherwise, ascending order (default: false) | ||
| 86 | !> @param[in] sort_fun _(optional)_ A function to determine the sorting order | ||
| 87 | module subroutine sort_real(array, inverse, sort_fun) | ||
| 88 | real(wp), intent(inout) :: array(:) | ||
| 89 | logical, intent(in), optional :: inverse | ||
| 90 | procedure(user_function_1d_interface), optional :: sort_fun | ||
| 91 | end subroutine | ||
| 92 | |||
| 93 | !> @brief Compute the union of two sets of integers | ||
| 94 | !> | ||
| 95 | !> @param[out] union The resulting union of the two sets | ||
| 96 | !> @param[in] set1 The first set of integers | ||
| 97 | !> @param[in] set2 The second set of integers | ||
| 98 | module subroutine set_union(union, set1, set2) | ||
| 99 | integer, allocatable, intent(out) :: union(:) | ||
| 100 | integer, intent(in) :: set1(0:), set2(0:) | ||
| 101 | end subroutine | ||
| 102 | |||
| 103 | !> @brief Compute the prime factorization of an integer | ||
| 104 | !> | ||
| 105 | !> @param[in] n The integer to factorize | ||
| 106 | !> @param[out] factors Array containing the prime factors (with repetition) | ||
| 107 | module pure subroutine prime_factors(n, factors) | ||
| 108 | integer, intent(in) :: n | ||
| 109 | integer, allocatable, intent(out) :: factors(:) | ||
| 110 | end subroutine | ||
| 111 | |||
| 112 | !> @brief Split an integer into two factors that are as close as possible | ||
| 113 | !> | ||
| 114 | !> Given an integer N, find n1 and n2 such that N = n1 * n2 and |n1 - n2| is minimized. | ||
| 115 | !> The result satisfies n1 <= n2. | ||
| 116 | !> | ||
| 117 | !> @param[in] n The integer to split | ||
| 118 | !> @param[out] n1 The smaller factor | ||
| 119 | !> @param[out] n2 The larger factor | ||
| 120 | module pure subroutine balanced_split_2(n, n1, n2) | ||
| 121 | integer, intent(in) :: n | ||
| 122 | integer, intent(out) :: n1, n2 | ||
| 123 | end subroutine | ||
| 124 | |||
| 125 | !> @brief Split an integer into three factors that are as close as possible | ||
| 126 | !> | ||
| 127 | !> Given an integer N, find n1, n2, and n3 such that N = n1 * n2 * n3 and the factors are as balanced as possible. | ||
| 128 | !> The result satisfies n1 <= n2 <= n3. | ||
| 129 | !> | ||
| 130 | !> @param[in] n The integer to split | ||
| 131 | !> @param[out] n1 The smallest factor | ||
| 132 | !> @param[out] n2 The middle factor | ||
| 133 | !> @param[out] n3 The largest factor | ||
| 134 | module pure subroutine balanced_split_3(n, n1, n2, n3) | ||
| 135 | integer, intent(in) :: n | ||
| 136 | integer, intent(out) :: n1, n2, n3 | ||
| 137 | end subroutine | ||
| 138 | |||
| 139 | |||
| 140 | !> @brief Find a root of a function using Brent's method | ||
| 141 | !> | ||
| 142 | !> @param[in] fun The function for which to find a root | ||
| 143 | !> @param[in] xa_in The lower bound of the interval | ||
| 144 | !> @param[in] xb_in The upper bound of the interval | ||
| 145 | !> @param[in] x_tol The tolerance for convergence | ||
| 146 | !> @param[in] max_iter The maximum number of iterations to perform | ||
| 147 | !> @param[in] fa_in _(optional)_ The value of the function at xa_in (if not provided, it will be computed) | ||
| 148 | !> @param[in] fb_in _(optional)_ The value of the function at xb_in (if not provided, it will be computed) | ||
| 149 | !> @param[out] success _(optional)_ A logical flag indicating whether the method converged successfully | ||
| 150 | !> | ||
| 151 | !> @return The estimated root of the function | ||
| 152 | module real(wp) function brent(fun, xa_in, xb_in, x_tol, max_iter, fa_in, fb_in, success) result(xsol) | ||
| 153 | implicit none | ||
| 154 | |||
| 155 | procedure(user_function_1d_interface) :: fun | ||
| 156 | real(wp), intent(in) :: xa_in, xb_in, x_tol | ||
| 157 | integer, intent(in) :: max_iter | ||
| 158 | real(wp), intent(in), optional :: fa_in, fb_in | ||
| 159 | logical, intent(out), optional :: success | ||
| 160 | |||
| 161 | end function brent | ||
| 162 | end interface | ||
| 163 | |||
| 164 | !> @brief Generic interface for sorting arrays | ||
| 165 | interface sort | ||
| 166 | module procedure sort_int, sort_real | ||
| 167 | end interface sort | ||
| 168 | |||
| 169 | !> @brief Generic interface for balanced splitting of integers | ||
| 170 | interface balanced_split | ||
| 171 | module procedure balanced_split_2, balanced_split_3 | ||
| 172 | end interface balanced_split | ||
| 173 | |||
| 174 | contains | ||
| 175 | |||
| 176 | !> @brief A function that always returns zero | ||
| 177 | !> | ||
| 178 | !> @param[in] xp The x-coordinate | ||
| 179 | !> @param[in] yp The y-coordinate | ||
| 180 | !> @param[in] zp The z-coordinate | ||
| 181 | !> | ||
| 182 | !> @return Always returns 0.0 | ||
| 183 | 43909152 | pure real(wp) function zero_function(xp, yp, zp) result(val) | |
| 184 | real(wp), intent(in) :: xp, yp, zp | ||
| 185 | |||
| 186 | val = 0._wp | ||
| 187 | 43909152 | end function zero_function | |
| 188 | |||
| 189 | !> @brief Compute the factorial of a non-negative integer n | ||
| 190 | !> | ||
| 191 | !> @param[in] n The non-negative integer for which to compute the factorial | ||
| 192 | !> | ||
| 193 | !> @return The factorial of n | ||
| 194 | 16 | pure integer function factorial(n) result(ans) | |
| 195 | implicit none | ||
| 196 | |||
| 197 | integer, intent(in) :: n | ||
| 198 | |||
| 199 | integer :: i | ||
| 200 | |||
| 201 | ans = 1 | ||
| 202 |
2/2✓ Branch 3 → 4 taken 120 times.
✓ Branch 3 → 5 taken 16 times.
|
136 | do i = 2, n |
| 203 | 136 | ans = ans * i | |
| 204 | end do | ||
| 205 | 16 | end function | |
| 206 | |||
| 207 | !> @brief Compute the Levi-Civita symbol for three indices | ||
| 208 | !> | ||
| 209 | !> The Levi-Civita symbol is defined as: | ||
| 210 | !> +1 if (i, j, k) is an even permutation of (1, 2, 3) | ||
| 211 | !> -1 if (i, j, k) is an odd permutation of (1, 2, 3) | ||
| 212 | !> | ||
| 213 | !> @param[in] i The first index | ||
| 214 | !> @param[in] j The second index | ||
| 215 | !> @param[in] k The third index | ||
| 216 | !> | ||
| 217 | !> @return The value of the Levi-Civita symbol | ||
| 218 | 98223 | pure integer function levi_civita(i, j, k) | |
| 219 | !$acc routine seq | ||
| 220 | implicit none | ||
| 221 | |||
| 222 | integer, intent(in) :: i, j, k | ||
| 223 | |||
| 224 |
6/6✓ Branch 2 → 3 taken 97974 times.
✓ Branch 2 → 24 taken 249 times.
✓ Branch 3 → 4 taken 65316 times.
✓ Branch 3 → 24 taken 32658 times.
✓ Branch 4 → 5 taken 32658 times.
✓ Branch 4 → 24 taken 32658 times.
|
98223 | if (i == j .or. j == k .or. i == k) then |
| 225 | levi_civita = 0 | ||
| 226 | else if ((i == 1 .and. j == 2 .and. k == 3) .or. & | ||
| 227 |
15/18✓ Branch 5 → 6 taken 12204 times.
✓ Branch 5 → 8 taken 20454 times.
✓ Branch 6 → 7 taken 6124 times.
✓ Branch 6 → 8 taken 6080 times.
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 24 taken 6124 times.
✓ Branch 8 → 9 taken 10205 times.
✓ Branch 8 → 11 taken 16329 times.
✓ Branch 9 → 10 taken 5080 times.
✓ Branch 9 → 11 taken 5125 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 24 taken 5080 times.
✓ Branch 11 → 12 taken 10249 times.
✓ Branch 11 → 14 taken 11205 times.
✓ Branch 12 → 13 taken 5125 times.
✓ Branch 12 → 14 taken 5124 times.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 24 taken 5125 times.
|
32658 | & (i == 2 .and. j == 3 .and. k == 1) .or. & |
| 228 | & (i == 3 .and. j == 1 .and. k == 2)) then | ||
| 229 | levi_civita = 1 | ||
| 230 | else if ((i == 1 .and. j == 3 .and. k == 2) .or. & | ||
| 231 |
11/18✓ Branch 14 → 15 taken 6080 times.
✓ Branch 14 → 17 taken 10249 times.
✓ Branch 15 → 16 taken 6080 times.
✗ Branch 15 → 17 not taken.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 23 taken 6080 times.
✓ Branch 17 → 18 taken 5125 times.
✓ Branch 17 → 20 taken 5124 times.
✓ Branch 18 → 19 taken 5125 times.
✗ Branch 18 → 20 not taken.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 23 taken 5125 times.
✓ Branch 20 → 21 taken 5124 times.
✗ Branch 20 → 24 not taken.
✓ Branch 21 → 22 taken 5124 times.
✗ Branch 21 → 24 not taken.
✓ Branch 22 → 23 taken 5124 times.
✗ Branch 22 → 24 not taken.
|
16329 | & (i == 2 .and. j == 1 .and. k == 3) .or. & |
| 232 | & (i == 3 .and. j == 2 .and. k == 1)) then | ||
| 233 | levi_civita = -1 | ||
| 234 | else | ||
| 235 | ! error stop "invalid indices for Levi-Civita symbol" | ||
| 236 | levi_civita = 0 | ||
| 237 | end if | ||
| 238 | |||
| 239 | 98223 | end function levi_civita | |
| 240 | |||
| 241 | !> @brief Convert a string to uppercase | ||
| 242 | !> | ||
| 243 | !> @param[in] in The input string to be converted | ||
| 244 | !> | ||
| 245 | !> @return The uppercase version of the input string | ||
| 246 | ✗ | function uppercase(in) result(out) | |
| 247 | implicit none | ||
| 248 | |||
| 249 | character(len=*), intent(in) :: in | ||
| 250 | character(len=len(in)) :: out | ||
| 251 | integer :: i, j | ||
| 252 | |||
| 253 | ✗ | do i = 1, len(in) | |
| 254 | ✗ | j = iachar(in(i:i)) | |
| 255 | ✗ | if (j >= iachar("a") .and. j <= iachar("z")) then | |
| 256 | ✗ | out(i:i) = achar(iachar(in(i:i)) - 32) | |
| 257 | else | ||
| 258 | ✗ | out(i:i) = in(i:i) | |
| 259 | end if | ||
| 260 | end do | ||
| 261 | |||
| 262 | ✗ | end function uppercase | |
| 263 | |||
| 264 | !> @brief Convert a string to lowercase | ||
| 265 | !> | ||
| 266 | !> @param[in] in The input string to be converted | ||
| 267 | !> | ||
| 268 | !> @return The lowercase version of the input string | ||
| 269 | 26 | pure function to_lowercase(in) result(out) | |
| 270 | implicit none | ||
| 271 | |||
| 272 | character(len=*), intent(in) :: in | ||
| 273 | character(len=len(in)) :: out | ||
| 274 | integer :: i, j | ||
| 275 | |||
| 276 |
2/2✓ Branch 3 → 4 taken 126 times.
✓ Branch 3 → 8 taken 26 times.
|
152 | do i = 1, len(in) |
| 277 | 126 | j = iachar(in(i:i)) | |
| 278 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 126 times.
|
152 | if (j >= iachar('A') .and. j <= iachar('Z')) then |
| 279 | ✗ | out(i:i) = achar(j + 32) | |
| 280 | else | ||
| 281 | 126 | out(i:i) = in(i:i) | |
| 282 | end if | ||
| 283 | end do | ||
| 284 | 26 | end function to_lowercase | |
| 285 | |||
| 286 | !> @brief Compute the cumulative sum of an array | ||
| 287 | !> | ||
| 288 | !> @param[out] out The output array containing the cumulative sums | ||
| 289 | !> @param[in] in The input array | ||
| 290 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 319 times.
|
319 | subroutine cumsum(out, in) |
| 291 | implicit none | ||
| 292 | |||
| 293 | real(wp), intent(in) :: in(:) | ||
| 294 | real(wp), intent(out) :: out(size(in) + 1) | ||
| 295 | |||
| 296 | integer :: i | ||
| 297 | |||
| 298 | 319 | out(1) = 0 | |
| 299 |
2/2✓ Branch 5 → 6 taken 6905 times.
✓ Branch 5 → 7 taken 319 times.
|
7224 | do i = 1, size(in) |
| 300 | 7224 | out(i + 1) = out(i) + in(i) | |
| 301 | end do | ||
| 302 | 319 | end subroutine cumsum | |
| 303 | |||
| 304 | !> @brief Convert multi-dimensional Cartesian indices to a linear index (0-based) | ||
| 305 | !> | ||
| 306 | !> @param[in] indices The array of Cartesian indices | ||
| 307 | !> @param[in] dims The dimensions of the multi-dimensional array | ||
| 308 | !> | ||
| 309 | !> @return The corresponding linear index | ||
| 310 |
2/4✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 2554 times.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 2554 times.
|
2554 | pure integer function cartesian_to_linear(indices, dims) result(linear_index) |
| 311 | implicit none | ||
| 312 | |||
| 313 | integer, intent(in) :: indices(:) | ||
| 314 | integer, intent(in) :: dims(:) | ||
| 315 | integer :: i, stride | ||
| 316 | |||
| 317 | linear_index = 0 | ||
| 318 | stride = 1 | ||
| 319 |
2/2✓ Branch 7 → 8 taken 3831 times.
✓ Branch 7 → 9 taken 2554 times.
|
6385 | do i = size(dims), 1, -1 |
| 320 | 3831 | linear_index = linear_index + indices(i) * stride | |
| 321 | 6385 | stride = stride * dims(i) | |
| 322 | end do | ||
| 323 | |||
| 324 | 2554 | end function cartesian_to_linear | |
| 325 | |||
| 326 | !> @brief Initialize the random seed with a fixed value for reproducibility | ||
| 327 | !> | ||
| 328 | !> @param[in] random_seed_value An optional integer | ||
| 329 | !> @note This functionality is normally available by means of the random_init intrinsic, but NVHPC does not support it on a single | ||
| 330 | !> image | ||
| 331 | 21 | subroutine init_random_seed(random_seed_value) | |
| 332 | implicit none | ||
| 333 | integer, optional, intent(in) :: random_seed_value | ||
| 334 | integer :: seed_size | ||
| 335 | 21 | integer, allocatable :: seed(:) | |
| 336 | |||
| 337 | 21 | call random_seed(size=seed_size) | |
| 338 |
7/14✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 21 times.
✓ Branch 5 → 4 taken 21 times.
✗ Branch 5 → 6 not taken.
✓ Branch 6 → 7 taken 21 times.
✗ Branch 6 → 8 not taken.
✓ Branch 8 → 9 taken 21 times.
✗ Branch 8 → 10 not taken.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 21 times.
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 21 times.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 21 times.
|
63 | allocate (seed(seed_size)) |
| 339 |
1/2✓ Branch 16 → 17 taken 21 times.
✗ Branch 16 → 18 not taken.
|
21 | if (present(random_seed_value)) then |
| 340 | ✗ | seed = random_seed_value | |
| 341 | else | ||
| 342 |
2/2✓ Branch 21 → 22 taken 168 times.
✓ Branch 21 → 23 taken 21 times.
|
189 | seed = 42 ! default seed value |
| 343 | end if | ||
| 344 | 21 | call random_seed(put=seed) | |
| 345 |
1/2✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 21 times.
|
21 | deallocate (seed) |
| 346 | 21 | end subroutine init_random_seed | |
| 347 | |||
| 348 | end module m_common | ||
| 349 |