Fortran - Calling functions from functions from functions ...

Author

Calling functions from functions from functions ...

Simon Geard

2005-04-22, 8:57 am

Lets say I have a function g(p,x) where x and p are of some defined type.

I also have a function NewtonR which I can use like so:

x = NewtonR(g,p,x0)

There is also another function I want to call ImplicitTrap which can be
done as:

u = ImplicitTrap(g,p,u0,step,n)

This function defines a function of g, and it is this function that I'd
like to pass to NewtonR. Can it be done?

Basically the code is implementing the implicit trapezoidal rule for
numerical integration with the NR solver being used at each step with
automatic differentiation being used to handle the differentiation. My
current version implements NR inside IT but I'd like to use the external
NR if possible.

If you want to see the actual code that's not a problem, I just didn't
want to throw lines of code to the group unnecessarily.

Thanks,

Simon Geard

Arjen Markus

2005-04-22, 8:57 am

Simon Geard wrote:
>
> Lets say I have a function g(p,x) where x and p are of some defined type.
>
> I also have a function NewtonR which I can use like so:
>
> x = NewtonR(g,p,x0)
>
> There is also another function I want to call ImplicitTrap which can be
> done as:
>
> u = ImplicitTrap(g,p,u0,step,n)
>
> This function defines a function of g, and it is this function that I'd
> like to pass to NewtonR. Can it be done?
>
> Basically the code is implementing the implicit trapezoidal rule for
> numerical integration with the NR solver being used at each step with
> automatic differentiation being used to handle the differentiation. My
> current version implements NR inside IT but I'd like to use the external
> NR if possible.
>
> If you want to see the actual code that's not a problem, I just didn't
> want to throw lines of code to the group unnecessarily.
>
> Thanks,
>
> Simon Geard

Hm, you might do it like this:

x = NewtonR(g,ImplicitTrap,p,x0)

with:
real function NewtonR(func, trapf, p, x)

external :: func, trapf

u = trapf(func,p,u0,step,n)

...

(you need to add a few syntactical details of course)

The only thing you can not do (at least not portably) is store the
function argument in some variable for later use. You will need
Fortran 2003 for that.

Regards,

Arjen

Simon Geard

2005-04-22, 3:59 pm

Arjen Markus wrote:
> Simon Geard wrote:
>
>
>
> Hm, you might do it like this:
>
> x = NewtonR(g,ImplicitTrap,p,x0)
>
> with:
> real function NewtonR(func, trapf, p, x)
>
> external :: func, trapf
>
> u = trapf(func,p,u0,step,n)
>
> ...
>
> (you need to add a few syntactical details of course)
>
> The only thing you can not do (at least not portably) is store the
> function argument in some variable for later use. You will need
> Fortran 2003 for that.
>
> Regards,
>
> Arjen

That was my starting point. The problem is then that NewtonR has been
redefined and now has an extra argument which is what I want to avoid.
I've seen implementations of this in python and C++ so I'm now thinking
about the fortran90 and tcl implementations.

Simon

Arjen Markus

2005-04-22, 3:59 pm

Simon Geard wrote:
>

>
> That was my starting point. The problem is then that NewtonR has been
> redefined and now has an extra argument which is what I want to avoid.
> I've seen implementations of this in python and C++ so I'm now thinking
> about the fortran90 and tcl implementations.
>
> Simon

In Tcl this would not be a problem - you can always use [interp alias],
for instance or set a namespace variable to the name of your specific
function.

In Fortran 90 you might achieve the effect via an interface statement
and
modules ... Just thinking out loud ... Is the function ImplicitTrap()
going to
vary per application of NewtonR?

Can you explain the importance of keeping the interface to NewtonR
constant?
(If you add the extra argument, things are much simpler, but I think you
may be able to use some tricks to achieve what you want ...)

Regards,

Arjen

Walter Spector

2005-04-22, 3:59 pm

Simon Geard wrote:
> ... The problem is then that NewtonR has been
> redefined and now has an extra argument which is what I want to avoid.

Perhaps the use of OPTIONAL arguements would help?

Walt
-...-
Walt Spector
(w6ws at earthlink dot net)

NuclearWizard

2005-04-22, 4:00 pm

Simon Geard wrote:
> Lets say I have a function g(p,x) where x and p are of some defined
type.
>
> I also have a function NewtonR which I can use like so:
>
> x = NewtonR(g,p,x0)
>
> There is also another function I want to call ImplicitTrap which can
be
> done as:
>
> u = ImplicitTrap(g,p,u0,step,n)
>
> This function defines a function of g, and it is this function that
I'd
> like to pass to NewtonR. Can it be done?
>
> Basically the code is implementing the implicit trapezoidal rule for
> numerical integration with the NR solver being used at each step with

> automatic differentiation being used to handle the differentiation.
My
> current version implements NR inside IT but I'd like to use the
external
> NR if possible.
>
> If you want to see the actual code that's not a problem, I just
didn't
> want to throw lines of code to the group unnecessarily.
>
>
> Thanks,
>
> Simon Geard

Yo Simon,
I ran into this problem as well. This is the "fix".

