[SOLVED] [Fortran] Array-valued function leads to segfault

ProgrammingThis forum is for all programming questions.
The question does not have to be directly related to Linux and any language is fair game.

Notices

Welcome to LinuxQuestions.org, a friendly and active Linux Community.

You are currently viewing LQ as a guest. By joining our community you will have the ability to post topics, receive our newsletter, use the advanced search, subscribe to threads and access many other special features. Registration is quick, simple and absolutely free. Join our community today!

Note that registered members see fewer ads, and ContentLink is completely disabled once you log in.

I'm trying to implement a solver for a system of differential equations in Fortran. The solver contains a number of functions which are supposed take real values, 1D arrays of real values or both as arguments and return arrays of real numbers, all of which cause the program to segfault.

Example:

Code:

function y_exakt(t)
implicit none
real::t, pi
real,dimension(4)::y_exakt
write(*,*) t
pi = acos(-1.)
y_exakt = (/ cos(10.*t), -10.*sin(10.*t), sin(10.*t), 10.*cos(10.*t) /)
end function

In the main program, the function is declared as real::y_exakt. As soon as the program gets to the last line of the function and the function tries to return the array to the main program, it will crash.

Another function, which takes reals and arrays as arguments, is behaving even more weirdly:

Code:

function eulercauchy(t, y, dt)
implicit none
real::t, dt, f
real,dimension(4)::y, dy, eulercauchy
write(*,*) t !only for
write(*,*) y !debugging
write(*,*) dt !purposes
dy = dt*f(t,y)
write(*,*) dy
eulercauchy = dy
end function

When I pass
t = 0, y = (/ 1, 0, 0, 10 /) and dt = 3.75000015E-02
the function takes
t = 1, y = (/ 3.75000015E-02, 0, 0, 5.60519386E-45 /) and dt = 1
(the last number in the array seems to change randomly). Then the program crashes either when f(t,y) is called or when dy is returned (after removing the call to f).

What could cause these (memory?) problems and/or what could I do to identify the problem? Increasing the maximum stack size with ulimit or compiling the program with -fno-automatic has had no effect.

I'm using gfortran (gcc 4.4.3) on a 64-bit Ubuntu Lucid machine. The complete program can be found at http://pastebin.com/gCUTJmUH. Thanks for your time.

Last edited by dfoell; 06-09-2010 at 05:22 AM.
Reason: problem solved

No, I didn't even know about FORTRAN interfaces until now I'm not exactly a FORTRAN expert, sorry...
Would interface blocks help to identify or solve the problem?
An interface block seemingly contains the headline of the function, the declaration part and the end of the function and enables the compiler to check if the data types match. Is that correct?

Edit: I added an interface block:

Code:

interface
function f(y)
implicit none
real::B
real,dimension(4)::y, f
end function
function y_exakt(t)
implicit none
real::t, pi
real,dimension(4)::y_exakt
end function
function eulercauchy(t, y, dt)
implicit none
real::t, dt, f
real,dimension(4)::y, dy, eulercauchy
end function
function runge_kutta(t, y, dt)
implicit none
real::t, dt, f
real,dimension(4)::y, dy, ya, yb, yc, dya, dyb, dyc, runge_kutta
end function
end interface

... which solved the weird problem with the wrong arguments being passed. But unfortunately, the program will still segfault.
I don't really see what's the right way of declaring a function returning a vector, because e.g. real,dimension(4) function xyz(a,b,c) will not work - every time I try to call the function the compiler obviously thinks I'm trying to access an item in the array ("Error: Rank mismatch in array reference"). Does declaring the type of the function inside the function itself really make any sense?

I removed the "real::f, runge_kutta, eulercauchy, y_exakt" declaration from the main program and the real::f declaration from the runge_kutta and eulercauchy functions.
Then I put all functions behind a CONTAINS inside the main program instead, and voila, no more segfaults.

So I guess the correct way to declare/define a vector function would be e.g.

Code:

contains
function xyz(a)
real::a
real,dimension(3)::xyz
xyz = (/ a, a+1, a+2 /)
end function

(I also defined the function f as f(y) and called it with f(t,y) - DUH )

LinuxQuestions.org is looking for people interested in writing
Editorials, Articles, Reviews, and more. If you'd like to contribute
content, let us know.