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

Help
7

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
1 Like