1. This is what NewtonR should look like:
FUNCTION NewtonR( g , p , x0 , u ) RESULT(x)
!g,p,x0 kind-type declarations
OPTIONAL :: u
INTERFACE
FUNCTION u(g,p,u0,step,n)
!g,p,u0,step,n,u kind-type declarations
ENDFUNCTION
ENDINTERFACE
I am just typing, I have not compiled the previous statements. I could
be wrong but I am pretty sure all the above statements will compile on
Compaq Visual Fortran 6.6.

2. All your integration functions, like ImplicitTrap, should have the
same arguments and kind-type declarations for those arguments to
facilitate passing to NewtonR. There is no runtime-polymorphism in
Fortran! Good for speed but bad because it requires some extra
compiler logic to catch the passing of a function to a function.
Basically the compiler must be required to allow interface statements
to contain kinds and types defined in the scope of the interface
statement. As an example, consider the following situation which I
think should compile but DID NOT.

FUNCTION NewtonR( g , p , x0 , u ) RESULT(x)
!#PURPOSE#
! This function solves ODEs.

!#METHOD#
! Newton-Raphson.
! <insert some math here>

!#PARAMETERS#
! @ real kind for the left side matrix [KND_g]
! @ real kind for the right side matrix [KND_p]
! @ real kind for the solution [KND_x]
INTEGER,PARAMETER :: KND_g = 8
INTEGER,PARAMETER :: KND_p = 8
INTEGER,PARAMETER :: KND_x = 4

!#REQUIRED INPUT#
! @ left hand side matrix [g]
! @ right hand side matrix [p]
! @ initial value [x0]
REAL(KND_g),INTENT(IN) :: g (:,;)
REAL(KND_p),INTENT(IN) :: p (1:SIZE(g,1),1:SIZE(g,2))
REAL(KND_x),INTENT(IN) :: x0( 1:SIZE(g,2))

!#OPTIONAL INPUT#
! @ integration function [u]
OPTIONAL :: u
INTERFACE
FUNCTION u(g,p,u0,step,n)
REAL(KND_g) ,INTENT(IN) :: g(:,:)
REAL(KND_p) ,INTENT(IN) :: p (1:SIZE(g,1),1:SIZE(g,2))
REAL(KND_x) ,INTENT(IN) :: u0( 1:SIZE(g,2))
INTEGER,OPTIONAL,INTENT(IN) :: step,n
ENDFUNCTION
ENDINTERFACE

!#REQURIED OUTPUT#
! @ solution vector [x]
REAL(KND_x) :: x (1:SIZE(g,2))

!#BEGIN#

<some setup code>

<some solver loop with the following lines>
IF( PRESENT(u) )THEN
<execute integration function passed into this routine>
ELSE
<execute default integration function>
ENDIF

<post processing>

!#END#
ENDFUNCTION

The compiler flags on the inclusion of KND_x, KND_g, etc. in the
interface. (If someone uses a compiler that actually understands
this---could you post it?) So the way you would have to do it is
having a module that contains a shit ton of kinds and types like KND_x
and KND_g and countless other kinds that all your other functions
need---I hate this way---it seems very messy and KND_x to one function
may be the dependent variable and KND_x to another function may be
independent! This would lead to the tendency to create your own kinds
to add to the "master kind module" and then a user of your codes must
go into this kind file and set and choose the kinds they need for
specific tasks. The better approach is to just have overloaded
functions able to accept multiple kinds and types---all known AT
COMPILE TIME---which then an interface selects the specific procedure
it needs from the generic procedure the user called. In this way,
users have a straightforward way of dictating which procedures are
used---it is the procedure that does what the generic name says
(ImplicitTrap) with the specific algorithm designed to best handle the
arguments they provided. (I would name the ImplicitTrap procedure
something like IntegrateITrap.) Parameters are some values that are
SET at compile time. The compiler should just replace all instances of
<kind-variable> with <kind-value> and then compile the resulting code.
With this type of procedure overloading, you need a procedure for every
type-kind combination and then an interface inside the preamble of a
module:
INTERFACE GenericName
MODULE PROCEDURE SpecificTypeKindCombo1
MODULE PROCEDURE SpecificTypeKindCombo2
MODULE PROCEDURE SpecificTypeKindCombo3
MODULE PROCEDURE SpecificTypeKindCombo4
...
ENDINTERFACE

CONTAINS
and then in the body:
<all the procedures for the specific types/kinds>

THERE are no templates in Fortran90/95! But Copy/Paste/Replace all
instances of KIND1 with KIND2 works okay---it just looks ugly to have
all that repetition. When I code an algorithm, I like to keep it in a
seperate file and then INCLUDE it in functions where various types are
declared, as long as the operators in the algorithm are defined for the
actual types, you will get multiple uses out of this single piece of
algorithm code.

Hope this helps,
-Will

Simon Geard

2005-04-25, 8:58 am

>>Simon Geard

<snip stuff="original post"/>

>
>
> Yo Simon,
> I ran into this problem as well. This is the "fix".
>

<snip stuff="lots of good example code"/>

Will,

Thanks for this. It'll be a few days before I can try it out - I'll let
you know how I get on.

Regards,

Simon

Sponsored Links