Fortran - Passing array valued functions as argument to function.

Author

Passing array valued functions as argument to function.

andy2o

2004-11-19, 3:58 pm

Hi all,

I've spent several years working with Fortran 95, yet somehow I've
never needed to write functions which use other functions as arguments
until today. I've struggled to make sense of this from the books I
have to hand (no M&R just now).

The purpose of the code will be to compute a definite integral of the
passed function. Here's a toy code example of what I naively expected
to work:

MODULE one
CONTAINS
FUNCTION f1(x_in) RESULT (y_out)
DOUBLE PRECISION :: x_in(:)
DOUBLE PRECISION :: y_out(SIZE(x_in))
y_out = 2*x_in
END FUNCTION f1
END MODULE one

PROGRAM test
USE one
DOUBLE PRECISION :: a(3)=(/ 1.0d0,2.0d0,3.0d0 /)
PRINT*, f1(a)
PRINT*, eval(a,f1)

CONTAINS
FUNCTION eval(x_in,func) RESULT (y_out)
DOUBLE PRECISION :: x_in(:)
DOUBLE PRECISION :: y_out(SIZE(x_in))
DOUBLE PRECISION, EXTERNAL :: func

y_out = func(x_in)
END FUNCTION eval
END PROGRAM test

but I get:

bash$ f95 test.f90 #(Note f95=Nag's compiler on linux.)
Error: test.f90, line 14: Array supplied for scalar argument FUNC (no.
2) of EVAL
[f95 error termination]

Could someone explain how this is meant to be done in Fortran 95.
Many thanks,
Yours,
andy.

PS:
I'd also much appreciate an explanation of why when I replace the
function f1 with

FUNCTION f1(x_in) RESULT (y_out)
DOUBLE PRECISION :: x_in(:)
DOUBLE PRECISION :: y_out
y_out = 2*x_in(1)
END FUNCTION f1

the code segfaults instead.

Jan Vorbrüggen

2004-11-19, 3:58 pm

CVF 6.6C says:
Error: The shape matching rules of actual arguments and
dummy arguments have been violated. [F1]
PRINT*, eval(a,f1)
-----------------^

which is because you declare the second argument of eval to be

> DOUBLE PRECISION, EXTERNAL :: func

while f1 returns a rank-one double-precision array. If, instead, you
replace this line by

INTERFACE
FUNCTION func(x_in) RESULT (y_out)
DOUBLE PRECISION :: x_in(:)
DOUBLE PRECISION :: y_out(SIZE(x_in))
END FUNCTION func
END INTERFACE

at least CVF compiles the code and produces the correct result:

2.00000000000000 4.00000000000000 6.00000000000000
2.00000000000000 4.00000000000000 6.00000000000000

Jan

Jugoslav Dujic

2004-11-19, 3:58 pm

andy2o wrote:
| Hi all,
|
| I've spent several years working with Fortran 95, yet somehow I've
| never needed to write functions which use other functions as arguments
| until today. I've struggled to make sense of this from the books I
| have to hand (no M&R just now).
|
| The purpose of the code will be to compute a definite integral of the
| passed function. Here's a toy code example of what I naively expected
| to work:
|
| MODULE one
| CONTAINS
| FUNCTION f1(x_in) RESULT (y_out)
| DOUBLE PRECISION :: x_in(:)
| DOUBLE PRECISION :: y_out(SIZE(x_in))
| y_out = 2*x_in
| END FUNCTION f1
| END MODULE one
|
| PROGRAM test
| USE one
| DOUBLE PRECISION :: a(3)=(/ 1.0d0,2.0d0,3.0d0 /)
| PRINT*, f1(a)
| PRINT*, eval(a,f1)
|
| CONTAINS
| FUNCTION eval(x_in,func) RESULT (y_out)
| DOUBLE PRECISION :: x_in(:)
| DOUBLE PRECISION :: y_out(SIZE(x_in))
| DOUBLE PRECISION, EXTERNAL :: func
|
| y_out = func(x_in)
| END FUNCTION eval
| END PROGRAM test
|
| but I get:
|
| bash$ f95 test.f90 #(Note f95=Nag's compiler on linux.)
| Error: test.f90, line 14: Array supplied for scalar argument FUNC (no.
| 2) of EVAL
| [f95 error termination]
|
| Could someone explain how this is meant to be done in Fortran 95.

