other/m_common.F90 - GCC Code Coverage Report

Directory:	src/
Coverage:	low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%

	Coverage	Exec / Excl / Total
Lines:	85.2%	52 / 0 / 61
Functions:	88.9%	8 / 0 / 9
Branches:	66.3%	57 / 0 / 86

    other/m_common.F90
    
        Line
        Branch
        Exec
        Source
      
        !> @brief Common module for the BSpline FEEC package
      
        !>
      
        !> This module contains common definitions, constants, and interfaces used throughout the package.
      
        module m_common
      
          private
      
          public :: wp, pi, factorial, levi_civita, set_union, sort, sort_perm, uppercase, to_lower, zero_function, cumsum, prime_factors, balanced_split
      
          public :: user_function_1d_interface, user_function_3d_interface, user_function_3d_vec_interface, user_function_3d_mat_interface
      
          public :: cartesian_to_linear, brent
      
          public :: init_random_seed
      
          public :: linear_solve_3x3_transpose
      
          public :: VERBOSITY_WARN_NEVER, VERBOSITY_WARN_ON_FAILURE, &
      
                    VERBOSITY_WARN_ALWAYS, VERBOSITY_ERROR_ON_FAILURE
      
          !> Single precision real numbers
      
          integer, parameter :: single = selected_real_kind(6)
      
          !> Double precision real numbers
      
          integer, parameter :: double = selected_real_kind(15)
      
        !> Default value which unused if CMake is used
      
        #ifndef CMAKE_BSPLINE_FEEC_PRECISION
      
        #define CMAKE_BSPLINE_FEEC_PRECISION double
      
        #endif
      
          !> The working precision for real-valued numbers
      
          integer, parameter :: wp = CMAKE_BSPLINE_FEEC_PRECISION
      
          !> The value of pi (3.14159...)
      
          real(wp), parameter :: pi = 3.1415926535897932384626433832795028841971_wp
      
          !> Interface for user-defined functions
      
          abstract interface
      
            pure function user_function_1d_interface(x) result(ans)
      
              import wp
      
              implicit none
      
              real(wp), intent(in) :: x
      
              real(wp) :: ans
      
            end function
      
            pure function user_function_3d_interface(x, y, z) result(ans)
      
              import wp
      
              implicit none
      
              real(wp), intent(in) :: x, y, z
      
              real(wp) :: ans
      
            end function
      
            pure function user_function_3d_vec_interface(ind, x, y, z) result(ans)
      
              import wp
      
              implicit none
      
              integer, intent(in) :: ind
      
              real(wp), intent(in) :: x, y, z
      
              real(wp) :: ans
      
            end function
      
            pure function user_function_3d_mat_interface(s1, s2, x, y, z) result(ans)
      
              import wp
      
              implicit none
      
              integer, intent(in) :: s1, s2
      
              real(wp), intent(in) :: x, y, z
      
              real(wp) :: ans
      
            end function
      
          end interface
      
          ! TODO use enums or some other type for verbosity levels
      
          !> Verbosity levels for warnings and errors
      
          integer, parameter :: VERBOSITY_WARN_NEVER = 0
      
          integer, parameter :: VERBOSITY_WARN_ON_FAILURE = 1
      
          integer, parameter :: VERBOSITY_WARN_ALWAYS = 2
      
          integer, parameter :: VERBOSITY_ERROR_ON_FAILURE = 3
      
          ! Module procedure interfaces for submodules
      
          interface
      
            !> @brief Sort an array of integers using the quicksort algorithm
      
            !>
      
            !> @param[inout] array The array of integers to be sorted
      
            !> @param[in] inverse _(optional)_ If true, sort in descending order; otherwise, ascending order (default: false)
      
            module pure subroutine sort_int(array, inverse)
      
              integer, intent(inout) :: array(:)
      
              logical, intent(in), optional :: inverse
      
            end subroutine
      
            !> @brief Sort an array of real numbers using the quicksort algorithm
      
            !>
      
            !> @param[inout] array The array of real numbers to be sorted
      
            !> @param[in] inverse _(optional)_ If true, sort in descending order; otherwise, ascending order (default: false)
      
            !> @param[in] sort_fun _(optional)_ A function to determine the sorting order
      
            module subroutine sort_real(array, inverse, sort_fun)
      
              real(wp), intent(inout) :: array(:)
      
              logical, intent(in), optional :: inverse
      
              procedure(user_function_1d_interface), optional :: sort_fun
      
            end subroutine
      
            !> @brief Compute the union of two sets of integers
      
            !>
      
            !> @param[out] union The resulting union of the two sets
      
            !> @param[in] set1 The first set of integers
      
            !> @param[in] set2 The second set of integers
      
            module subroutine set_union(union, set1, set2)
      
              integer, allocatable, intent(out) :: union(:)
      
              integer, intent(in) :: set1(0:), set2(0:)
      
            end subroutine
      
            !> @brief Compute the prime factorization of an integer
      
            !>
      
            !> @param[in] n The integer to factorize
      
            !> @param[out] factors Array containing the prime factors (with repetition)
      
            module pure subroutine prime_factors(n, factors)
      
              integer, intent(in) :: n
      
              integer, allocatable, intent(out) :: factors(:)
      
            end subroutine
      
            !> @brief Split an integer into two factors that are as close as possible
      
            !>
      
            !> Given an integer N, find n1 and n2 such that N = n1 * n2 and |n1 - n2| is minimized.
      
            !> The result satisfies n1 <= n2.
      
            !>
      
            !> @param[in] n The integer to split
      
            !> @param[out] n1 The smaller factor
      
            !> @param[out] n2 The larger factor
      
            module pure subroutine balanced_split_2(n, n1, n2)
      
              integer, intent(in) :: n
      
              integer, intent(out) :: n1, n2
      
            end subroutine
      
            !> @brief Split an integer into three factors that are as close as possible
      
            !>
      
            !> Given an integer N, find n1, n2, and n3 such that N = n1 * n2 * n3 and the factors are as balanced as possible.
      
            !> The result satisfies n1 <= n2 <= n3.
      
            !>
      
            !> @param[in] n The integer to split
      
            !> @param[out] n1 The smallest factor
      
            !> @param[out] n2 The middle factor
      
            !> @param[out] n3 The largest factor
      
            module pure subroutine balanced_split_3(n, n1, n2, n3)
      
              integer, intent(in) :: n
      
              integer, intent(out) :: n1, n2, n3
      
            end subroutine
      
            !> @brief Find a root of a function using Brent's method
      
            !>
      
            !> @param[in] fun The function for which to find a root
      
            !> @param[in] xa_in The lower bound of the interval
      
            !> @param[in] xb_in The upper bound of the interval
      
            !> @param[in] x_tol The tolerance for convergence
      
            !> @param[in] max_iter The maximum number of iterations to perform
      
            !> @param[in] fa_in _(optional)_ The value of the function at xa_in (if not provided, it will be computed)
      
            !> @param[in] fb_in _(optional)_ The value of the function at xb_in (if not provided, it will be computed)
      
            !> @param[out] success _(optional)_ A logical flag indicating whether the method converged successfully
      
            !>
      
            !> @return The estimated root of the function
      
            module real(wp) function brent(fun, xa_in, xb_in, x_tol, max_iter, fa_in, fb_in, success) result(xsol)
      
              implicit none
      
              procedure(user_function_1d_interface) :: fun
      
              real(wp), intent(in)      :: xa_in, xb_in, x_tol
      
              integer, intent(in)   :: max_iter
      
              real(wp), intent(in), optional :: fa_in, fb_in
      
              logical, intent(out), optional :: success
      
            end function brent
      
          end interface
      
          !> @brief Generic interface for sorting arrays
      
          interface sort
      
            module procedure sort_int, sort_real
      
          end interface sort
      
          interface
      
            !> @brief Compute a sort permutation of an array of real numbers (descending, high to low)
      
            !>
      
            !> @param[in] array The array of real numbers to rank
      
            !> @param[out] perm The permutation indices such that array(perm) is sorted descending
      
            !> @param[in] sort_fun _(optional)_ Applied to each element before comparison
      
            module subroutine sort_perm(array, perm, sort_fun)
      
              real(wp), intent(in) :: array(:)
      
              integer, intent(out) :: perm(:)
      
              procedure(user_function_1d_interface), optional :: sort_fun
      
            end subroutine
      
          end interface
      
          !> @brief Generic interface for balanced splitting of integers
      
          interface balanced_split
      
            module procedure balanced_split_2, balanced_split_3
      
          end interface balanced_split
      
        contains
      
          !> @brief A function that always returns zero
      
          !>
      
          !> @param[in] xp The x-coordinate
      
          !> @param[in] yp The y-coordinate
      
          !> @param[in] zp The z-coordinate
      
          !>
      
          !> @return Always returns 0.0
      
        37888032
          pure real(wp) function zero_function(xp, yp, zp) result(val)
      
            real(wp), intent(in) :: xp, yp, zp
      
            val = 0._wp
      
        37888032
          end function zero_function
      
          !> @brief Compute the factorial of a non-negative integer n
      
          !>
      
          !> @param[in] n The non-negative integer for which to compute the factorial
      
          !>
      
          !> @return The factorial of n
      
        16
          pure integer function factorial(n) result(ans)
      
            implicit none
      
            integer, intent(in) :: n
      
            integer :: i
      
            ans = 1
      
          2/2✓ Branch 3 → 4 taken 120 times.
