Abstract interfaces with deferred shape arrays

Help
1

I’m trying to build some simple numerical functions as a way of getting used to Fortran, and I’m starting with derivatives. I want to write a function that returns me the Jacobian of any function that I pass, independent of the number of inputs and outputs. I’ve built the following module:

module solvers
use kinds, only: wp => dp
implicit none

abstract interface
    function func(x) result (y)
    ! y = f(x), scalar function
        import :: wp
        real(wp), intent(in) :: x
        real(wp) :: y
    end function func

    function func2(x) result (y)
    ! function with arbitrary number of inputs and outputs
    ! y(1) = f(x1, x2, ...), y(2) = g(x1, x2, ...), etc...
        import :: wp
        real(wp), intent(in), dimension(:) :: x
        real(wp), dimension(:) :: y
    end function func2
end interface

contains

function diff(x, F, opt_h, opt_method) result (J)
    procedure(func1) :: F
! code to calculate derivatives...
end function diff

This code doesn’t compile because of func2 with the error: Error: Array 'y' at (1) cannot have a deferred shape. It does compile if I declare y as allocatable on the abstract interface. But then, all functions that I pass to the function diff need to have an allocatable output, otherwise i get the error Error: Interface mismatch in dummy procedure 'f' at (1): ALLOCATABLE attribute mismatch in function result.

Ideally I’d like my generic diff function to work with whatever function I pass to it, regardless if its result is allocatable or not. Is this possible?

Additionally, I’m defining the diff function like this:

    function diff(x, F, opt_h, opt_method) result(J)
        real(wp), intent(in), dimension(:)            :: x
        real(wp), intent(in), optional              :: opt_h
        procedure(func2) :: F
        character(len=20), optional     :: opt_method
        real(wp)                        :: h
        real(wp), allocatable, dimension(:)          :: F_x
        integer                         :: i, n_vars, n_funcs
        character(len=20)               :: method
        real, allocatable, dimension(:,:)            :: J

        n_vars = size(x)
        F_x = F(x)
        n_funcs = size(F_x)
        ! code to calculate derivatives...
    end function diff

Basically, I’m calling F to find out the size of the array F returns. This works in this situation, but is there any way of evaluating the size of a function return array without the function call?

Thanks in advance!

1 Like
2

Do you have a small working example of code that compiles and produces that error ?
Maybe you could give func2() another argument that tells it what size y(:) should be, like

 function func2(x, ny) result (y)
        import :: wp
        real(wp), intent(in), dimension(:) :: x
        integer, intent(in) :: ny
        real(wp), dimension(ny) :: y
end function func2
1 Like
3

@DCLopes ,

You seek Generics facility in Fortran, but unfortunately that is a couple of decades away, if not much longer or never.

In the meantime, have you looked at the parameterized derived type (PDT) feature, introduced starting with Fortran 2003, over 20 years ago, but which is still very buggy or dysfunctional in most Fortran processor (compiler) implementations :sob:

It is unclear to me yet from the original post the design you seek - take a look at the following, is this along the lines what you would prefer?

module diff_m
   type, abstract :: diff_t(n)
      integer, len :: n = 2
   contains
      procedure(IFunc), deferred :: Func
      procedure(IJacobian), deferred :: Jacobian
   end type
   abstract interface
      function Ifunc( this, x ) result( f )
         import :: diff_t
         ! Argument list
         class(diff_t(n=*)), intent(in) :: this
         real, intent(in) :: x( this%n )
         ! Function result
         real :: f( this%n )
      end function 
      function IJacobian( this, x ) result( df )
         import :: diff_t
         ! Argument list
         class(diff_t(n=*)), intent(in) :: this
         real, intent(in) :: x( this%n )
         ! Function result
         real :: df( this%n, this%n )
      end function 
   end interface
end module
module ex_m
   use diff_m, only : diff_t, Ifunc, IJacobian
   type, extends(diff_t) :: ex_t
   contains
      procedure :: Func => ex_func
      procedure :: Jacobian => ex_Jacobian
   end type
contains
   function ex_func( this, x ) result( f )
      ! Argument list
      class(ex_t(n=*)), intent(in) :: this
      real, intent(in) :: x( this%n )
      ! Function result
      real :: f( this%n )
      f(1) = x(1)*x(2)
      f(2) = 5.0*x(1) + sin( x(2) )
   end function 
   function ex_Jacobian( this, x ) result( df )
      ! Argument list
      class(ex_t(n=*)), intent(in) :: this
      real, intent(in) :: x( this%n )
      ! Function result
      real :: df( this%n, this%n )
      df(1,1) = 2.0*x(1)*x(2)
      df(2,1) = 5.0
      df(1,2) = x(1)*x(2)
      df(2,2) = cos( x(2) )
   end function 
