CODEWITHSUNDEEP
functionfortran
user1543042
user1543042

Reputation: 3440

Fortran passing statement function to function

I am looking to pass a statement function to be evaluated in another function, is this possible? I am trying to have a function which can numerically integrate a general function. There are other values in my real problem so I don't want to have an external function.

Edit 2014/11/1 02:22 GMT

There were a couple of bugs in the code @IanH put in. This is the code which works after I take care of the bugs.

program trapezoidZ
  implicit none

  ! main program global variables and constants here.

  call main
contains
subroutine main
    real :: u, integral1
    integer :: i
    integer, parameter :: n = 5000
    ! delete declaration of trapezoid_integration and f.

    integral1 = trapezoid_integration(TrigSin, n, 0., 3.14159 / 2.0)

    write (*,*) "I - Trap = ",1 - integral1
end subroutine main

function my_f(u)
    real :: my_f
    real, intent(in) :: u

    my_f = (u**4)*EXP(u)/((EXP(u)-1.0)**2)
end function my_f

function TrigSin(u) result(my_f)
    real :: my_f
    real, intent(in) :: u

    my_f = sin(u)
end function TrigSin

function trapezoid_integration (f, n, start_val, end_val) result (integral)
    ! delete declaration of integrand (assuming you have it).
    procedure(my_f) :: f

    integer :: n
    real :: start_val, end_val
    real :: integral,u,h
    integer :: i

    integral = 0.0

    do i=0,n
        u = start_val + ((end_val - start_val)*i)/n

        ! implement Eqn G.4
        if ((i.eq.0).or.(i.eq.n)) then
            integral = integral+integrand(f, u)
        else
            integral = integral+(2.0*integrand(f, u))
        end if

    end do

    h=(end_val - start_val)/n
    integral = (h/2.0)*integral

end function trapezoid_integration

function integrand(f, x) result (value)
    implicit none
    real :: x
    real :: value
    real :: f

    if (x .lt. 0.00001) then
        x = 0.00001
    end if

    value = f(x)
end function integrand

end program trapezoidZ

Below this is the original code. DO NOT USE

program trapezoidZ
  implicit none
  integer, parameter :: n = 10
  real :: u, integral
  integer :: i
  real :: f, trapezoid_integration

  f(u) = (u**4)*EXP(u)/((EXP(u)-1.0)**2)

  integral = trapezoid_integration(f, n, 0.2, 3.0)
  write (*,*) '#trapezoidal integration = ',integral

end program trapezoidZ

function trapezoid_integration(f, n, start_val, end_val)
  implicit none
  integer :: n
  real :: start_val, end_val
  real :: f
  real :: integral,u,h
  integer :: i

  integral = 0.0

  do i=0,n
     u = start_val + ((end_val - start_val)*i)/n

     ! implement Eqn G.4
     if ((i.eq.0).or.(i.eq.n)) then
        integral = integral+integrand(f, u)
     else
        integral = integral+(2.0*integrand(f, u))
     end if

  end do

  h=(end_val - start_val)/n
  integral = (h/2.0)*integral

end subroutine trapezoid_integration

function integrand(f, x) result (value)
  implicit none
  real :: x
  real :: value
  real :: f

  if (x .lt. 0.00001) then
     x = 0.00001
  end if

  value = f(x)
end function integrand

Upvotes: 0

Views: 377

Answers (2)

High Performance Mark
High Performance Mark

Reputation: 78364

Paragraph 12.6.4.4 of the draft of the 2008 standard that I have to hand states:

A statement function shall not be supplied as an actual argument.

so I think you cannot do what you want to do in the way you want to do it. I wonder what your compiler told you when you tried to compile the code above. When I corrected the mismatch between function trapezoid_integration and end subroutine trapezoid_integration in your post my compiler (a recent edition of gfortran) complained, pointing at the use of f in the call to trapezoid_integration:

Error: Statement function 'f' at (1) is not allowed as an actual argument

It complained about several other things too. It might have pointed out, but didn't, that statement functions have been obsolescent since the 1990 standard.

There is a variety of ways to fix your program, one of which is to make the function to be integrated internal, that is, declare it inside a contains section within the program. In fact, I'd contain all the functions you've shown inside such a section within the program. As you've written the code the two procedures trapezoid_integration and integrand are both external and the compiler can't check their interfaces when they are called. Rewrite f as a function and insert contains where you have end program and move end program to the end of the source file, then tidy up the remaining errors.

Upvotes: 1

IanH
IanH

Reputation: 21451

It is not permitted in standard Fortran - a statement function is not in the list of permitted procedure types in constraint C1235 in Fortran 2008.

However, as that constraint indicates, you are not restricted to external functions. Since Fortran 90 the actual argument can be a module procedure. As of Fortran 2008 it can be an internal procedure.

That offers the possibility of compact and much more robust source:

program trapezoidZ
  implicit none

  ! main program global variables and constants here.

  call main
contains
  subroutine main
    real :: u, integral
    integer :: i
    integer, parameter :: n = 10
    ! delete declaration of trapezoid_integration and f.

    integral = trapezoid_integration(my_f, 0.2, 3.0)
    write (*,*) '#trapezoidal integration = ',integral
  end subroutine main

  function my_f(u)
    real :: my_f
    real, intent(in) :: u

    my_f = (u**4)*EXP(u)/((EXP(u)-1.0)**2)
  end function my_f

  function trapezoid_integration(f, n, start_val, end_val)
    ! delete declaration of integrand (assuming you have it).
    procedure(my_f) :: f
    ...
  end function trapezoid_integration

  function integrand(f, x) result (value)
    ...
  end function integrand
end program trapezoidZ

Upvotes: 4

Related Questions