Fortran How can I properly use function calls in Fortran 90?

[JellyKing]

Hi all,

I'm new to Physics Forums, but not to physics, since I'm now doing my MSc. I've had a persistent problem in a segment of my Fortran code that I've worked around, but would like to fix.

The long and short is that I'm integrating a six-dimensional ODE with a fourth-order Runge-Kutta method.* I've tried to do it with function calls as follows. In the main program body, there will be a line to iterate the solution one step, using some function I call rk. I'd like the line of code to be something like

output = rk(input, parameters)

Sure enough, after the END program, I have the function, initiated by

FUNCTION rk(inputx,parametersx)
REAL(6) :: rk,...
...
rk = ...
END FUNCTION

In the opening lines, I'll declare that rk is a six-dimensional real (REAL(6)). That's fine. However, when I try to call it as above, I get an error like

Error: FUNCTION 'rk' is of rank (1) at {function declaration} and rank (0) and {function call}.

Okay, so I figure I need to tell the program that rk is going to be a REAL(6). But now the compiler tries to parse the arguments as array indices, so I'll get an error like


v = rk(w,E,L,Q,a,dtau)
...1
Error: Array index at (1) must be of INTEGER type


I haven't had this problem with functions that return scalars. For now, I've worked around the function calls by using routines, so that instead of a line like

output = rk(input,parameters)

I use

CALL rk(input,parameters,output)

Sorry this is a long post, and if there are any other details that will help find a solution just ask. Thanks in advance for any help you can provide. This has been bugging me for a while.

JK

*I can't for the life of me work out how one is supposed to import other people's modules, but that's another story.

 

[alphysicist]

Hi JellyKing,

There's not much program in your code to go on, but I guess my first question is do you have an explicit interface to the function? That could be the source of the error.

If not, probably the easiest thing is to show you a working case in case it's helpful; here is a trivial program that uses an array-valued function as an internal procedure (just takes an array and adds 1 to each element):

PROGRAM array_valued_functions
implicit none
integer::i
integer,dimension(5)::a,b

b=(/ (i,i=1,5) /)
a=addone(b)
print*,a

contains

function addone(x)
implicit none
integer,dimension(:),intent(in)::x
integer,dimension(size(x))::addone

addone=x+1
end function addone

end program array_valued_functions

Here is the same example with an interface block which might be better for your case:

PROGRAM array_valued_functions
implicit none
integer::i
integer,dimension(5)::a,b


interface
function addone(x)
integer,dimension(:),intent(in)::x
integer,dimension(size(x))::addone
end function addone
end interface


b=(/ (i,i=1,5) /)
a=addone(b)
print*,a

end program array_valued_functions


function addone(x)
implicit none
integer,dimension(:),intent(in)::x
integer,dimension(size(x))::addone

addone=x+1
end function addone

If you still have trouble, could you perhaps write a toy program that shows the behavior you're talking about and post it here?

 

[JellyKing]

Okay, here's something that reproduces the error.

PROGRAM toy
	REAL, DIMENSION(2) :: z
	
	z(1) = 1
	z(2) = 2
	
	z = array_valued_function(z)
	PRINT *, z
	
END PROGRAM toy

FUNCTION array_valued_function(x)
	REAL, DIMENSION(2) :: array_valued_function,x
	
	array_valued_function = x+1
END FUNCTION

I already notice you making use of a couple of structures I had no idea about. I've never seen interface or contains before. I am fairly new to Fortran, from a smattering of Java, C++, a bit more Python, and a lot of MATLAB. If I include IMPLICIT NONE, I'm told array_valued_function has no implicit type. If I supplant

REAL :: array_valued_function(2)

I get back to my error of the compiler wanting the argument of the function to be an index of the array (described in my first post).

I'll look up those interface and contains structures when I get back to this particular problem. Thanks for the help though! :)

JK

 

[alphysicist]

