Use pointer to refer function or array in Fortran

Question

The problem is to write code that will allow to use the same name of variable that can refer array or function basing on conditions.

In more details: I have real*8 function that takes 4 integer indices. Calculation of such function is quite expensive. The values that return this function are used many millions times. The most obvious solution - to calculate them once and use 4-dimension array instead of function. It saves huge amount of time.

But when the size of the task increases, there is no possibility to store such array in memory. So to be hardware-independent, it is required to turn off storing in the memory and use function instead.

The only one solution that came to me is to use abstract interface with "virtual function" that do nothing but return the values of array by four indices.

First of all, It is necessary to define abstract interface:

abstract interface
   real(kind=rglu) function integral(a,b,c,d)
   use glob, only: iglu,rglu
   integer(kind=iglu), intent(in) :: a,b,c,d
   end function integral
end interface
procedure (integral), pointer :: R=>null()

next, write the function:

real(kind=rglu) function stored_int(a,b,c,d) result(ret)
implicit none
integer(kind=iglu), intent(in) :: a,b,c,d

ret=hR(a,b,c,d); return
end function stored_int

and then we will use it in the following way:

       if (storeIntegrals) then
            do i = 1,N
            do j = 1,N
            do a = 1,N
            do b = 1,N
                hR(i,j,a,b)=my_integral(i,j,a,b)
            enddo
            enddo
            enddo
            enddo
            R=>stored_int
       else
            R=>my_integral
       endif

Function my_integral is that we need to replace with an array. Such approach unfortunately shows very bad performance. Compiled with ifort -O3 -fpp -openmp on four cores it gives twice worse result then the same code but with an array (without virtual function).

Are any other variants to solve this problem?

BTW, there is no point putting return into every function in Fortran. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Feb 2, 2018 at 23:50
I used to put return and if it does not affect performance I will continue writing it. — MuKeP
– MuKeP, Commented Feb 3, 2018 at 22:09

user5713492 · Accepted Answer · 2018-02-03 10:59:49Z

2

One thing you could try is to declare integral, stored_int, and my_integral to be BIND(C) and have their 4 arguments passed by value. This might result in stored_int being invoked a little faster.

Failing this there are still some things that might work for you. You could try creating a user-defined type Tarray that contains hR as component R and type Tfun that contains my_integeral as component R. Then the syntax for accessing an element of hR would be identical to that for invoking function my_integral. You need only maintain one code base, which would be moved to an INCLUDE file. Then you can invoke one or the other via a generic name. Here is such an INCLUDE file:

! sub.i90
subroutine sub1(T1,k,mess)
   implicit none
   type(T) T1
   integer k
   character(*) mess
   write(*,'(a)') mess
   write(*,'(*(g0:1x))') T1%R(k)
end subroutine sub1

And the stuff needed to set up the generic machinery:

! funarray.f90
module funmod
   use ISO_FORTRAN_ENV, only: wp => REAL64
   implicit none
   private
   type, public :: Tfun
      contains
      procedure, NOPASS :: R => fun
   end type Tfun
   contains
      function fun(i) bind(C)
         integer, value :: i
         real(wp) fun
         fun = i ! or whatever
      end function fun
end module funmod

module arraymod
   use ISO_FORTRAN_ENV, only: wp => REAL64
   implicit none
   private
   integer, parameter :: N = 10
   type, public :: Tarray
      real(wp) R(N)
   end type Tarray
end module arraymod

module genfun
   use funmod, only: T => Tfun
   implicit none
   private
   public sub
   interface sub
      module procedure sub1
   end interface sub
   contains
include 'sub.i90'
end module genfun

module genarray
   use arraymod, only: T => Tarray
   implicit none
   private
   public sub
   interface sub
      module procedure sub1
   end interface sub
   contains
include 'sub.i90'
end module genarray

module gencombine
   use genfun
   use genarray
   implicit none
end module gencombine

program gentest
   use funmod
   use arraymod
   use gencombine
   implicit none
   type(Tfun) T1
   type(Tarray) T2
   integer i
   do i = 1, size(T2%R)
      T2%R(i) = T1%R(i)
   end do
   call sub(T1,3,'Invoked for function')
   call sub(T2,4,'Invoked for array')
end program gentest

Output with ifort:

Invoked for function
3.000000000000000
Invoked for array
4.000000000000000

answered Feb 3, 2018 at 10:59

user5713492

9846 silver badges11 bronze badges

Sign up to request clarification or add additional context in comments.

7 Comments

francescalus Over a year ago

Regarding bind(c), why is it necessary to use that for pass-by-value, and what if real(rglu) and integer(igul) aren't C interoperable?

Vladimir F Героям слава Over a year ago

@francescalus Only with bind(c) one gets the true C-like pass-by-value. Maybe in this function-pointer case, when the compiler can't do the usual optimization machinery, it will make some difference. Normally it doesn't. The difference will probably vary between the 32-bit x86 and 64-bit x86 calling conventions.

Vladimir F Героям слава Over a year ago

Not sure about the other solution, I assume runtime switching is necessary. Compile-time switching between a function and an array is trivial and no derived type is necessary.

user5713492 Over a year ago

@francescalus Intel Fortran sort of broke Fortran pass by value except for BIND(C) and some usages of !DEC$ ATTRIBUTES VALUE. With gfortran BIND(C) would not be necessary to enforce true pass by value in this context. The companion processor might be the Fortran compiler itself in which case there would be no question of interoperability of numeric types. Anyhow, one could set rglu = C_DOUBLE and igul=C_INT or igul=C_INT64_T if necessary.

user5713492 Over a year ago

@VladimirF runtime switching is possible even without the derived type and I did have some misgivings about putting it in there after posting, but in defense of my approach I can envision instances where the derived type could lead to faster code. What if the shape of the array hR (they've gotta enable MathJax on this site; comments look so primitive without it) were known at compile time but array or function were a runtime decision? Then 'hR' must be ALLOCATABLE but with derived type the compiler can still hardwire in its shape for possibly more efficient code.

|

Vladimir F Героям слава · Accepted Answer · 2018-02-02 23:49:04Z

1

It is completely impossible to refer with the same pointer to a function and to an array. Even in C data pointers (void *) and function pointers are different. Even in the hardware in CPU's the addresses for data and for code can be implemented differently.

So yes, the wrapper "virtual" function you used seems to be the only way I can see to use the same pointer.

answered Feb 2, 2018 at 23:49

Vladimir F Героям слава

60.7k4 gold badges82 silver badges131 bronze badges

2 Comments

MuKeP Over a year ago

I understand, that there is no such possibility. The title of the question is quite provocative. I just thought that there is some solution obvious for everybody except me.

Vladimir F Героям слава Over a year ago

Please ask what you really want to ask then. That is what we normally do on this site.

Collectives™ on Stack Overflow

Use pointer to refer function or array in Fortran

2 Answers 2

7 Comments

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

7 Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related