✓ Branch 3 → 5 taken 16 times.

        136
            do i = 2, n
      
        136
              ans = ans * i
      
            end do
      
        16
          end function
      
          !> @brief Compute the Levi-Civita symbol for three indices
      
          !>
      
          !> The Levi-Civita symbol is defined as:
      
          !>    +1 if (i, j, k) is an even permutation of (1, 2, 3)
      
          !>    -1 if (i, j, k) is an odd permutation of (1, 2, 3)
      
          !>
      
          !> @param[in] i The first index
      
          !> @param[in] j The second index
      
          !> @param[in] k The third index
      
          !>
      
          !> @return The value of the Levi-Civita symbol
      
        74223
          pure integer function levi_civita(i, j, k)
      
            !$acc routine seq
      
            implicit none
      
            integer, intent(in) :: i, j, k
      
          6/6✓ Branch 2 → 3 taken 73974 times.
✓ Branch 2 → 24 taken 249 times.
✓ Branch 3 → 4 taken 49316 times.
✓ Branch 3 → 24 taken 24658 times.
✓ Branch 4 → 5 taken 24658 times.
✓ Branch 4 → 24 taken 24658 times.

        74223
            if (i == j .or. j == k .or. i == k) then
      
              levi_civita = 0
      
            else if ((i == 1 .and. j == 2 .and. k == 3) .or. &
      
          15/18✓ Branch 5 → 6 taken 8204 times.
✓ Branch 5 → 8 taken 16454 times.
✓ Branch 6 → 7 taken 4124 times.
✓ Branch 6 → 8 taken 4080 times.
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 24 taken 4124 times.
✓ Branch 8 → 9 taken 8205 times.
✓ Branch 8 → 11 taken 12329 times.
✓ Branch 9 → 10 taken 4080 times.
✓ Branch 9 → 11 taken 4125 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 24 taken 4080 times.
✓ Branch 11 → 12 taken 8249 times.
✓ Branch 11 → 14 taken 8205 times.
✓ Branch 12 → 13 taken 4125 times.
✓ Branch 12 → 14 taken 4124 times.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 24 taken 4125 times.

        24658
              &      (i == 2 .and. j == 3 .and. k == 1) .or. &
      
              &      (i == 3 .and. j == 1 .and. k == 2)) then
      
              levi_civita = 1
      
            else if ((i == 1 .and. j == 3 .and. k == 2) .or. &
      
          11/18✓ Branch 14 → 15 taken 4080 times.
✓ Branch 14 → 17 taken 8249 times.
✓ Branch 15 → 16 taken 4080 times.
✗ Branch 15 → 17 not taken.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 23 taken 4080 times.
✓ Branch 17 → 18 taken 4125 times.
✓ Branch 17 → 20 taken 4124 times.
✓ Branch 18 → 19 taken 4125 times.
✗ Branch 18 → 20 not taken.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 23 taken 4125 times.
✓ Branch 20 → 21 taken 4124 times.
✗ Branch 20 → 24 not taken.
✓ Branch 21 → 22 taken 4124 times.
✗ Branch 21 → 24 not taken.
✓ Branch 22 → 23 taken 4124 times.
✗ Branch 22 → 24 not taken.

        12329
              &      (i == 2 .and. j == 1 .and. k == 3) .or. &
      
              &      (i == 3 .and. j == 2 .and. k == 1)) then
      
              levi_civita = -1
      
            else
      
              ! error stop "invalid indices for Levi-Civita symbol"
      
              levi_civita = 0
      
            end if
      
        74223
          end function levi_civita
      
          !> @brief Solve the 3x3 linear system A^T x = b using a provided determinant
      
          !>
      
          !> @param[in] A The 3x3 matrix
      
          !> @param[in] detA The determinant of A
      
          !> @param[in] b The right-hand-side vector
      
          !> @param[out] x The solution vector
      
        1002
          pure subroutine linear_solve_3x3_transpose(A, detA, b, x)
      
            !$acc routine seq
      
            implicit none
      
            real(wp), intent(in) :: A(3, 3)
      
            real(wp), intent(in) :: detA
      
            real(wp), intent(in) :: b(3)
      
            real(wp), intent(out) :: x(3)
      
            real(wp) :: c11, c12, c13
      
            real(wp) :: c21, c22, c23
      
            real(wp) :: c31, c32, c33
      
          2/2✓ Branch 2 → 3 taken 1 time.
