Is it possible to define the product of 2 procedure pointers at runtime?

Tavi007 · February 22, 2024, 11:26am

Hello!

Assume that I have 2 procedure pointers, that can vary at runtime. Both procedures have the same interface. The interface takes in x,y,z real values in and returns also a real value. Now I would like to have procedure pointer to the product. Is this even possible?

I know that this is possible in matlab:

function func = get_linear_volume_basis(id)
    switch id
        case 1
            func = @(x) 1-x(:,1)-x(:,2)-x(:,3);
        case 2
            func = @(x) x(:,1);
        case 3
            func = @(x) x(:,2);
        case 4
            func = @(x) x(:,3);
    end
end

function prod = get_product(i,j)
      func_i = get_linear_volume_basis(i)
      func_j = get_linear_volume_basis(j)
      prod = @(x) func_i(x) * func_j(x)
end

My current fortran code looks like this:

module BasisVolumeFunction_Module
   implicit none

   public :: get_linear_volume_basis_function

   abstract interface
      function integrand(x, y, z) result(val)
         real*8, intent(in) :: x, y, z
         real*8 :: val
      end function integrand
   end interface

   private
contains

   function get_product(i,j) result(ptr)
      integer :: i,j
      procedure(integrand), pointer :: ptr
      procedure(integrand), pointer :: func_i, func_j
      func_i = get_linear_volume_basis_function(i)
      func_j = get_linear_volume_basis_function(j)

      ! what to do here?
   end function get_product

   function get_linear_volume_basis_function(id) result(ptr)
      integer :: id
      procedure(integrand), pointer :: ptr
      select case (id)
      case (1)
         ptr => linear_volume_basis_1
      case (2)
         ptr => linear_volume_basis_2
      case (3)
         ptr => linear_volume_basis_3
      case (4)
         ptr => linear_volume_basis_4
      case default
         print *, 'Error: incorrect id for linear volume basis used:', id
      end select
   end function get_linear_volume_basis_function

   function linear_volume_basis_1(x, y, z) result(val)
      real*8, intent(in) :: x, y, z
      real*8 :: val
      val = 1.0d0 - x - y - z
   end function linear_volume_basis_1

   function linear_volume_basis_2(x, y, z) result(val)
      real*8, intent(in) :: x, y, z
      real*8 :: val
      val = x
   end function linear_volume_basis_2

   function linear_volume_basis_3(x, y, z) result(val)
      real*8, intent(in) :: x, y, z
      real*8 :: val
      val = y
   end function linear_volume_basis_3

   function linear_volume_basis_4(x, y, z) result(val)
      real*8, intent(in) :: x, y, z
      real*8 :: val
      val = z
   end function linear_volume_basis_4

end module BasisVolumeFunction_Module

FedericoPerini · February 22, 2024, 12:05pm

The product interface will look something like this:

function get_product(proca,procb,x,y,z) result(val)
   procedure(integrand) :: proca,procb
   real*8, intent(in) :: x,y,z
   real*8 :: val
   val = proca(x,y,z)*procb(x,y,z)
end function

this is independent of what you have associated proca and procb with. The best approach to store those pointers if you’re going to deal with many options, is within a derived type.

Tavi007 · February 22, 2024, 12:32pm

This is not entirely what i want, i think. your get_product function does not fit the integrand interface, but that is what i acutally want returned from my get_product function.

Arjen · February 22, 2024, 12:45pm

You could do something like:

function get_result( this, x, y, z) 
     if ( .not. associated(this%second) ) then
           get_result = this%first(x,y,z)
     else 
           get_result = this%first(x,y,z) * this%second(x,y,z)
     endif
end function

(leaving out a multitude of details)

ivanpribec · February 22, 2024, 12:54pm

Are i and j dummy arguments really runtime variables, or will they be filled with integer literals?