end module
   use ex_m, only : ex_t
   type(ex_t) :: ex
   real :: x( ex%n ), y( ex%n), dy( ex%n, ex%n ) 
   x = [ 1.0, 2.0 ]
   y = ex%Func( x )
   dy = ex%Jacobian( x )
   print *, "f = ", y
   print *, "Jacobian = ", dy
end

You can try NAG Fortran latest, if you have access, or Intel Fortran (IFORT) with above and tinker with the design:

C:\temp>ifort /free /standard-semantics /traceback /Qdiag-disable:10448 p.f
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.12.0 Build 20240222_000000
Copyright (C) 1985-2024 Intel Corporation.  All rights reserved.

Microsoft (R) Incremental Linker Version 14.36.32537.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:p.exe
-subsystem:console
-incremental:no
p.obj

C:\temp>p.exe
 f = 2.000000 5.909297
 Jacobian = 4.000000 5.000000 2.000000 -0.4161468
4

May be something like this could help you:

module example
use iso_fortran_env, only: wp => real64
implicit none
abstract interface
    function fun(x) result(y)
        import :: wp
        real(wp), intent(in) :: x(:)
        real(wp) :: y(size(x))
    end function
end interface

contains

function eval(f,x) result(y)
    procedure(fun) :: f
    real(wp), intent(in) :: x(:)
    real(wp) :: y(size(x))
    y = f(x)
end function

end module

program main
use example
use iso_fortran_env, only: wp => real64

real(wp) :: x(6)
real(wp) :: y(6)

x = 2._wp
y = eval(myfun,x)
print *, y !> will print 4 4 4 4 4 4

contains

function myfun(x) result(y)
    real(wp), intent(in) :: x(:)
    real(wp) :: y(size(x))
    y = x**2
end function

end program

For an nD Jacobian you could adapt with something like

function jacobian(f,x) result(J)
    procedure(fun) :: f
    real(wp), intent(in) :: x(:)
    real(wp) :: J(size(x),size(x))
    ! ...
end function
5

I’m curious, are there practical differences in this case, e.g. between stack and heap allocation of the result array?

6

Thank you all for the answers so far!

Maybe you could give func2() another argument that tells it what size y(:) should be, like

@banana-bred that’s what I did after posting, and it works. That’s also the approach that some packages such as MINPACK or GitHub - ivan-pi/stiff3: Adaptive solver for stiff systems of ODEs using semi-implicit Runge-Kutta method of third order use. Eventually if there is no solution I’ll do it, I just wanted to explore if it was possible to do this without explicitly declaring the function return array size.

In the meantime, have you looked at the parameterized derived type (PDT) feature, introduced starting with Fortran 2003, over 20 years ago, but which is still very buggy or dysfunctional in most Fortran processor (compiler) implementations :sob:

@FortranFan I’ll take a look; I’ll confess I don’t really understand yet your code, but you are analytically defining the jacobian, while my objective is to design a function that returns the numerical Jacobian based on an input function with an arbitrary number of variables

May be something like this could help you:

@hkvzjal this is precisely what I need, except for the real(wp) :: y(size(x)) part; ideally, I’d like the function to return the Jacobian of a function like this:

    function myfunc(x) result (y)
        real(wp), intent(in), dimension(2) :: x
        real(wp), dimension(3) :: y

        y(1) = x(1)**2 + x(2)
        y(2) = exp(x(2))
        y(3) = sin(x(2))/x(1)
    
    end function myfunc

Since y doesn’t have the same size as x, right now the options that I have are either passing the size of y explicitly like @banana-bred suggests, or declaring the y as allocatable on the interface:

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

This works, but then every function I try to pass to my jacobian function will need to have allocatable output. Eventually if there’s no possible solution I’ll gravitate towards one of those, but having a flexible interface that works with whatever function I pass it would be ideal.

7

I haven’t followed the thread very well (so maybe missing the point…), but another approach may be to use a subroutine rather than a function, so that either allocatable or non-allocatable arrays can be passed as “result” variable (and defined in-place in that routine)? Indeed, I use subroutines more often than functions (except the result variable is rather simple) because we can control the behavior of variables more explicitly. But functions may also be appealing in that they can be used directly in expressions (so more like math).

Yet another approach might be to use derived types to handle data like x and y (with related info) and calculate their values via user-provided routines (though I’m not sure if this is convenient for the problem…)

8

I saw @septc typing and waited for the answer to be published. In fact, suggesting the same thing I was about to write overlapped. For more flexibility, you might have a look at the subroutined version of your interfaces.

1 Like