CODEWITHSUNDEEP
fortranparameter-passinggfortran
Volker
Volker

Reputation: 210

Passing an internal procedure as argument

I want to solve a differential equation lots of times for different parameters. It is more complicated than this, but for the sake of clarity let's say the ODE is y'(x) = (y+a)*x with y(0) = 0 and I want y(1). I picked the dverk algorithm from netlib for solving the ODE and it expects the function on the right hand side to be of a certain form. Now what I did with the Intel Fortran compiler is the following (simplified):

subroutine f(x,a,ans)
implicite none
double precision f,a,ans,y,tol,c(24),w(9)
...
call dverk(1,faux,x,y,1.d0,tol,ind,c,1,w)
...
contains
    subroutine faux(n,xx,yy,yprime)
           implicite none
           integer n
           double precision xx,yy(n),yprime(n)
           yprime(1) = (yy(1)+a)*xx
    end subroutine faux
end subroutine f

This works just fine with ifort, the sub-subroutine faux sees the parameter a and everything works as expected. But I'd like the code to be compatible with gfortran, and with this compiler I get the following error message:

Error: Internal procedure 'faux' is not allowed as an actual argument at (1)

I need to have the faux routine inside f, or else I don't know how to tell it the value of a, because I can't change the list of parameters, since this is what the dverk routine expects.

I would like to keep the dverk routine and understand how to solve this specific problem without a workaround, since I feel it will become important again when I need to integrate a parameterized function with different integrators.

Upvotes: 6

Views: 7361

Answers (3)

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

Reputation: 60008

Passing an internal procedure is allowed by Fortran 2008. That means that your code is correct Fortran and just a more recent compiler is needed.

At the time when this question was posted the support for that in existing compilers was patchy. Some did support it, some did not. This changed considerably and at the time of writing this revision of the answer one can normally use it without needing to think about compiler support any more.

You are not limited to sharing the data in a module. You just have to make sure, that the containing procedure is still running. It is not possible to save the pointer and make a so called 'closure' known from LISP and many modern languages. In languages like Fortran or C the enclosing environment for future calls is not saved and pointers to internal functions become invalid.

The advantage of internal procedure over module procedure is that you can share the data without sharing it in the whole module.

It works well will reasonably old versions of GCC (gfortran), Intel Fortran and Cray Fortran (personally tested) and possibly others. I remember not being able to compile by PGI and Oracle Fortran. As shown by Daniel, it can be checked at http://fortranwiki.org/fortran/show/Fortran+2008+status (search Internal procedure as an actual argument). The most recent version of this table is periodically published in ACM SIGPLAN Fortran Forum.

There is some interesting info in this article https://stevelionel.com/drfortran/2009/09/02/doctor-fortran-in-think-thank-thunk/

Be aware that the implementation of this feature normally requires executable stack due to the use of trampolines. One does not normally think about that but if the OS forbids executable stack, the executable might fail (there were issues in initial versions of Windows Subsystem for Linux, for example).

Upvotes: 3

bdforbes
bdforbes

Reputation: 1546

You could put this all in a module, and make a a global module variable. Make faux a module procedure. That way, it has access to a. 

module ode_module
    double precision::a

    contains

    subroutine f(x,a,ans)
        implicit none
        double precision f,ans,y,tol,c(24),w(9)

        call dverk(1,faux,x,y,1.d0,tol,ind,c,1,w)

    end subroutine 

    subroutine faux(n,xx,yy,yprime)
       implicite none
       integer n
       double precision xx,yy(n),yprime(n)
       yprime(1) = (yy(1)+a)*xx
    end subroutine faux

end module

Upvotes: 5

Daniel_H
Daniel_H

Reputation: 105

You can check here which compilers support this F2008 feature:

http://fortranwiki.org/fortran/show/Fortran+2008+status

On that page, search for "Internal procedure as an actual argument" .

Upvotes: 4

Related Questions