Function Returning an array in Fortran

Question

It is my understanding that you can return an array from a function in Fortran, but for some reason my code is only returning the first value in the array I am asking it to return. This is the function:

function polynomialMult(npts,x,y)
    integer npts
    double precision x(npts), results(npts + 1), y(npts,npts)

    polynomialMult =  x(1:npts) + 1

end function

and this is where I'm calling it

 C(1:numPoints) = polynomialMult(numPoints,x,f)

print *, C(1:numPoints)`

right now it doesn't do anything useful because I am trying to understand the syntax before I write the logic. I saw some stuff about specifying types for functions, but when I write

integer function polynomialMult(npts,x,y)

or whatever I get a compilation error.

Solution 1

To define a function which returns an array include the function declaration inside the function, like this:

function polynomialMult(npts,x,y)
    integer npts
    double precision x(npts), results(npts + 1), y(npts,npts)

! Change the next line to whatever you want
    double precision, dimension(npts) :: polynomialMult

    polynomialMult =  x(1:npts) + 1

end function

Your declaration

integer function polynomialMult(npts,x,y)

declares that the function returns an integer. An integer, not an array of integers. I don't think the standard allows function declarations such as:

integer, dimension(10) function polynomialMult(npts,x,y)

but I could be wrong. I always use the form I showed you above.

If you have an up to date Fortran compiler you can do clever things such as return an allocated array. And I suggest you figure out array syntax. For example, your statement:

polynomialMult =  x(1:npts) + 1

could more concisely be written:

polynomialMult =  x + 1

since Fortran will map the scalar addition to all elements of the array x which you have declared to have only npts elements.

Passing the sizes of arrays into subroutines is very FORTRAN77 and almost always unnecessary now. Generally you either want to operate on every element in an array (as in the array syntax example) or you should let the subprogram figure out the size of the array it is dealing with.

Solution 2

I agree with the previous responder that the following works:

polynomialMult = x + 1

However, without knowing that polynomialMult and x are arrays, one might assume it is a scalar operation. I prefer to be obvious and do it this way:

polynomialMult(:) = x(:) + 1

I have even insisted that the coders in my group do it this way. I don't like to work hard to understand someone's code--I want it to be obvious what they are doing.

Comments

[High Performance Mark]

Haha, returned to this question when it got another up vote. There is a subtle 'issue' with this answer: polynomialMult is an array and polynomialMult(:) is an array section and the two are not always interchangeable, though they are here.

[astay13]

@jvriesem, it depends on how the indices of polynomialMult are set. Sometimes there are good reasons to use non-standard indices such as double precision polynomialMult(0:npts-1). In that case you risk bugs with polynomialMult(0:npts-1) = x(0:npts-1) + 1 since you might accidently mistype the indices in the expression.