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?

Answer 1

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

Answer 2

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.