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

dfoell · 06-08-2010, 12:11 PM

Hello everyone,

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.

paulsm4 · 06-08-2010, 01:59 PM

Hi -

Q: Did you declare an interface for your function?

dfoell · 06-08-2010, 02:22 PM

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?

dfoell · 06-09-2010, 05:22 AM

Okay, I solved the problem myself

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

)