✓ Branch 2 → 4 taken 1001 times.

        1002
            if (abs(detA) == 0._wp) then
      
        1
              x = 0._wp
      
        1
              return
      
            end if
      
            ! Cofactors C_ij of A; (A^T)^(-1) = C / det(A)
      
        1001
            c11 = A(2, 2) * A(3, 3) - A(2, 3) * A(3, 2)
      
        1001
            c12 = A(2, 3) * A(3, 1) - A(2, 1) * A(3, 3)
      
        1001
            c13 = A(2, 1) * A(3, 2) - A(2, 2) * A(3, 1)
      
        1001
            c21 = A(1, 3) * A(3, 2) - A(1, 2) * A(3, 3)
      
        1001
            c22 = A(1, 1) * A(3, 3) - A(1, 3) * A(3, 1)
      
        1001
            c23 = A(1, 2) * A(3, 1) - A(1, 1) * A(3, 2)
      
        1001
            c31 = A(1, 2) * A(2, 3) - A(1, 3) * A(2, 2)
      
        1001
            c32 = A(1, 3) * A(2, 1) - A(1, 1) * A(2, 3)
      
        1001
            c33 = A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)
      
        1001
            x(1) = (c11 * b(1) + c12 * b(2) + c13 * b(3)) / detA
      
        1001
            x(2) = (c21 * b(1) + c22 * b(2) + c23 * b(3)) / detA
      
        1001
            x(3) = (c31 * b(1) + c32 * b(2) + c33 * b(3)) / detA
      
          end subroutine linear_solve_3x3_transpose
      
          !> @brief Convert a string to uppercase
      
          !>
      
          !> @param[in] in The input string to be converted
      
          !>
      
          !> @return The uppercase version of the input string
      
        ✗
          function uppercase(in) result(out)
      
            implicit none
      
            character(len=*), intent(in) :: in
      
            character(len=len(in)) :: out
      
            integer :: i, j
      
        ✗
            do i = 1, len(in)
      
        ✗
              j = iachar(in(i:i))
      
        ✗
              if (j >= iachar("a") .and. j <= iachar("z")) then
      
        ✗
                out(i:i) = achar(iachar(in(i:i)) - 32)
      
              else
      
        ✗
                out(i:i) = in(i:i)
      
              end if
      
            end do
      
        ✗
          end function uppercase
      
          !> @brief Convert a string to lowercase
      
          !>
      
          !> @param[in] in The input string to be converted
      
          !>
      
          !> @return The lowercase version of the input string
      
        28
          pure function to_lower(in) result(out)
      
            implicit none
      
            character(len=*), intent(in) :: in
      
            character(len=len(in)) :: out
      
            integer :: i, j
      
          2/2✓ Branch 3 → 4 taken 132 times.
