internal procedure as an actual argument
qsc
Suppose we need to integrate a function f(x,t) about the x variable,
for some time t. So we can view t as a parameter. We are going to, and
we have to, make use of an integrator, say integration_routine, which
takes a function in one variable as an argument. Hence we need to
somehow hide the parameter t from the integrator, but we still want to
be able to supply the information when needed, for we want to be able
to integrate the function for different times.
I thought of internal procedures. I thought that writing the funtion f
as an internal procedure should solve the problem. And I even thought
"finally I come to a situation when an internal function is absolutely
necessary". Sadly, as you may know, this did not work out. My g95
compiler told me that an internal procedure can't be an actual
argument.
I will post a copy of the program that I described above, so that you
will know what I am trying to do. And you may be able to help me
out ...
!!!!!!!!!!!Fortran code !!!!!!
subroutine compute_integral( )
real :: result
real :: t
t = 0.5
call integration_routine( f, result )
....
t = 1.5
call integration_routine( f, result )
....
contains
function f( x )
real :: x
real :: f
f = exp(x*t)
end function
end subroutine compute_integral
!!!!!!!!!!!!!!! Fortran code ends !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Code like this won't work. What would you recommend? The only way I
can think of
is to make the function f as an external procedure, and make t a
global variable, then
put t, compute_integral, and f in one module.I am not sure if this
approach works
or not. I would rather not to take this approach if there is any
other better and safer
ways to accomplish the same task.
e p chandler
Put T in a module by itself. Then USE that module where needed:
--- example program ---
module arg
implicit none
real t
end module
implicit none
real,external :: f
integer n
do n=1,5
call compute_integral(f,1.0,2.0,n)
end do
end
subroutine compute_integral(f,x1,x2,n)
use arg
implicit none
real f,x1,x2,u,v
integer n
t=0.5
call integration_routine(f,x1,x2,n,u)
t=1.0
call integration_routine(f,x1,x2,n,v)
print *,n,u,v
end
real function f(x)
use arg
implicit none
real x
f=exp(x*t)
end
subroutine integration_routine(f,x1,x2,n,v)
implicit none
real f,x1,x2,v
integer n
real x,dx,hx,gx,s
integer i
dx=(x2-x1)/n
hx=dx/2.0
gx=hx/sqrt(3.0)
s=0
do i=0,n-1
x=x1+(i+0.5)*dx
s=s+f(x-gx)+f(x+gx)
end do
s=s*hx
v=s
end
--- end example ---
Note: This sample program uses 2 point Gaussian quadrature. (2 points
per cell) This technique can produce remarkably good results for some
integrands even with low N and a low number of points per cell.
-- e
qsc
Hi, thanks for the reply, and thanks for the elegant sample code. Your
approach seems better thank what I can thought of. I think I will try
it out, to see if it can make my program better. Thank you!
glen herrmannsfeldt
> Suppose we need to integrate a function f(x,t) about the x variable,
> for some time t. So we can view t as a parameter. We are going to, and
> we have to, make use of an integrator, say integration_routine, which
> takes a function in one variable as an argument. Hence we need to
> somehow hide the parameter t from the integrator, but we still want to
> be able to supply the information when needed, for we want to be able
> to integrate the function for different times.
The traditional method is through COMMON.
> I thought of internal procedures. I thought that writing the funtion f
> as an internal procedure should solve the problem. And I even thought
> "finally I come to a situation when an internal function is absolutely
> necessary". Sadly, as you may know, this did not work out. My g95
> compiler told me that an internal procedure can't be an actual
> argument.
The problems with internal procedures as actual arguments has
been discussed before. It gets much more complicated when you
allow recursion. PL/I, which does allow it, passes enough
information to find the block active at the time of the call
(or assignment to an ENTRY variable).
http://publibz.boulder.ibm.com/cgi-bin/bookmgr_OS390/BOOKS/ibm3lr30/6.10
There are complications in adding this to Fortran related
to back compatibility.
-- glen
Steve Lionel
> I thought of internal procedures. I thought that writing the funtion f
> as an internal procedure should solve the problem. And I even thought
> "finally I come to a situation when an internal function is absolutely
> necessary". Sadly, as you may know, this did not work out. My g95
> compiler told me that an internal procedure can't be an actual
> argument.
Intel Fortran supports this as an extension - the standard even has
wording suggesting how it ought to work if a compiler implements it.
There may be other compilers that support it as well.
Steve
relaxmike
on function pointers (for example the cray pointers of gfortran)
included in a OO integration class.
The following object-based module is a sample integration class, which
allows
to store various settings such as the solution, plus a function
pointer
to store the adress of the function to integrate.
(I am not sure that the source compiles, but this is the idea)
module integration
implicit none
public :: integration_new
public :: integration_free
public :: integration_set_integration_function
public :: integration_solve
public :: integration_get_solution
!
! Derived-type to make integration !
!
type, public :: T_INTEGRATION
private
! Address of the function to integrate
integer :: fonc_address
! Initial guess, solution, etc... comes here
double precision :: solution
end type T_INTEGRATION
contains
!
! Set the basic parameters of the optimizer.
!
subroutine integration_new ( this )
implicit none
type(T_INTEGRATION), intent(inout) :: this
!
! Initialize the fonc pointer
!
this % fonc_address = 0
end subroutine integration_new
!
! Frees memory for pointers in derived types
!
subroutine integration_free ( this )
implicit none
type(T_INTEGRATION), intent(inout) :: this
end subroutine integration_free
!
! Set the objective function to optimize
!
subroutine integration_set_integration_function ( this , new_fonc )
implicit none
type ( T_INTEGRATION ), intent(inout) :: this
interface new_fonc_interface
double precision function new_fonc ( x )
implicit none
double precision, intent(in) :: x
end function new_fonc
end interface new_fonc_interface
this % fonc_address = loc ( new_fonc )
end subroutine integration_set_integration_function
!
! Integration solver
!
subroutine integration_solve ( this )
implicit none
type ( T_INTEGRATION ), intent(inout) :: this
! Solver ...
! Here, the function can be computed with the subroutine
integration_fonc
end subroutine integration_solve
!
! Accessor to the solution
!
subroutine integration_get_solution ( this , solution )
implicit none
type ( T_INTEGRATION ), intent(inout) :: this
double precision, intent(out) :: solution
solution = this % solution
end subroutine integration_get_solution
!
! This fonction allows the integrator to compute the function
!
double precision function integration_fonc ( this , x )
implicit none
type ( T_INTEGRATION ), intent(inout) :: this
double precision, intent(in) :: x
double precision :: fonc_pointee
external fonc_pointee
pointer ( fonc_pointer , fonc_pointee )
!
! Evaluate the function
!
fonc_pointer = this % fonc_address
optim_fonc = fonc_pointee ( x )
end function integration_fonc
end module integration
The client program has to set the function by the way of the
integration_set_integration_function accessor.
The following example shows two functions :
* f1 only uses a constant parameter
* f2 get that parameter from a module.
program client_integration
implicit none
use integration, only : integration_new, &
integration_free, &
integration_set_integration_function, &
integration_solve, &
integration_get_solution
type ( T_INTEGRATION ) :: myintegration
double precision :: solution
interface
double precision function f1 ( x )
implicit none
double precision, intent(in) :: x
end function f1
double precision function f2 ( x )
implicit none
double precision, intent(in) :: x
end function f2
end interface
call integration_new ( myintegration )
call integration_set_integration_function ( myintegration , f1 )
call integration_solve ( myintegration )
call integration_get_solution ( myintegration , solution )
call integration_free ( myintegration )
end program client_integration
function f1( x )
real :: x
real, parameter :: f = 2.d0
f = exp(x*t)
end function f1
function f2( x )
use myparam, only : f
real :: x
f = exp(x*t)
end function2
module myparam
implicit none
private
real, parameter, public :: f = 2.d0
end module myparam
Indeed, these problems really are good examples for that way
of programming. It includes optimization solvers, Newton
methods, ode solvers, statistics analysis etc...
Best regards,
Michaël
qsc
It is very interesting! However I couldn't get your code compiled with
gfortran ... There are many language mistakes that I am not able to
fix ...
After roughly going through your code, I would like to ask
1. In function f1, the first f really should be t, and the second f
should be f1? And in function f2, f should be f2? In module myparam, f
should be t?
2. Have you merely given a skeleton of the program for integration, or
your program really can do integration after the bugs are killed?
Mirko....@gmail.com
This was one USEFUL thread (among many). Thank you for the examples.
Mirko
relaxmike
But it should be clear that I did not create an integration method
from scratch in 2 minutes... It is just the structure into which an
integration method can be put. The question here is the fortran,
not the mathematics, isn't it ?
The goal is to consider that example to apply OO methods
to numerical algorithms in fortran.
The other point is to show an application of the cray pointers
which allow to uncouple the numerical algorithm from the particular
function to evaluate.
All in all, this is the routine for C/C++ developpers (and for many
years),
but it is not so common in fortran.
These are the fixed files.
integration.f90
integration_fonc = fonc_pointee ( x )
end function integration_fonc
end module integration
myparam.f90 :
module myparam
implicit none
private
real, parameter, public :: t = 2.d0
end module myparam
client_integration.f90 :
program client_integration
use integration, only : integration_new, &
integration_free, &
integration_set_integration_function, &
integration_solve, &
integration_get_solution, &
T_INTEGRATION
implicit none
type ( T_INTEGRATION ) :: myintegration
double precision :: solution
interface
double precision function f1 ( x )
implicit none
double precision, intent(in) :: x
end function f1
double precision function f2 ( x )
implicit none
double precision, intent(in) :: x
end function f2
end interface
call integration_new ( myintegration )
call integration_set_integration_function ( myintegration , f1 )
call integration_solve ( myintegration )
call integration_get_solution ( myintegration , solution )
call integration_free ( myintegration )
end program client_integration
double precision function f1( x )
implicit none
real :: x
real, parameter :: t = 2.d0
f1 = exp(x*t)
end function f1
double precision function f2 ( x )
use myparam, only : t
implicit none
real :: x
f2 = exp(x*t)
end function f2
These are the compile lines :
gfortran -fcray-pointer -c integration.f90
gfortran -fcray-pointer -c myparam.f90
gfortran -fcray-pointer -c client_integration.f90
gfortran -fcray-pointer client_integration.o integration.o myparam.o -
o integr.exe
When you execute it, well nothing happens obviously....
Best regards,
Michaël
Arjen Markus
Fortran 2003 does allow this sort of things in a standard
way. See for instance the book Fortran 95/2003 explained
by Metcalf, Reid and Cohen.
The problem with previous versions of Fortran is
that the type checking prevents you from passing along arbitrary
data. (In C, this is easy: just use void *).
An alternative is to use the transfer() function, but that
tends to scare people off ;).
(It is actually one subject in an article I wrote for the
next issue of Fortran Forum - Fortran 2003 and design patterns.)
Regards,
Arjen
Steve Lionel
> OK, this is a fixed version which compiles.
I will point out, though, that this code uses an extension of a Cray
pointee being a procedure. Many compilers that support Cray pointers
(themselves an extension) don't support procedure pointees.
Steve
relaxmike
but the good news is that I can learn !
Arjen, can you show an example of how the use of the "transfer" method
allows to solve the problem ?
If I understand that page :
http://macresearch.org/advanced_fortran_90_callbacks_with_the_transfer_function
the transfer function allows to create generic callbacks by passing
the
function as an argument, isn'it ?
What is not standard fortran 90 is the previous is
the use of the cray pointers.
Instead, the fortran 2003 includes the "procedure pointer"
which is similar in the use, although it has nothing
in common to cray pointers in the form.
(I see that here : p. 119 of http://www.idris.fr/data/cours/lang/fortran/f2003/Fortran_2003.pdf)
Considering the fact that the "this" derived type argument
being the first argument of each subroutine, I suppose that we all
agree that this is standard and will be available as the "type bound
procedure" method in the fortran 2003 standard :
(see p. 151 of http://www.idris.fr/data/cours/lang/fortran/f2003/Fortran_2003.pdf)
All in all, what matters here is the fact that pointers allow to
uncouple
the numerical method associated with the integration method from the
specific function which is to be integrated.
Another solution would be to take the function to integrate as an
argument
of the integrator.
The advantage of the cray pointers over the previous solution is that
the
function can then be evaluated directly, at any level of the call
tree, without
explicitely passing the argument-function at all levels of the call
tree.
By the way, it is not clear to me that "procedure pointers" allow that
(may
be we should define that in an example and test it against a f03
compiler).
Anyway, all this is very nice and clean (and fun...), comparing to the
old-school fortran 77 !
Best regards,
Michaël