Then that would be the error; array-valued function require an explicit interface, which can be given in three ways: putting the function in a module and using a corresponding USE statement, making the function an internal subprogram with a CONTAINS statement, or using an interface block in the main (calling) program.

To get your toy program running you might use:

PROGRAM toy
	REAL, DIMENSION(2) :: z

interface
     FUNCTION array_valued_function(x)
     implicit none
     REAL, DIMENSION(2) :: array_valued_function,x
     end function array_valued_function
end interface


	
	z(1) = 1
	z(2) = 2
	
	z = array_valued_function(z)
	PRINT *, z
	
END PROGRAM toy

FUNCTION array_valued_function(x)
	REAL, DIMENSION(2) :: array_valued_function,x
	
	array_valued_function = x+1
END FUNCTION

 

[matTOfort]

Hello and thanks for your info about use of INTERFACE in fortran.

I have a question:

I have a main code within which I have used an external extrinsic function. I have coped and pasted both here:
THIS IS THE MAIN CODE
program matrix
implicit none
double precision, dimension(3,3) :: A
double precision, dimension(3,1) :: B
double precision :: X(3)

integer :: m ! SEE HERE PLEASE

interface
function X_Part_Piv(A,B,n)
double precision, allocatable :: AUG(:,:),AUGP(:,:),AUX1(:,:),AUX2(:,:)
integer :: n,p,k,r,m,o,s
double precision :: big,factor,matrmultip,X_Part_Piv(n),A(n,n),B(n,1)
end function X_Part_Piv
end interface

A(1,1)=9
A(1,2)=33
A(1,3)=-3
A(2,1)=0
A(2,2)=9
A(2,3)=11
A(3,1)=.2
A(3,2)=23
A(3,3)=15
B(1,1)=0
B(2,1)=4
B(3,1)=6

m=size(A,2)! SEE HERE PLEASE

X=X_Part_Piv(A,B,m)

end program matrix

THIS IS THE FUNCTION

function X_Part_Piv(A,B,n)
implicit none

double precision, allocatable :: AUG(:,:),AUGP(:,:),AUX1(:,:),AUX2(:,:)
integer :: n,p,k,r,m,o,s
double precision :: big,factor,matrmultip,X_Part_Piv(n),A(n,n),B(n,1)

n=size(A,2)
p=n+1
allocate (AUG(n,n+1))
do m=1,n
AUG(m,n+1)=B(m,1)
do o=1,n
AUG(m,o)=A(m,o)
end do
end do

allocate (AUGP(n,n+1))

do k=1,n-1

big=maxval(abs(AUG(k:n,k)))
s=maxloc(abs(AUG(k:n,k)),dim=1)!dim=1 enforces the max loc of a an (n,1) matrix to be found
r=s+k-1

if (r/=k) then
AUGP=AUG
AUG(k,:)=AUGP(r,:)
AUG(r,:)=AUGP(k,:)
end if

do m=k+1,n

factor=AUG(m,k)/AUG(k,k)
AUG(m,k:p)=AUG(m,k:p)-factor*AUG(k,k:p)

end do

end do

X_Part_Piv(n)=AUG(n,p)/AUG(n,n)

do m=n-1,1,-1
allocate (AUX1(1,n-m))
allocate (AUX2(n-m,1))
AUX1(1,1:n-m)=AUG(m,m+1:n)
AUX2(1:n-m,1)=X_Part_Piv(m+1:n)
matrmultip=0
do o=1,n-m
matrmultip=matrmultip+AUX1(1,o)*AUX2(o,1)
end do
X_Part_Piv(m)=(AUG(m,p)-matrmultip)/AUG(m,m)
deallocate(AUX1)
deallocate(AUX2)
end do
end function X_Part_Piv


My function has three arguments A,B, and n. n is basically the size of A, so it depends on A.

I cannot omit n from my function`s argument. I have done many try and errors to be able to have X_Part_Piv(A,B) instead of X_Part_Piv(A,B,n)

Can anybody help me with that?