✓ Branch 3 → 8 taken 28 times.

        160
            do i = 1, len(in)
      
        132
              j = iachar(in(i:i))
      
          1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 132 times.

        160
              if (j >= iachar('A') .and. j <= iachar('Z')) then
      
        ✗
                out(i:i) = achar(j + 32)
      
              else
      
        132
                out(i:i) = in(i:i)
      
              end if
      
            end do
      
        28
          end function to_lower
      
          !> @brief Compute the cumulative sum of an array
      
          !>
      
          !> @param[out] out The output array containing the cumulative sums
      
          !> @param[in] in The input array
      
          1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 249 times.

        249
          subroutine cumsum(out, in)
      
            implicit none
      
            real(wp), intent(in) :: in(:)
      
            real(wp), intent(out) :: out(size(in) + 1)
      
            integer :: i
      
        249
            out(1) = 0
      
          2/2✓ Branch 5 → 6 taken 5435 times.
✓ Branch 5 → 7 taken 249 times.

        5684
            do i = 1, size(in)
      
        5684
              out(i + 1) = out(i) + in(i)
      
            end do
      
        249
          end subroutine cumsum
      
          !> @brief Convert multi-dimensional Cartesian indices to a linear index (0-based)
      
          !>
      
          !> @param[in] indices The array of Cartesian indices
      
          !> @param[in] dims The dimensions of the multi-dimensional array
      
          !>
      
          !> @return The corresponding linear index
      
          2/4✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 2134 times.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 2134 times.

        2134
          pure integer function cartesian_to_linear(indices, dims) result(linear_index)
      
            implicit none
      
            integer, intent(in) :: indices(:)
      
            integer, intent(in) :: dims(:)
      
            integer :: i, stride
      
            linear_index = 0
      
            stride = 1
      
          2/2✓ Branch 7 → 8 taken 3201 times.
