Variable-length array output :D 2

[NickFort]

Hi all,

I'm sure this is obvious to many of you, but it wasn't to me! I just had a little revelation regarding variable-length arrays as the output of subroutines/functions, and I have to share this with those of you who don't know, because I can't contain my excitement.

I was just sitting here wondering how one would have a variable-length array (vector in my case) as an output to a function. I knew allocatable arrays could be allocated to different sizes, but what I didn't know was that if the dummy argument of a subroutine was allocatable and allocated within the subroutine, the argument of the actual procedure call (which also has the "allocatable" attribute) gets allocated accordingly.

Now, I don't recall reading this in any of the Fortran books I've read, but I must admit that I tend to skip to relevant sections and skim over others. I'd been wondering for aaaaaages how it could be done, and I decided to play around today, expecting the worst. This post is here for the benefit of those who -- like me before about 20 minutes ago -- don't know how to do it.

To this demonstrate how this works is the code below. It simply generates a random-length vector, whose elements are random reals, sends it to the driver program, prints the result, and finishes. It works on gfortran, but not on g95 (which doesn't seem to have a working "system_clock" feature to generate the "random" seed). If you want to use in with g95, you'll have to use a different method for the random seed.

program var_length
implicit none

interface
    subroutine random_vector(vector)
    implicit none
    real, dimension(:), allocatable, intent(out) :: vector
    end subroutine
end interface

integer :: length_call,k
real, dimension(:), allocatable :: vector_call

call random_vector(vector_call)

length_call = size(vector_call)

print *, "length_call: ", length_call
do k=1,size(vector_call)
    print *, vector_call(k)
end do

deallocate(vector_call)

end program var_length

!==============================================

subroutine random_vector(vector)
implicit none
real, dimension(:), allocatable, intent(out) :: vector
integer :: length
real :: length_real

if (allocated(vector)) deallocate(vector)
call init_random_seed()
call random_number(length_real)
length = int(length_real*10)+1
allocate (vector(length))
call random_number(vector)

end subroutine random_vector

!==============================================

subroutine init_random_seed()
    integer :: i, n, clock
    integer, dimension(:), allocatable :: seed
    call random_seed(size = n)
    allocate(seed(n))
    call system_clock(count=clock)
    seed = clock + 37 * (/ (i - 1, i = 1, n) /)
    call random_seed(put = seed)
    deallocate(seed)
end subroutine

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[xwb]

Note that this will only work from F95 upwards. It won't work on F90.

 

[NickFort]

Hmmm, interesting. Practically, that's not really an issue nowadays, is it, as basically all F90 compilers will support at least F95?

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[xwb]

Not really: some older compilers are F90 only. Took me ages to figure that out. It really depends: some sites have a Fortran compiler that was paid for on some budget or other and that is the only one they know. Someone asked about this about 4 years back. I submitted an example similar to the one above but he couldn't get it to build or execute. Took a while to figure out that he was using an F90 instead of F95 compiler.

 

[NickFort]

So how was it done pre-F95? Was there a way of doing it at all? Or did you use a big array, populate only what was needed, and output an integer telling you how many elements of the array were part of the desired output? I can think of a few other ways of doing it, but all seem very convoluted...

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[xwb]

Pre-95, it couldn't be done within a routine. All the allocatable stuff had to be done outside the calling routine. You basically did one call to find out how big it was meant to be, passed it back to the caller, allocated the array and passed it in.

 

[NickFort]

Yeah, that's one of the "I can think of a few other ways of doing it, but all seem very convoluted..." I mentioned above.

In fact, that's exactly what I was trying to do, when I compiled, got the right result, and then realised I hadn't actually allocated the array explicitly in the calling program. That's how I stumbled upon this little gem (and I do consider it a gem, because it really makes life easier!).

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[xwb]

The great thing about making mistakes like these is that you never forget them.

 

[GerritGroot]

Thanks Nick,

I've been dealing with this problem many times and didn't know this solution existed. My previous commercial compiler didn't let me either, but in spite of the clumsy debugging I switched to gfortran since a year, so this will do.

I also used the solution xwb mentioned, calling the array-filling routine twice with some flag, but it's really annoying to program that.

Before, I also have been using modules to do this and I made a separate routine inside the same module as where the array was allocated and called it "FreeMyStuff" or so to DEALLOCATE everything again.

 

[NickFort]

I'm glad this helped someone! I think it's a fantastically useful addition to Fortran.

Regarding gfortran and clunky debugging, I agree with you on that. That's why I use both g95 and gfortran for debugging -- they often give different error messages, and one will be more useful than the other.

At some point, I want to write a detailed comparison of gfortran and g95. I've run into quite a few problems with the latter, despite its better error messages and more advanced implementation of F2003 (or at least, that was the case last time I checked.) But that's another matter altogether, and doesn't belong in this thread.

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[GerritGroot]

I usually include the option -Wall in gfortran.

http://gcc.gnu.org/onlinedocs/gcc-4.1.2/gfortran/Warning-Options.html#Warning-Options

 

[NickFort]

Yeah, I should probably use -Wall a little more often than I do at the moment; it's a good idea.

Regarding my claim earlier regarding g95:

despite its better error messages and more advanced implementation of F2003 (or at least, that was the case last time I checked.)

It turns out that I'm wrong; according this this document:

http://www.fortranplus.co.uk/resources/fortran_2003_2008_compiler_support.pdf

gfortran is winning on both the F2003 and F2008 fronts the last time that was updated.

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[mikrom]

Could it be possible to use for this purpose a function instead of subroutine?
I mean a function without argument which would return random vector of random length.

 

[NickFort]

I've just tried to make a function that does this, but at the function call, there is an "Exception: Access Violation" according to g95, suggesting that the "vector_call = random_vector()" requires "vector_call" to be explicity allocated.

I would have expected that to work, but it doesn't; maybe someone else can shed some light on this.

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[NickFort]

Hmmm... It's supposed to be a F2003 feature. In the gfortran 4.6 manual, we have:

TR 15581:

- ALLOCATABLE dummy arguments
- ALLOCATABLE function results
- ALLOCATABLE components of derived types

All three of those are useful, and supposedly implemented by gfortran.

However, it still crashes, even if I force -std=f2003.

I'm about to flip through "Fortran 95/2003 Explained", to see if Messrs Metcalf, Reid and Cohen have any answers for us...

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[NickFort]

After looking through "Fortran 95/2003 Explained", "Fortran 95/2003 for Scientists and Engineers" and "The Fortran 2003 Handbook", I've come to the following conclusion, although it's not stated specifically in any of those books:

The return value of a function can be allocatable. However, what it's assigned to cannot. So, for example, you can define

.
.
real, dimension(100) :: vector_call
.
.
vector_call(:size(random_number()) = random_number()

However, I find this a bit clumsy and limiting myself. If "vector_call" could be allocatable as well, it would be a lot more convenient.

If I've missed something and it can be done, I would love to stand corrected.

Does anyone know if there are plans to implement this in the F2008 standard if it isn't in the F2003 one? If there are, we'll just have to wait until 2020 or so until the full F2008 is implemented in most compilers...

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[NickFort]

Oops, there's a missing parenthesis in the above code. It should be:

vector_call(:size(random_number())) = random_number()

Apologies!

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[NickFort]

xwb, regarding my original post: it appears that this is a F2003 feature, not F95. The reasons I say this:

(1) "Allocatable dummy arguments" are included in the F2003 part of Metcalf et al.'s "Fortran 95/2003 Explained".
(2) If I compile the code with gfortran's -std=f95 to force the F95 standard, I get the error message "Error: Fortran 2003: ALLOCATABLE attribute with DUMMY attribute"

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[NickFort]

OK, to settle it once and for all. Straight from Tobias Burnus:

Allocatable function results and dummy arguments indeed work - however,
you use a related feature, which is not yet supported:

real, dimension), allocatable :: vector_call
vector_call = random_vector()


Here, vector_call is unallocated and is supposed to get allocated upon
assignment. That's allowed in Fortran 2003 but not yet implemented in
gfortran. Cf. "Assignment to an allocatable array" at the F2003 status at

-

http://fortranwiki.org/fortran/show/Fortran+2003+status

-

http://www.fortranplus.co.uk/resources/fortran_2003_2008_compiler_support.pdf

-

http://gcc.gnu.org/wiki/Fortran2003Status

(if you use other compilers: note, that some require some special option
to enable (re)allocate on assignment.)

and the description of the items at
- ftp://ftp.nag.co.uk/sc22wg5/N1601-N1650/N1648.pdf

The gfortran bug is:
-

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=35810

Like with all software development, especially with volunteer-based
ones, it is difficult to give an estimate when it will be implemented;
to my knowledge no one currently works on this features - but as one of
the few missing F2003 features, its priority raises as other F2003
features get implemented.

So it is in the standard, but hasn't been implemented in either gfortran or g95 yet. I do find it odd that assigning to an allocatable array works from a subroutine; I've asked Tobias why this is, and I'll keep you guys posted.

So, mikrom, to answer your question: it will work one day with open source compilers. Currently, it should work on Cray, Intel, IBM, and NAG.

If anyone has any of these compilers, I'd appreciate it if you could confirm that the following code works:

program var_length
implicit none

interface
    function random_vector() result(vector)
    implicit none
    real, dimension(:), allocatable :: vector
    end function
end interface

integer :: length_call,k
real, dimension(:), allocatable :: vector_call

vector_call = random_vector()

length_call = size(vector_call)

print *, "length_call: ", length_call
do k=1,size(vector_call)
    print *, vector_call(k)
end do

deallocate(vector_call)

end program var_length

!==============================================

function random_vector() result(vector)
implicit none
real, dimension(:), allocatable :: vector
integer :: length
real :: length_real

if (allocated(vector)) deallocate(vector)
call init_random_seed()
call random_number(length_real)
length = int(length_real*10)+1
allocate (vector(length))
call random_number(vector)

end function random_vector

!==============================================

subroutine init_random_seed()
    integer :: i, n, clock
    integer, dimension(:), allocatable :: seed
    call random_seed(size = n)
    allocate(seed(n))
    call system_clock(count=clock)
    seed = clock + 37 * (/ (i - 1, i = 1, n) /)
    call random_seed(put = seed)
    deallocate(seed)
end subroutine

--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.

 

[mikrom]

Hi NickFort,
Thanks for valuable information. You earned star from me.