That was a nice (almost good) try. However, as you probably know:

- functions returning arrays require an explicit interface
- funcions with assumed-shape dummies require an explicit interface.

Your f1 above fulfills both conditions.

You don't have an explicit interface there -- only DOUBLE PRECISION,
EXTERNAL, which is not good enough. Even if you did USE one, and
"Eval" is internal to the program, "Eval" still does not know that
its dummy argument "Func" is associated with "F1". Compilers are
not required to look that deep into your thoughts.

Since func is a dummy argument, the only means to achieve the explicit
interface is the INTERFACE block. Instead of

DOUBLE PRECISION, EXTERNAL :: func

use
INTERFACE
FUNCTION func(x_in) RESULT(y_out)
DOUBLE PRECISION :: x_in(:)
DOUBLE PRECISION :: y_out(SIZE(x_in))
END FUNCTION
END INTERFACE

and there you are.

The second sigsegv is again consequence of lack of explicit interface --
without it, "Eval" "presumes" it has assumed-size rather than
assumed-shape dummy argument, and that mismatch causes the sigsegv.

--
Jugoslav
___________
www.geocities.com/jdujic

Please reply to the newsgroup.
You can find my real e-mail on my home page above.

andy2o

2004-11-22, 9:02 am

"Jugoslav Dujic" <jdujic@yahoo.com> wrote in message news:<306egrF2t1046U1@uni-berlin.de>...
[snip]
> You don't have an explicit interface there -- only DOUBLE PRECISION,
> EXTERNAL, which is not good enough. Even if you did USE one, and
> "Eval" is internal to the program, "Eval" still does not know that
> its dummy argument "Func" is associated with "F1". Compilers are
> not required to look that deep into your thoughts.

Dear Jan and Jugoslav,

Thanks for your replies, they solved my problem. There's a lot more
Fortran 95 left for me to learn, and I will read up on INTERFACE
blocks!

Thanks,
andy

Jan Vorbrüggen

2004-11-22, 9:02 am

> Thanks for your replies, they solved my problem. There's a lot more
> Fortran 95 left for me to learn, and I will read up on INTERFACE blocks!

I believe the example we gave you is the 1 in the 1.5 reasons to actually
use INTERFACE blocks. The other 0.5 is when interfacing with library code
you do not control, and for which no module file are available. All other
documented uses should not be used, because better alternatives are available.
So you might as well read up other things, but not the gory details of
INTERFACE blocks.

Jan

Richard E Maine

2004-11-22, 3:58 pm

Jan Vorbrüggen <jvorbrueggen-not@mediasec.de> writes:

> I believe the example we gave you is the 1 in the 1.5 reasons to actually
> use INTERFACE blocks....

> So you might as well read up other things, but not the gory details of
> INTERFACE blocks.

Yes. Some textbooks have been known to introduce interface bodies
too early, with the result that people try to overuse them (in
my opinion). I really recommend that you *NOT* use them much;
you just happen to have found one of the few places where you
do need them. Some people end up trying to use them in places
where they not only aren't needed, but aren't even allowed.

I'd agree with Jan's 1.5 reasons, but to me this is the 1/2. Almost
all of my interface bodies are for external library routines,
particularly C ones. So I'd label that as the 1 and dummy procedures
as the 1/2, but we are close enough.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Jan Vorbrüggen

2004-11-23, 8:57 am

> I'd agree with Jan's 1.5 reasons, but to me this is the 1/2. Almost
> all of my interface bodies are for external library routines,
> particularly C ones. So I'd label that as the 1 and dummy procedures
> as the 1/2, but we are close enough.

Hmmm - I put the dummy procedures as counting fully because in this case
there is a requirement to use them in scenarios as posited by the OP.
The case of the external libraries count, IMO, as the 0.5 because they
are a nice-to-have, additional-safety-net, and-so-on, feature, but not a
requirement to get compileable code.

Jan

Sponsored Links