Smart questions
Smart answers
Smart people
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Member Login




Remember Me
Forgot Password?
Join Us!

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips now!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

Join Tek-Tips
*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.
Jobs from Indeed

Link To This Forum!

Partner Button
Add Stickiness To Your Site By Linking To This Professionally Managed Technical Forum.
Just copy and paste the
code below into your site.

Variable-length array output :DHelpful Member!(2) 

NickFort (TechnicalUser) (OP)
23 Jul 10 8:46
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.

CODE

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 (Programmer)
23 Jul 10 16:54
Note that this will only work from F95 upwards.  It won't work on F90.
Helpful Member!(2)  NickFort (TechnicalUser) (OP)
24 Jul 10 9:04
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 (Programmer)
24 Jul 10 13:47
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 (TechnicalUser) (OP)
24 Jul 10 14:02
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 (Programmer)
24 Jul 10 16:32
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 (TechnicalUser) (OP)
25 Jul 10 4:49
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 (Programmer)
25 Jul 10 16:30
The great thing about making mistakes like these is that you never forget them.
GerritGroot (TechnicalUser)
27 Jul 10 4:03
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 (TechnicalUser) (OP)
27 Jul 10 4:23
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 (TechnicalUser)
27 Jul 10 6:20
NickFort (TechnicalUser) (OP)
27 Jul 10 12:18
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:

Quote:

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 (Programmer)
27 Jul 10 13:06
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 (TechnicalUser) (OP)
27 Jul 10 13:26
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 (TechnicalUser) (OP)
27 Jul 10 13:42
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 (TechnicalUser) (OP)
27 Jul 10 14:27
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

CODE

.
.
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 (TechnicalUser) (OP)
27 Jul 10 14:29
Oops, there's a missing parenthesis in the above code. It should be:

CODE

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

Apologies!

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

NickFort (TechnicalUser) (OP)
27 Jul 10 14:42
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 (TechnicalUser) (OP)
28 Jul 10 4:12
OK, to settle it once and for all. Straight from Tobias Burnus:

Quote:

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:

CODE

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 (Programmer)
28 Jul 10 15:15
Hi NickFort,
Thanks for valuable information. You earned star from me.
NickFort (TechnicalUser) (OP)
28 Jul 10 15:51
Thanks, mikrom, I appreciate it. Glad I could help!

While I'm here, an explanation from Tobias regarding why the subroutine works, but the function doesn't:

Quote:

in the main program:

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


Here, you pass (effectively) the address of the array descriptor to the
subroutine "random_vector" - and that variable gets evaluated.

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


Whereas in the function case, you allocate the result variable "vector"
- but not the variable "vector_call". Fortran works such that first all
on the RHS of an assignment is evaluated and then the result assigned to
the LHS. Thus, the RHS is allocated, but gfortran assumes (as does
Fortran 90/95) that the LHS is allocated with the correct size - and
does not touch the allocation status or the size (rather: shape) of the LHS.

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

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Back To Forum

Close Box

Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close