Fortran function return value polymorphism

[zoziha]

Problem background

Suppose I need to develop a multi-precision math library and define a polymorphic function eye.
Is there any magic in Fortran to make x=eye(2,2) return the value of the eye function according to the type of x?

Change to subroutine?

subroutines have different advantages from functions. Functions can be written in formulas.

Overload equal sign(assignment(=))？

Fortran Class (*) in Function Result

module hello_bob
  interface assignment(=)
    module procedure int_equal_func_class
  end interface

contains

  subroutine int_equal_func_class(a,b)
    integer, intent(out) :: a(:)
    class(*), intent(in) :: b(:)

    select type (b)
      type is (integer)
        a = b
    end select
  end subroutine int_equal_func_class

   function func(a)
     class(*), intent(in) :: a(:)
     class(*), allocatable :: func(:)

     ! No intrinsic assignment supported, also see note about bounds
     allocate(func, source=a)
   end function func

end module hello_bob

program bob
  use hello_bob
  integer i(4)

  i=func([1,2,3,4])
  print*, i

end program bob

Is there any downside to customizing such overloads?
For example, print *, eye(2,2) will be invalid?

Or does Fortran have other magic?

(I’m just interested, and I don’t really need a solution. Thank you! )

[ivanpribec]

There is a proposal open for this:

Function kind specified by optional input

IMO, having this would allow programmers to create and use derived-types (multiprecision integers, higher precision reals, reals supporting auto-differentiation) that are fully compatible with all of the intrinsic functions. Currently, if you a have program that relies on explicit conversions like real(val,sp) (assume val is a double-precision), you cannot replace val with a custom numeric type.

The best workaround I can think of currently is to pass a dummy variable to the procedure.

function eye(n,i) result(res)
  integer, intent(in) :: n
  integer(kind=...) :: i   ! just here to get right output kind
  integer(kind=kind(i)) :: res(n,n)
end function

[veryreverie]

You can do this with deferred evaluation

module m
  implicit none
  
  integer, parameter :: dp = selected_real_kind(15,300)
  
  type :: Eye
    integer :: m
    integer :: n
  end type
  
  interface assignment(=)
    module procedure assign_integers_Eye
    module procedure assign_reals_Eye
  end interface
contains

subroutine assign_integers_Eye(output,input)
  integer, allocatable, intent(out) :: output(:,:)
  type(Eye),            intent(in)  :: input
  
  integer :: i
  
  allocate(output(input%m,input%n))
  output = 0
  do i=1,min(input%m,input%n)
    output(i,i) = 1
  enddo
end subroutine

subroutine assign_reals_Eye(output,input)
  real(dp), allocatable, intent(out) :: output(:,:)
  type(Eye),             intent(in)  :: input
  
  integer :: i
  
  allocate(output(input%m,input%n))
  output = 0
  do i=1,min(input%m,input%n)
    output(i,i) = 1
  enddo
end subroutine
end module

program p
  use m
  implicit none
  
  integer,  allocatable :: int_array(:,:)
  real(dp), allocatable :: real_array(:,:)
  
  integer :: i
  
  int_array = Eye(2,2)
  
  write(*,*)
  do i=1,2
    write(*,*) int_array(i,:)
  enddo
  
  real_array = Eye(2,2)
  write(*,*)
  do i=1,2
    write(*,*) real_array(i,:)
  enddo
end program

The general form for this kind of deferred evaluation is that a function foo(a,b,c) returns a type foo which simply contains a, b and c, and then the assignment bar = foo is overloaded to do what you’d expect the function foo(a,b,c) to do.

I wouldn’t recommend overloading integer foo(:) = class(*) bar(:), as this will conflict with all integer foo(:) assignments, including foo(:) = [1,2,3].

[zoziha]

Your example is not runnable using the existing Fortran syntax, is it?

integer :: sp=4, dp = 8
print *, kind(sp), kind(dp)
   !! kind(4)=4, kind(8)=4

[ivanpribec]

Apologies. Click below for a full example of the workaround I had in mind:

Click here

module eye_mod

  implicit none

  integer, parameter :: sp = kind(1.0e0)
  integer, parameter :: dp = kind(1.0d0)

  interface eye
    module procedure eye_default
    module procedure eye_sp
    module procedure eye_dp
  end interface

contains

  function eye_default(n) result(res)
    integer, intent(in) :: n
    real(sp) :: res(n,n)
    integer :: k
    res = 0
    do k = 1, n
      res(k,k) = 1
    end do
  end function

  function eye_sp(n,i) result(res)
    integer, intent(in) :: n
    real(sp), intent(in) :: i
    real(kind(i)) :: res(n,n)
    integer :: k
    res = 0
    do k = 1, n
      res(k,k) = 1
    end do
  end function

  function eye_dp(n,i) result(res)
    integer, intent(in) :: n
    real(dp), intent(in) :: i
    real(kind(i)) :: res(n,n)
    integer :: k
    res = 0
    do k = 1, n
      res(k,k) = 1
    end do
  end function

