Most routines are unlikely to be called with every type possible, and the user can always cast input values with intrinsics like INT and REAL or declare compatible types and make copies and copy them back to original type. Not ideal, and definitely has several types of overhead, as making a copy of data is not free. Even with a pre-processor that supports templating the combinations can get out of
hand but in practise I find that making a few specific routines for the common cases and then using something like the concept in the following case to cast everything else to one of those cases works reasonably well, and I do not think it is unreasonable to have the user handle converting to a common denominator themselves instead of having to accomodate every possible permutation, but I do sometimes wish for a a two-way function that would cast everything to a supported case and then return a specific type back to the original, but never found the lack of that to be a show-stopper.
program demo_anyscalar_to_double
use, intrinsic :: iso_fortran_env, only : int8, int16, int32, int64
use, intrinsic :: iso_fortran_env, only : real32, real64, real128
implicit none
! call same function with many scalar input types
write(*,*)squarei(2_int8)
write(*,*)squarei(2_int16)
write(*,*)squarei(2_int32)
write(*,*)squarei(2_int64)
write(*,*)squarei(2.0_real32)
write(*,*)squarei(2.0_real64)
write(*,*)squarei(2.0_real128)
contains
function squarei(invalue) result (dvalue)
use M_anything, only : anyscalar_to_double
class(*),intent(in) :: invalue
doubleprecision :: invalue_local
doubleprecision :: dvalue
invalue_local=anyscalar_to_double(invalue)
dvalue=invalue_local*invalue_local
end function squarei
pure elemental function anyscalar_to_double(valuein) result(d_out)
use, intrinsic :: iso_fortran_env, only : error_unit !! ,input_unit,output_unit
implicit none
!$@(#) M_anything::anyscalar_to_double(3f): convert integer or real parameter of "any" kind to doubleprecision
class(*),intent(in) :: valuein
doubleprecision :: d_out
doubleprecision,parameter :: big=huge(0.0d0)
select type(valuein)
type is (integer(kind=int8)); d_out=dble(valuein)
type is (integer(kind=int16)); d_out=dble(valuein)
type is (integer(kind=int32)); d_out=dble(valuein)
type is (integer(kind=int64)); d_out=dble(valuein)
type is (real(kind=real32)); d_out=dble(valuein)
type is (real(kind=real64)); d_out=dble(valuein)
Type is (real(kind=real128))
!!if(valuein.gt.big)then
!! write(error_unit,*)'*anyscalar_to_double* value too large ',valuein
!!endif
d_out=dble(valuein)
type is (logical); d_out=merge(0.0d0,1.0d0,valuein)
type is (character(len=*)); read(valuein,*) d_out
class default
d_out=0.0d0
!!stop '*M_anything::anyscalar_to_double: unknown type'
end select
end function anyscalar_to_double
end program demo_anyscalar_to_double
For simplicity I just called anyscalar_to_double(3f) from a function with one parameter, but the idea scales to any number of parameters. I have a collection of the anyto functions that I do use, but always with an eye to what the overhead can be.
Doing the promoting and copying back and forth is ugly, and BLOCK statements can make the code much cleaner looking now-adays by making it possible to declare the temporary variables close to the call, but that has potentially high overhead if the compiler is not smart enough to not allocate and deallocate the scratch variables when written in a heavily used inner loop or something like that. I do wish there was some syntax to say “promote these to the next supported type and if used as output return it back to the original as the same type as the original” but even that is not a perfect solution, and other solutions like some languages just saying there is one “real” type that is usually a big one definitely has draw-backs too. So
unless Fortran gets a new arbitrary precision type I don’t really see a solution that would handle any type being allowed as an INOUT parameter. It would make an interesting read if proposed as a new Fortran feature. I bet I would learn something from the replies.
tostring in Fortran Wiki is so tmewhat related, in that it shows a function that can take just about any scalar types as arguments of from any of one to nine arguments. I use something like that in several messaging and logging utilities where it allows me to create a message almost like WRITE does; and I use the above type of promotion for creating my own operators, where it auto-promotes much like something like A +B +C would where it is standard that A, B, and C can be virtually any numeric type.
Assuming one parameter is always real and the other is always integer and as in your example
local kind names are used (just like the function anything_to_double there is another that converts anything to several different integer types) there would be no code change at all if the values are passed as input to the procedure, and if the conversion was done by the user no code for a passed value. Perhaps I am missing something; plus my point was the techniques are useful for subsets of the problem, but I agree there is no elegant solution for all cases. But in this case instead of having twenty type-specific routines you have one. An example procedure might clear it up.
I am picturing you have a module with routines in it like “M_anything” that exists as a standard library. You write a routine that takes two numeric values of any standard types, and convert the one to REAL and the second one to INTEGER and use them to calculate R**I and return the value, and in this case it being a function and Fortran being able to do explicit casting to the LHS the variable MYVALUE can be any numeric type big enough to hold the result, as an example of the “best case” that is a subcase of what you mentioned as I understood it.
module M_anything
use, intrinsic :: iso_fortran_env, only : int8, int16, int32, int64
use, intrinsic :: iso_fortran_env, only : real32, real64, real128
implicit none
private
public :: anyscalar_to_double, anyscalar_to_int64
contains
pure elemental function anyscalar_to_double(valuein) result(d_out)
use, intrinsic :: iso_fortran_env, only : error_unit !! ,input_unit,output_unit
!$@(#) M_anything::anyscalar_to_double(3f): convert integer or real parameter of "any" kind to doubleprecision
class(*),intent(in) :: valuein
doubleprecision :: d_out
select type(valuein)
type is (integer(kind=int8)); d_out=dble(valuein)
type is (integer(kind=int16)); d_out=dble(valuein)
type is (integer(kind=int32)); d_out=dble(valuein)
type is (integer(kind=int64)); d_out=dble(valuein)
type is (real(kind=real32)); d_out=dble(valuein)
type is (real(kind=real64)); d_out=dble(valuein)
Type is (real(kind=real128))
d_out=dble(valuein)
type is (logical); d_out=merge(0.0d0,1.0d0,valuein)
type is (character(len=*)); read(valuein,*) d_out
class default
d_out=0.0d0
!stop '*M_anything::anyscalar_to_double: unknown type'
end select
end function anyscalar_to_double
impure elemental function anyscalar_to_int64(valuein) result(ii38)
use, intrinsic :: iso_fortran_env, only : error_unit !! ,input_unit,output_unit
class(*),intent(in) :: valuein
integer(kind=int64) :: ii38
integer :: ios
character(len=256) :: message
select type(valuein)
type is (integer(kind=int8)); ii38=int(valuein,kind=int64)
type is (integer(kind=int16)); ii38=int(valuein,kind=int64)
type is (integer(kind=int32)); ii38=valuein
type is (integer(kind=int64)); ii38=valuein
type is (real(kind=real32)); ii38=int(valuein,kind=int64)
type is (real(kind=real64)); ii38=int(valuein,kind=int64)
Type is (real(kind=real128)); ii38=int(valuein,kind=int64)
type is (logical); ii38=merge(0_int64,1_int64,valuein)
type is (character(len=*)) ;
read(valuein,*,iostat=ios,iomsg=message)ii38
if(ios.ne.0)then
write(error_unit,*)'*anyscalar_to_int64* ERROR: '//trim(message)
stop 2
endif
class default
write(error_unit,*)'*anyscalar_to_int64* ERROR: unknown integer type'
stop 3
end select
end function anyscalar_to_int64
end module M_anything
program demo
use, intrinsic :: iso_fortran_env, only : int8, int16, int32, int64
use, intrinsic :: iso_fortran_env, only : real32, real64, real128
! this could be any type big enough to hold the answer
use, intrinsic :: iso_fortran_env, only : ik1=>int32
implicit none
integer(kind=ik1) :: myvalue
! call same function with many scalar input types
myvalue=real_and_int(2.0_real32,4_int8)
write(*,*)myvalue
myvalue=real_and_int(3.0_real32,2_int16)
write(*,*)myvalue
myvalue=real_and_int(4.0_real64,3_int32)
write(*,*)myvalue
myvalue=real_and_int(5.0_real32,7_int64)
write(*,*)myvalue
myvalue=real_and_int(6.0_real32,6_int8)
write(*,*)myvalue
myvalue=real_and_int(7.0_real64,4_int64)
write(*,*)myvalue
myvalue=real_and_int(8.0_real128,2_int32)
write(*,*)myvalue
contains
function real_and_int(invalue1,invalue2) result (dvalue)
use M_anything, only : anyscalar_to_double, anyscalar_to_int64
class(*),intent(in) :: invalue1, invalue2
doubleprecision :: dvalue
dvalue= anyscalar_to_double(invalue1)**anyscalar_to_int64(invalue2)
end function real_and_int
end program demo
So this is a cherry-picked case but shows that at least for some simple cases (and not counting the module) it is possible to write a multiple-parameter procedure in as little as six lines that will take two arguments of practically any numeric type. Not saying it would always be efficient but it is possible. Any you could add invalue3, invalue4, invalue5, … and add little or no new lines, assuming DVALUE was a function of {invalueN}.
If we added this to the program
function lots(in1,in2,in3,in4,in5) result (dvalue)
use M_anything, only : tod=>anyscalar_to_double, toi=>anyscalar_to_int64
class(*),intent(in) :: in1, in2, in3, in4, in5
doubleprecision :: dvalue
dvalue= tod(in1)**toi(in2)+tod(in3)/tod(in4)-tod(in5)
end function lots
it could be called with pretty much any five numeric values the formula could return a valid answer for:
You could call that function like
```fortran
write(*,*)lots(1.0,3,1234,456.0d0,989)
write(*,*)lots(3,4,5,6,7)
you would reasonably want to add checks in the procedure for ranges and types and such, but I do consider it a useful Fortran feature.
@urbanjost , your examples are using kinds that are not guaranteed to be supported. For example, I frequently use a processor that doesn’t support real128, so your examples above won’t even compile, and thus aren’t portable. Also, you’re any_to_* functions only support a limited number of kinds. If a user calls one with, for example, an 80 bit real, your code doesn’t work.
This is my point. With current Fortran, it is not possible to write code that both:
- can compile on any system
- supports any particular floating point (or integer) implementation
You can only choose one. If you choose 1), you can only use real or double precision (or kind(1.0) and kind(1.0d0). If you choose 2), you can use 2 or 3, maybe 4 if you’re feeling lucky, of the kinds from iso_fortran_env and hope they are supported everywhere your users want to run, and correspond to the kinds your users want.
What we really need is the ability to write something like the following:
pure function mean<T>(samples)
requirements<T>
sumable ! where this is appropriately defined somewhere
end requirements
T, intent(in) :: samples(:)
T :: mean
mean = sum(samples) / size(samples)
end function
Yep, there’s just no good solution for this. It’s been a problem for years. And it’s not just reals, it’s ints, logicals, characters, and user-defined types and classes. Hopefully the future generics feature will fix this and not be too verbose. We need a way for every function to be written so that it works with any types with zero run-time overhead. Will we ever get such a thing, I wonder?
The samples address the permutation question more than the type question in that a function like LOTS() can be called with any numeric type supported by the module. For the example DOUBLEPRECISION was chosen as the output type to keep the discussion simpler because the standard does require it.
Especially with today’s machines it is much rarer to need types to reduce memory usage than in the past, but can still be important, as some recent processors attest to. In the past many codes specifically used the smallest type possible; looking through recent codes I see very little use of multiple kinds. So I do not see it as a huge problem. The only solution that would completely eliminate it would be if Fortran were all arbitrary precision that I can see.
In fact adding and eliminating REAL128 support is one of the main reasons I need preprocessors so another possible solution would be for the standard to specifically define the set of supported kinds but I do not see that happening.
Fortran brings this upon itself by being very agnostic about types, hence the need for something like KIND in the language in the first case. This is the penalty paid for the language being flexible but still trying to remain efficient. There are several wonderful arbitrary precision libraries but the fact they are used mostly in niche cases and that they definitely incur a performance cost for high-end users has not made their use ubiquitous by any means. Allowing many kinds also has the penalty of making sure your types match, hence we get here. When many Fortran implementations really only did support a few real types the need for something like an interface block or module was less because a lot of routines only took one type so there was less chance for a mistake.
So I agree it is a problem, but I think it is a reasonable one to pay for higher performance and flexibiity. The biggest thing I hate is having something like REAL128 versions of routines in general libraries that I have to change to compile on an Apple or if using nvfortran or whatnot.
So I do wish there was some way to write a TYPE statement where the compiler would not fail if there were a case statement for a type, but would only fail if the TYPE were used.
In practice if you define all your own kinds by the precision they require in a module and then
create your own kind name instead of using the pre-defined ones you can change to a new type by changing a few lines and recompiling. That is a feature that is very flexible but looking at real-world codes and what other languages provide it looks like if Fortran had just standardized on REAL*4, REAL*8, … (the almost universal non-standard extensions particularly popular pre Fortran-90) it looks like 95% of people would have been happy with that.
I 100% agree with you, that I want things to be performing. In other words, I absolutely do not want runtime polymorphism or runtime type dispatch (for this use case). Rather, I want this resolution to happen at compile time and then I want compilers to inline and optimize things out.
So at least my goal is to figure out generics in a way that are as performing as writing the 20 type-specific routines by hand.
I don’t define separate versions of procedures for single and double precision reals. Instead a module that defines the KIND is USEd throughout the program:
module kind_mod
public :: wp
integer, parameter :: wp = selected_real_kind(15, 307)
end module kind_mod
Is the main problem with this that one may want to use floats of different precisions in a program, even in a single procedure? One could define single, double, and quadruple precision kinds in the module for procedures that use several kinds of float.
If you are writing an application, then it is typically not a problem, you just choose whatever precision you need, typically double and be done with it. You can even use double for most things, and single precision for some hot spots that need to run fast and where the accuracy allows that.
That is the main reason why Fortran does not have the generic feature, because it was not needed that much. Most Fortran codes are (were) applications. And if people provided a library for some physical computation, they would typically use the minimal precision needed to get the job done, and then you would use it like that.
However, if you want to write a library with generic code, say linear algebra, sorting, almost anything in stdlib, … , then you want to provide functions for all accuracies (kinds), because you want to make your library useful in all user’s applications, no matter which precision they end up using. It is this use case that we are designing a good solution for.
I agree, and thought I had noted that there can be a cost to the brevity allowed, but looking at something like the TO_STRING(3f) function mentioned above, that can quite like a list-directed WRITE statement take any argument types class(*) and type is has specified can be worth that cost.
The intrinsics have the advantage that they are built as part of a particular compiler and therefore can be aware of all the allowed types, but even they often require options to be of a matching type to control the permutation problem. So assuming templates/pre-processors and generic interfaces and class(*) do not meet the ideals, what change to the language (ignoring arbitrary precision, which have been mentioned) is a solution?
Automatic promotion, such as in VALUE=INTEGER_VALUE+REAL_VALUE+DOUBLEPRECISION_VALUE is already part of the language but also has costs. Overloading the plus operator to allow CLASS(*) (or using arbitrary precision) can make for nice neat syntax but seems also to not meet the requirements as expanded on. I am wondering how feasible to implement something like a multiple KIND would be. Instead of a KIND specification resulting in matching one KIND, what if I could declare a variable to allow all the kinds satisfying a range?
Even if still limited to a particular type from a coding standpoint if I could say a variable is all the types having a particular specification, and allow * to designate all kinds of that type, it would cover a lot of cases from the coders viewpoint (not so sure compiler writers would think the same). So instead of a template and a preprocessor and specifying explicit types I have to obtain externally I could write something like
subroutine example(a,b,c)
integer(kind=*) :: a
real(kind=*) :: b,c
c=b**a
end subroutine example
And then it would be up to the compiler whether it promoted or internally expanded this; and in most real cases a restriction saying on any given call that all INTEGER parameters and REAL parameters must be of a single type to restrict permutation heck would seem to cover a lot of cases. Fortran leaves a lot up to the compiler, so it would be free to do something like use the equivalent of class(*) at low optimization levels, but something more sophisticated when higher performance levels are required. There would still be trade-offs but it seems that would cover a lot of turf and provide a way to write a very generic routine that works with all supported types intuitively. That still won’t make a request for a 100-digit integer work if the programming environment does not support it, but that is true of the intrinsics themselves. I cannot call an intrinsic with a 128-bit real if that is not a supported type.
well, constants like 3_* seem a little odd, so maybe always having to give a name to the kind
and to get a range, not sure if an array or colon syntax would be better, like in
selected_real_kind(p=4:10)
I would suggest to be pragmatic and start with a very limited set, e.g.
integer, parameter :: rs = IEEE_elected_real_kind(6, 37)
integer, parameter :: rl = IEEE_selected_real_kind(15, 307)
integer, parameter :: is = selected_int_kind(9)
integer, parameter :: il = selected_int_kind(18)
Names are preliminary 
I would guess that this is sufficient for more than 90% of the use cases and it goes definitely beyond what typical C or C++ libraries offer.
People will complain if this selection is too limited and one can then start discussing the pros (support 5% more use cases) and cons (longer compilation times for all others, more complexity) of adding a new kind.
I believe that such a demand-driven approach is superior to “whiteboard engineering”. Start with a minimal viable product and extend based on the feedback of the customers. This is how python evolves (via PEP) and somehow in contrast to the definition of the Fortran standard. I wouldn’t be surprised if more people are interested in support for complex numbers (which have not been discussed at all) than in real128.
If this limited subset is working for stdlib and all/most packages in FPM in an automated way (via preprocessor and coding conventions), the automated extension to more types is theoretically straight forward.
Does anyone have experience with tolerances for varying precision? Can that be handled automatically, i.e. a tolerance of rl_tol = 5.e-15 for double precision follows automatically from a precision rs_tol = 5 e-7 for single precision? It would be a pity if multiple kinds are supported in theory but numerical performance is degenerated for some cases.
I agree with your approach to writing generic code. I have been thinking about this problem for a while and wanted to post a similar question here before seeing this post. I have personally decided to write minimal working libraries and extend them based on the demand. Regarding tolerance, wouldn’t epsilon() be enough to define a kind-agnostic tolerance?
That would be my thought. Even if you didn’t want your tolerance to be exactly epsilon, you could define it in terms of that (or tiny or huge, but those seem a bit less appropriate).
I’ve used things like 1e3_dp * epsilon(1._dp), but I usually use double precision. I don’t have much experience writing code like 1e3_wp * epsilon(1._wp) where wp (working precision) could then be changed to either single or double precision and things would work.
From my pragmatic perspective this would be a sufficient reason for focusing on these four types. Exit with error if il == -1 and give a warning if rs == rl and we have a solution that fits for more than 90% of the users. The remaining users can still complain (or change to a modern processor).
I don’t know if it matters to the discussion, but something to keep in mind is new hardware. There are hardware accelerators that are using different data types than single and double precision. do concurrent can sometimes take advantage of GPU acceleration now, and there are tensor processing units, tensor cores on GPUs that work with half precision, etc. It would be good to keep these new hardware and their optimal precisions/types in mind, in case the auto-parallelization constructs continue to be improved.
MPLAPACK solves the problem of having LAPACK support multiple kinds of reals (beyond single and double precision) by translating Fortran to C++ using FABLE. The current maintainers of LAPACK use little of modern Fortran. It’s still fixed source form, without argument intents, or modules etc. If there is ever to be a modern LAPACK that uses the features of modern Fortran, it will not come from them.
MPLAPACK supports binary64, binary128, FP80 (extended double), MPFR, GMP, and QD libraries (double-double and quad-double). Users can choose MPFR or GMP for arbitrary accurate calculations, double-double or quad-double for fast 32 or 64 decimal calculations.
As far as I can see, MLAPACK is designed to be used with C++. Calling routines in it from Fortran code may not work unless a C (rather than C++) interface is made available.