[aerosayan]

I want to write fortran 90 or fortran 2003 code and be able to select between either single/double precision at compile time.


In C++ it is extremely easy to do this using template meta-programming.

template<typename T>
T sum(T a, T b)
{
    return a+b;

}

I can use the sum function for adding any type that supports the + operator.


However I have very little idea on how to implement such a generic function in Fortran.


I can select between single and double precision for basic arithmetic using selected_real_kind method and using real(xp) to define the required precision.

! Define single & double precision
integer, parameter :: sp = selected_real_kind(6, 37)
integer, parameter :: dp = selected_real_kind(15, 307)


! The selected precision level
integer, parameter :: xp = selected_real_kind(precision(1.0_dp))

Using this, I can simply write code that would allow simple arithmetic operations in my preferred precision.

! Constant fractions
real(xp), parameter ::             &
    half        = 0.5_xp,          &
    third       = 1.0_xp / 3.0_xp, &
    fourth      = 1.0_xp / 4.0_xp, &
    fifth       = 1.0_xp / 5.0_xp, &
    sixth       = 1.0_xp / 6.0_xp

This is somewhat good, but I face an issue when I try to call fortran functions and subroutines. Since fortran functions and subroutines are not meta-programmable, we can only pass arguments of valid type. For example, we can only get correct square root for double precision arithmetic if I use dsqrt. When I change the value of xp to selected_real_kind(precision(1.0_sp)), the code obviously no longer works correctly when I use dsqrt.


So, I tried to define generic interfaces to mathematical functions like :

!==-------------------------------------------------------------------------==!
!{{{ > MODULE : NCX_GENERIC_MATH
!==-------------------------------------------------------------------------==!
! - Generic interfaces to mathematical functions
!==----                                                                -----==!
module ncx_generic_math
use ncx_precision

private
public :: sqrt_gn, abs_gn

!==-----------------                                       -----------------==!
! > INTERFACE : SQRT_GN
!==----                                                                -----==!
! - Generic interface for SQRT
!==----                                                                -----==!
interface sqrt_gn

pure function sqrt_dp(x)
use ncx_precision, only : dp
implicit none
real(dp), intent(in) :: x
real(dp)             :: sqrt_dp
end function

pure function sqrt_sp(x)
use ncx_precision, only : sp
implicit none
real(sp), intent(in) :: x
real(sp)             :: sqrt_sp
end function

end interface

Then I could define pure functions to use these interfaces:

pure function sqrt_dp(x)
use ncx_precision, only : dp
implicit none
real(dp), intent(in) :: x
real(dp)             :: sqrt_dp
sqrt_dp = dsqrt(x)
end function

pure function sqrt_sp(x)
use ncx_precision, only : sp
implicit none
real(sp), intent(in) :: x
real(sp)             :: sqrt_sp
sqrt_sp = sqrt(x)
end function

Now, my code compiles and works correctly.

I can use sqrt_gn(9.0_xp) to get the correct square root value for my preferred precision and the code works just fine.



But I have some questions....


QUESTIONS :

+ Is there a better way to do this? <-- MOST IMPORTANT REQUEST

+ Is there any performance loss due to using this indirection through the interface and then to the pure functions?