end module

The main program would look something like:

program test_eye

  use eye_mod, only: eye, wp => dp

  implicit none

  ! Dummy variable used to select precision
  real(wp), parameter :: dum_wp = 1.0_wp

  integer, parameter :: n = 3
  real(wp) :: I(n,n)
  integer :: k

  I = eye(n,dum_wp)

  do k = 1, n
    print *, I(k,:)
  end do

end program

[zoziha]

Thank you for providing a magic, eye has changed from a function to a type structure, but it achieves its function wonderfully!

Haha, If I have functions zeros(dim1), zeros(dim1,dim2), which have different dimensions, can I still use your method to construct such a polymorphic return value function?

[zoziha]

Thank you very much! Clean and tidy!

[veryreverie]

Haha, yeah, “magic” is the right word. See the C++ Armadillo library for this kind of nonsense writ large (complete with incomprehensible template error messages yelling about “meat glue” for some reason).

You can make zeros work; you can have two types, zeros1d which is constructed by zeros(dim1) and zeros2d which is constructed by zeros(dim1,dim2). Then you just need to overload integer(:) = zeros1d and integer(:,:) = zeros2d.

You can get pretty far with this approach, although I can’t promise that anyone reading your code in the future (probably including future you) will thank you for it!

[zoziha]

Can I achieve the effect of the official function, similar to real(x,8) instead of real(x, 8.0_wp)? The former may be more intuitive and concise. Is there such a method?

[ivanpribec]

@zoziha

Can I achieve the effect of the official function, similar to real(x,8) instead of real(x, 8.0_wp) ? The former may be more intuitive and concise. Is there such a method?

There is no method currently.

The workaround I presented relies on using an actual argument of given kind in place of its integer kind specifier.

Please avoid using integer literals like 4 and 8 as kind specifier at all costs. When you write real(x,8) it is not obvious that 8 is the kind specifier . One improvement is to be explicit with real(x,kind=8). This can still result in non-portability, as the integer kind specifiers are not universal across compilers!

When you write my_func(x,8.0_wp) it is again misleading, as the numerical value has no purpose other than to communicate the kind (e.g. you could pick -555.0_wp with the same end effect). That is why I defined a dummy parameter dum_wp in the example above.

[zoziha]

After practice, I found that it seems that we are too stubborn on Fortran's strong typing attributes, or underestimated the strength of Fortran?

example

program test

    real :: stime, etime
    integer, allocatable :: iA(:,:)
    real(kind=8), allocatable :: rA(:,:)
    complex :: cA(2000,2000)

    call cpu_time(stime)
    iA = ones(2000,2000)
    call cpu_time(etime)
    print "(A,ES10.3)", "etime-stime (seconds): ", etime-stime

    call cpu_time(stime)
    rA = ones(2000,2000)
    call cpu_time(etime)
    print "(A,ES10.3)", "etime-stime (seconds): ", etime-stime
    
    call cpu_time(stime)
    cA = ones(2000,2000)
    call cpu_time(etime)
    print "(A,ES10.3)", "etime-stime (seconds): ", etime-stime

contains

    pure function ones(ndim1, ndim2) result(result)

        integer, intent(in) :: ndim1, ndim2
        integer, allocatable :: result(:,:)
        allocate(result(ndim1, ndim2), source=1)

    end function ones

end program test

The example shows that it is possible to convert integer type to real/complex, and gfortran and ifort run at the same level of speed.
It is still flawed.
For detailed implementation, please see:
ones:

1. docs: stdlib/stdlib_linalg.md at add_exs/ones_and_zero · fortran-fans/stdlib (github.com)
2. code: stdlib/stdlib_linalg_expand.fypp at c6eaa283a8c3e8663b24a2c87ffd7f530c5b46e3 · fortran-fans/stdlib · GitHub

eye code: stdlib/stdlib_linalg.fypp at ebc61cd15278e270988341664856b873ba07c3a3 · fortran-fans/stdlib · GitHub

zeros: …

[nshaffer]

Relying on type numeric conversion can easily lead to unexpected results. For example, assuming ones is defined as in your post:

Real, allocatable :: a(:,:)

A = ones(4,4)/2
Print *, a
End

I expect most users would want to see a bunch of 0.5s printed out, rather than zeros.

(Note I have not actually compiled this program, so I apologize if I have made a silly error)

[zoziha]

Yes, this solution has this disadvantage and some others as well, because it is essentially an integer type of array, which has been converted here, so we should use

A = ones(4,4)/2.0

This is relatively predictable, and we should indicate such type conversions and drawbacks in the existing implementation docs, examples as well.
The situation you said also appears in:

real :: a
a = 1/2 !! a = 0.0

The Fortran committee is relatively friendly to compiler developers, and seems to have added a lot of shackles to language users. I hope that generics and the following proposals (see Function kind specified by optional input · Issue #91 · j3-fortran/fortran_proposals (github.com)) can be passed, and we can change the existing solution.
Let’s let this code implementation stay in the experimental version, which is our relatively better solution now.