✓ Branch 7 → 9 taken 2134 times.

        5335
            do i = size(dims), 1, -1
      
        3201
              linear_index = linear_index + indices(i) * stride
      
        5335
              stride = stride * dims(i)
      
            end do
      
        2134
          end function cartesian_to_linear
      
          !> @brief Initialize the random seed with a fixed value for reproducibility
      
          !>
      
          !> @param[in] random_seed_value An optional integer
      
          !> @note This functionality is normally available by means of the random_init intrinsic, but NVHPC does not support it on a single
      
          !>       image
      
        22
          subroutine init_random_seed(random_seed_value)
      
            implicit none
      
            integer, optional, intent(in) :: random_seed_value
      
            integer :: seed_size
      
        22
            integer, allocatable :: seed(:)
      
        22
            call random_seed(size=seed_size)
      
          7/14✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 22 times.
✓ Branch 5 → 4 taken 22 times.
✗ Branch 5 → 6 not taken.
✓ Branch 6 → 7 taken 22 times.
✗ Branch 6 → 8 not taken.
✓ Branch 8 → 9 taken 22 times.
✗ Branch 8 → 10 not taken.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 22 times.
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 22 times.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 22 times.

        66
            allocate (seed(seed_size))
      
          1/2✓ Branch 16 → 17 taken 22 times.
✗ Branch 16 → 18 not taken.

        22
            if (present(random_seed_value)) then
      
        ✗
              seed = random_seed_value
      
            else
      
          2/2✓ Branch 21 → 22 taken 176 times.
✓ Branch 21 → 23 taken 22 times.

        198
              seed = 42  ! default seed value
      
            end if
      
        22
            call random_seed(put=seed)
      
          1/2✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 22 times.

        22
            deallocate (seed)
      
        22
          end subroutine init_random_seed
      
        end module m_common

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