[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?

[sbaffini]

While Fortran doesn't have, for some good reasons, the same form of polymorphism used by C++, you are clearly misinterpreting how selected_xxx_kind types are meant to be used. Once you have the xp kind you can have on input to any function/subroutine that knows it something like:

REAL(xp), INTENT(IN) :: input

Because xp is just a placeholder for some already premapped type. Of course, this only works for basic types and not, say, for derived ones (still, they can be parameterized)

Also, you should not, in any case, use something like dsqrt. sqrt would have worked just fine independently from the declaration issue

To give you an example:

MODULE mysqrt_m
USE myprecision_m !where xp is defined
FUNCTION mysqrt(input) RESULT(res)
REAL(xp), INTENT(IN) :: input
REAL(xp) :: res
res = sqrt(input)
ENDFUNCTION mysqrt
ENDMODULE mysqrt_m

[agd]

Why not just use the compiler flag to specify precision? Most FORTRAN compilers have a flag (For intel you can use -r8). Of course, if you are building a module that you want to be portable and usable for either single or double precision then the correct answer is to create single and double precision routines and provide the explicit interface. With regard to performance issues, I have never seen any significant problems, even dealing with complex-valued functions and routines. The interface allows the compiler to know ahead of time which routine to use based on the data type being sent. So there is no type conversion hit to performance (or so I've been told).

[praveen]

agd gives good suggestion. I have always done it like this. Makes the code more readable also. Declare floats as "real :: x", set values like

x = 1.0

Use sqrt (not dsqrt), etc.

Then while compiling, say with gfortran, use -fdefault-real-8 if you want double precision version.

[sbaffini]

Unless there is a very specific business case, honestly, I would kind avoid approaches based on compiler flags. If switching precision is what you want, selected_xxx_kind are there also for this, and you should always add _xp to the reals so defined.

In my opinion code should always speak for itself.

Of course, if you still need to write subroutines for multiple types, you still need an interface etc., but that is a different issue

[aerosayan]

Originally Posted by sbaffini

Unless there is a very specific business case, honestly, I would kind avoid approaches based on compiler flags. If switching precision is what you want, selected_xxx_kind are there also for this, and you should always add _xp to the reals so defined.

In my opinion code should always speak for itself.

Of course, if you still need to write subroutines for multiple types, you still need an interface etc., but that is a different issue

Writing the interface becomes tedious due to the amount of repeated code. I found out a trick to remove repetition. We can #define MYTYPE and #include a template file which uses the MYTYPE macro. Then we can #undef MYTYPE for safety.


This is explained in this stackoverflow answer : https://stackoverflow.com/a/24070364

[sbaffini]

Originally Posted by aerosayan

Writing the interface becomes tedious due to the amount of repeated code. I found out a trick to remove repetition. We can #define MYTYPE and #include a template file which uses the MYTYPE macro. Then we can #undef MYTYPE for safety.


This is explained in this stackoverflow answer : https://stackoverflow.com/a/24070364

Honestly, it really depends from what you actually want to do. It isn't clear if you want to be able to compile your code in both single and double precision versions or if you want to write a set of functions and/or subroutines for both.

In the first case it is possible to use both the compiler flags and the selected kinds approaches. Still, as I wrote, I would prefer the latter.

In the second case, of course, there is no chance to use the two previous approaches and you need a different route. The tedious, formal one is the interface (which is a non solution, because you actually write a lot of identical code + boiler plate stuff for the interface).

The meta-language approach is another one. Personally, I am an integralist of formalisms, and the fact that you are using some meta language to solve an issue with your selected language, well... I know it is used a lot, but I feel pain every time I see such solutions.

Honestly, if a library is what you want to write, and you don't want to use Fortran interfaces and stuff, there is clear indication that you have selected the wrong language in the first place, at least for now (may change in the future, who knows).

I don't write code for others to be used in a library, so I might be skewed, but whenever I find myself in the need of the same subroutine for different types (and honestly, this only happened for linked lists and sort/search stuff) I first take the chance to analyze the thing for every different data type, in order to understand if I can rearrange stuff differently, reuse memory in some cases, etc. Then I just put the different version in a module and keep the names different and clear. Later in time, when everything is stabilized, I decide if the boilerplate code is actually worth the effort just for the sake of using the same name for calling the function/subroutine on different types... the answer is typically not, but there is certainly bias due to my use case.

[el_mojito]

Originally Posted by aerosayan

QUESTIONS :

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

Yes. You are doing this waaaaaa(...)aay to complicated.
Use modern Fortran and just use sqrt(), which will take either a single or double precision real and return a value of the same kind.





Second point:
Don't use selected_real_kind(). Use the REAL32/REAL64-Types from the ISO_FORTRAN_ENV module

use, intrinsic :: ISO_FORTRAN_ENV, only : INT8, INT16, INT32, INT64, &       ! Standardized Fortran 2008 kind definitions
                                          REAL32, REAL64, REAL128
integer, parameter :: rk = REAL64

real(rk) :: my_real

[Gerry Kan]

Dear Sayan:

If you just want to switch between single and double precision float points during compile time, you can declare them as REAL(KIND=4)'s (single precision) in the code, and then use the gfortran compiler switch -freal-4-real-8 to force them into REAL(KIND=8)'s (double precision). You can even convert them into "long double" with the -freal-4-real-16 switch if you want.

There used to be an -r8 switch which will force all REAL's into DOUBLE PRECISION's in g77, but it seems that it's not in gfortran anymore.

You might want to check for portability first, but I suspect there should be compiler options in other Fortran products that are equivalent.

Gerry.

P.S. - here is a list of gfortran compiler options

[sbaffini]

Originally Posted by el_mojito

Second point:
Don't use selected_real_kind(). Use the REAL32/REAL64-Types from the ISO_FORTRAN_ENV module

use, intrinsic :: ISO_FORTRAN_ENV, only : INT8, INT16, INT32, INT64, &       ! Standardized Fortran 2008 kind definitions
                                          REAL32, REAL64, REAL128
integer, parameter :: rk = REAL64

real(rk) :: my_real

Without context, this is just bad advice. None of them are actually required to exist by the standard.

[el_mojito]

Originally Posted by sbaffini

Without context, this is just bad advice. None of them are actually required to exist by the standard.

Um, yes, they are?


Page 400:

https://j3-fortran.org/doc/year/10/10-007r1.pdf

[sbaffini]

Originally Posted by el_mojito

Um, yes, they are?


Page 400:

https://j3-fortran.org/doc/year/10/10-007r1.pdf

Ok, lazy me... none of them are actually required to return usable values for any given platform.

If the processor supports no kind of a particular size, that constant shall be equal to −2 if the processor supports a kind with larger size and −1 otherwise.

[sbaffini]

Also, this https://stevelionel.com/drfortran/20...kes-all-kinds/