What seems to be the simpler option is to simply supply the combinations of callbacks (excluding the symmetric values):

    function bxx(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = x*x
    end function

    function bxy(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = x*y
    end function

    function bxz(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = x*z
    end function

    function byy(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = y*y
    end function

    function byz(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = y*z
    end function

    function bzz(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = z*z
    end function

    function bxq(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = x*(1 - x - y- z)
    end function

    function byq(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = y*(1 - x - y- z)
    end function

    function bzq(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = z*(1 - x - y- z)
    end function

    function bqq(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = (1 - x - y- z)**2
    end function

Tavi007 · February 22, 2024, 12:58pm

I would like to see these details. Because i still don’t get how your function ‘get_results’ fits this interface

   abstract interface
      function integrand(x, y, z) result(val)
         real*8, intent(in) :: x, y, z
         real*8 :: val
      end function integrand
   end interface

To give you some context, why i want a function that fits this interface exactly. I want to integrate any function, that fits this interface with various qudrature rules. For this i have

   abstract interface
      function integrate(func) result(val)
         import integrand
         procedure(integrand) :: func
         real*8 :: val
      end function integrate
   end interface

   type, public :: Quadrature
      real*8, allocatable :: weights(:)
      real*8, allocatable :: points(:, :)
   contains
      procedure :: integrate => integrate_volume
   end type Quadrature

   type(Quadrature), public :: rule_1, rule_2, rule_3 ! will be initialized later
contains

   function integrate_volume(this, func) result(val)
      class(Quadrature) :: this
      procedure(integrand) :: func
      real*8 :: val
      integer :: i

      val = 0.0d0
      do i = 1, size(this%weights,1)
         val = val + this%weights(i)*func(this%points(i, 1), &
                                          this%points(i, 2), &
                                          this%points(i, 3))
      end do
   end function integrate_volume
end module QuadratureModule

@ivanpribec
I’ll have a lot more basis function later (as well as their gradient function) and I will have to integrate each product combination. So yes, i could hardcode all of them, but i would rather not.

ivanpribec · February 22, 2024, 1:02pm

The other option I can suggest is to use a dispatch table:

type :: pp
    procedure(integrand), pointer :: f => null()
end type

type(pp), protected :: integrand_table(4,4)

contains

    subroutine init_integrand_table()

        integrand_table(1,1)%f => bxx
        integrand_table(1,2)%f => bxy
        integrand_table(1,3)%f => bxz
        integrand_table(1,4)%f => bxq
        integrand_table(2,1)%f => bxy   ! byx
        integrand_table(2,2)%f => byy 
        ! ...

    end subroutine

For the gradients you could just have another table in the table.

Click here to see the full example

module basis_functions
implicit none
public

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

abstract interface
    function integrand(x,y,z) result(f)
        import dp
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
    end function
end interface

type :: pp
    procedure(integrand), pointer, nopass :: f => null()
end type

type(pp), protected :: basis_table(4,4)

contains

    subroutine init_basis_table

        basis_table(1,1)%f => bxx
        basis_table(2,1)%f => bxy
        basis_table(3,1)%f => bxz
        basis_table(4,1)%f => bxq

        basis_table(1,2)%f => bxy
        basis_table(2,2)%f => byy
        basis_table(3,2)%f => byz
        basis_table(4,2)%f => byq
        
        basis_table(1,3)%f => bxz
        basis_table(2,3)%f => byz
        basis_table(3,3)%f => bzz
        basis_table(4,3)%f => bzq
        
        basis_table(1,4)%f => bxq
        basis_table(2,4)%f => byq
        basis_table(3,4)%f => bzq
        basis_table(4,4)%f => bqq

    end subroutine


    function bxx(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = x*x
    end function

    function bxy(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = x*y
    end function

    function bxz(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = x*z
    end function

    function byy(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = y*y
    end function

    function byz(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = y*z
    end function

    function bzz(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = z*z
    end function

    function bxq(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = x*(1 - x - y- z)
    end function

    function byq(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = y*(1 - x - y- z)
    end function

    function bzq(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = z*(1 - x - y- z)
    end function

    function bqq(x,y,z) result(f)
        real(dp), intent(in) :: x, y, z
        real(dp) :: f
        f = (1 - x - y- z)**2
    end function

end module

program use_table
use basis_functions, table => basis_table
implicit none

real(dp) :: x, y, z
integer :: i, j

call init_basis_table

x = 1.0_dp
y = 0.5_dp
z = 0.3_dp

do i = 1, 4
    do j = 1, 4
        print *, table(i,j)%f(x,y,z)
    end do
end do

end program

Arjen · February 22, 2024, 1:03pm

I can provide the details, bit it will not be right away

everythingfunctional · February 22, 2024, 2:26pm

Fortran doesn’t have what are called closures. Local variables disappear when you leave the scope in which they’re defined (except for saved variables and module variables (which inherently have the save attribute. So

Disappear after the get_product function returns. The simplest change to the code you have is to move that declaration to just before the contains statement of the module, add a procedure

function product_of_basis_functions(x, y, z) result(prod)
    real(dp), intent(in) :: x, y, z
    real(dp) :: prod
    prod = func_i(x, y, z) * func_j(x, y, z)
end function

And then the return from get_product is ptr => product_of_basis_functions.

Note, this means the func_i and func_j are shared between all of the pointers returned from get_product (i.e. subsequent calls to get_product will change func_i and func_j for all calls to product_of_basis_functions). If you want to have separate func_i and func_j for every call to get_product, things get much more complicated. Doable, but complicated.

Arjen · February 22, 2024, 2:52pm

What I am thinking of is less a closure type solution than a sort of composition, but I will try to work it tonight

eelis · February 22, 2024, 3:04pm

Thank you! I was wondering how to do this sort of thing in one of my hobby projects. But last time I tried, I ran into some strange runtime errors. Your example seems to run as intended.

ivanpribec · February 22, 2024, 3:06pm

Here is a second path:


module composite_function_mod
implicit none
private

public :: dp, integrand
public :: product_functor

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

abstract interface
    function integrand(x) result(y)
        import dp
        real(dp), intent(in) :: x(:,:)
        real(dp), allocatable :: y(:)
    end function
end interface

type :: product_functor
    private
    procedure(integrand), pointer, nopass :: a => null(), b => null()
contains
    procedure, public :: eval => eval_functor
end type

interface product_functor
    module procedure product_ij
end interface

contains

    function product_ij(i,j) result(this)
        integer, intent(in) :: i, j
        type(product_functor) :: this

        call set_basis(i,this%a)
        call set_basis(j,this%b)

    contains

        subroutine set_basis(i,f)
            integer, intent(in) :: i
            procedure(integrand), pointer :: f
            select case(i)
            case(1)
                f => b1
            case(2)
                f => b2
            case(3)
                f => b3
            case(4)
                f => b4
            case default
                error stop "invalid index"
            end select
        end subroutine

    end function

    function eval_functor(this,x) result(y)
        class(product_functor), intent(in) :: this
        real(dp), intent(in) :: x(:,:)
        real(dp), allocatable :: y(:)
        y = this%a(x) * this%b(x)
    end function

    function b1(x) result(f)
    real(dp), intent(in) :: x(:,:)
    real(dp), allocatable :: f(:)
    f = 1 - x(:,1) - x(:,2) - x(:,3)
    end function

    function b2(x) result(f)
    real(dp), intent(in) :: x(:,:)
    real(dp), allocatable :: f(:)
    f = x(:,1)
    end function

    function b3(x) result(f)
    real(dp), intent(in) :: x(:,:)
    real(dp), allocatable :: f(:)
    f = x(:,2)
    end function

    function b4(x) result(f)
    real(dp), intent(in) :: x(:,:)
    real(dp), allocatable :: f(:)
    f = x(:,3)
    end function

end module


program test

    use composite_function_mod
    implicit none

    real(dp) :: result(4,4)
    type(product_functor) :: pf

    integer :: i, j

    ! Integration rule
    real(dp), allocatable :: w(:), knots(:,:)

    ! Mock integration rule
    allocate(w(1000),knots(1000,3))
    call random_number(knots)
    w = 1.0_dp/1000


    do i = 1, 4
        do j = 1, 4
            pf = product_functor(i,j)
            result(i,j) = dot_product(w, pf%eval(knots))
        end do
    end do

    do i = 1, 4
        write(*,*) (result(i,j),j=1,4)
    end do

end program

themos · February 22, 2024, 3:07pm

Personal take: The essence of Fortran is that all executable code is present at program start. Fortran will not compose executable code at runtime, and if it ever does it shouldn’t be called Fortran anymore. In Fortran, data comes and goes but code is written in stone.

Tavi007 · February 22, 2024, 3:41pm

So what i get from this thread is, that it is not directly possible, but rather need some clever type bounded procedures. I like the approach from @ivanpribec in his composite_function_mod, so i will mark this one as the solution. For my own usecase, where i want to integrate various compositions (not only the product, but also vector-product and also chained functions), i will have to adjust his code, but that shouldn’t be the problem.

Beliavsky · February 22, 2024, 3:47pm

An exception that comes to mind is that it can be useful to set a format string at run time. It would be limiting if format strings were required to be character literals or named constants, and even worse if format strings were removed and one could only use FORMAT statements.

Nocaster60 · February 22, 2024, 5:00pm

You can use type bound procedures and that’s probably the way to do it but you can use separate procedures if there aren’t too many that you want at the same time - something along the lines of the following. This has two but you could go to quite a few before it gets too tedious.

module BasisVolumeFunction_Module
   implicit none

   public :: get_linear_volume_basis_function
   
   abstract interface
      function integrand(x, y, z) result(val)
         real*8, intent(in) :: x, y, z
         real*8 :: val
      end function integrand
   end interface
   
   type integrand_ptr_t
      procedure(integrand), nopass, pointer :: integrand_ptr
   end type
   
   type(integrand_ptr_t) :: func_i_ptrs(2),func_j_ptrs(2),product_ptrs(2)
   
   private
contains

   subroutine Init
   product_ptrs(1)%integrand_ptr => ProductProc1
   product_ptrs(2)%integrand_ptr => ProductProc2
   end subroutine Init
   

   function get_product(i,j,k) result(ptr)
      integer :: i,j,k
      procedure(integrand), pointer :: ptr
      procedure(integrand), pointer :: func_i, func_j
      func_i_ptrs(k)%integrand_ptr => get_linear_volume_basis_function(i)
      func_j_ptrs(k)%integrand_ptr => get_linear_volume_basis_function(j)
      
      ptr => product_ptrs(k)%integrand_ptr
   end function get_product
   
   
   function ProductProc1(x,y,z) result (val)
   real*8, intent(in) :: x, y, z
   real*8 :: val
   val = func_i_ptrs(1)%integrand_ptr(x,y,z)*func_j_ptrs(1)%integrand_ptr(x,y,z)
   end function ProductProc1
   

   function ProductProc2(x,y,z) result (val)
   real*8, intent(in) :: x, y, z
   real*8 :: val
   val = func_i_ptrs(2)%integrand_ptr(x,y,z)*func_j_ptrs(2)%integrand_ptr(x,y,z)
   end function ProductProc2
   

   function get_linear_volume_basis_function(id) result(ptr)
      integer :: id
      procedure(integrand), pointer :: ptr
      select case (id)
      case (1)
         ptr => linear_volume_basis_1
      case (2)
         ptr => linear_volume_basis_2
      case (3)
         ptr => linear_volume_basis_3
      case (4)
         ptr => linear_volume_basis_4
      case default
         print *, 'Error: incorrect id for linear volume basis used:', id
      end select
   end function get_linear_volume_basis_function

   function linear_volume_basis_1(x, y, z) result(val)
      real*8, intent(in) :: x, y, z
      real*8 :: val
      val = 1.0d0 - x - y - z
   end function linear_volume_basis_1

   function linear_volume_basis_2(x, y, z) result(val)
      real*8, intent(in) :: x, y, z
      real*8 :: val
      val = x
   end function linear_volume_basis_2

   function linear_volume_basis_3(x, y, z) result(val)
      real*8, intent(in) :: x, y, z
      real*8 :: val
      val = y
   end function linear_volume_basis_3

   function linear_volume_basis_4(x, y, z) result(val)
      real*8, intent(in) :: x, y, z
      real*8 :: val
      val = z
   end function linear_volume_basis_4

end module BasisVolumeFunction_Module

themos · February 22, 2024, 5:06pm

The output format specification is encoded in a character variable and is processed as data by the compiler’s runtime routine which is fixed in stone. But you are right in the sense that there is a small format interpreter embedded in the runtime.

For the OP, the “solution” is to write a small interpreter for whatever combinations are anticipated.

Arjen · February 22, 2024, 6:11pm

I have written a small illustration of what I meant. See the attached file.
combine.f90 (3.1 KB)

This can even be extended to full-scale expressions. And an alternative might be the package feq-parse, see Feq-parse v2.2.0 Release - Custom Function Support

Topic		Replies	Views
Procedure statement	6	776	March 19, 2021
Trying to understand why my pointer is 'losing' its value Help	20	556	November 30, 2023
Module procedures and abstract interfaces Help	8	1813	April 6, 2022
Selection of correct procedure/function at runtime (F2008)	8	794	September 11, 2021
Function `associated` bug for Intel Fortran compiler? Help	3	711	June 14, 2022

Is it possible to define the product of 2 procedure pointers at runtime?

Related Topics