Fortran - expanding an array

Author

expanding an array

Bart Vandewoestyne

2005-02-24, 8:58 am

I have a rank-1 array which I pass to a subroutine, and after the
subroutine the length of the array may be changed. I've written this
little test-example which i think is not quite OK, and I would like to
ask if it possible to point me to the cleanest way of doing this. The
example is:

program test_expand

implicit none

! Should this be a pointer???
integer, allocatable :: p(:)

allocate(p(3))

p = (/ 1, 2, 3 /)
print *, "p = ", p

call expand(p, 2)

! This should give p = 1 2 3 4 5
print *, "p = ", p

deallocate(p)

contains

subroutine expand(p, n)
integer, allocatable, intent(inout) :: p(:)
integer, intent(in) :: n

integer, allocatable :: x(:)
integer :: k
integer :: i

k = size(p, 1)
allocate(x(k+n))
x(1:k) = p
x(k+1:n) = (/ (k+i, i=1,n) /)
deallocate(p)
allocate(p(k+n))
p = x
deallocate(x)

end subroutine expand

end program test_expand

If I compile this with NagWare's f95 i get this message:

Extension: test_expand.f95, line 23: ALLOCATABLE dummy argument P
detected at ::@P
[f95 continuing despite warning messages]

but the program seems to work:

$ ./a.out
p = 1 2 3
p = 1 2 3 4 5

If I compile it with g95 i get wrong results:

$ g95 test_expand.f95
$ ./a.out
p = 1 2 3
p = 1 2 3 0 0

Is it cleaner to declare p as a pointer in my main program and pass that
pointer as an argument to my subroutine? And then allocate new storage in my
subroutine and make p point to that newly allocated space?

Regards,
Bart

--
"Share what you know. Learn what you don't."

Michael Metcalf

2005-02-24, 8:58 am

"Bart Vandewoestyne" <MyFirstName.MyLastName@telenet.be> wrote in message
news:1109241573.904867@seven.kulnet.kuleuven.ac.be...
> I have a rank-1 array which I pass to a subroutine, and after the

> x(k+1:n) = (/ (k+i, i=1,n) /)

This line should be

x(k+1:k+n) = (/ (k+i, i=1,n) /)

That aside, the use of allocatable arguments is a standarized extension to
f95 ("Fortran 90/95 Explained", Ch. 13), so you can use them, rather than
pointers, as long as your compilers allow them.

Regards,

Mike Metcalf

Bart Vandewoestyne

2005-02-24, 8:58 am

In article <cvkdnk$qhe$05$1@news.t-online.com>, Michael Metcalf wrote:
>
> This line should be
>
> x(k+1:k+n) = (/ (k+i, i=1,n) /)

Yes, sorry for that... i also spotted it after i did my post...

> That aside, the use of allocatable arguments is a standarized extension to
> f95 ("Fortran 90/95 Explained", Ch. 13), so you can use them, rather than
> pointers, as long as your compilers allow them.

OK, so the way I used them (with the above little mistake corrected) is
clean and ok?

I also tried to create a solution with pointers, this is what I've done:

program test_expand

implicit none

! I now use a pointer here
integer, pointer :: p(:)

allocate(p(3))

p = (/ 1, 2, 3 /)
print *, "p = ", p

call expand(p, 2)

! This should give p = 1 2 3 4 5
print *, "p = ", p

deallocate(p)

contains

subroutine expand(p, n)
integer, pointer :: p(:)
integer, intent(in) :: n

integer, pointer :: x(:)
integer :: k
integer :: i

k = size(p, 1)
allocate(x(k+n))
x(1:k) = p
x(k+1:k+n) = (/ (k+i, i=1,n) /)
deallocate(p)
allocate(p(k+n))
p = x
deallocate(x)

end subroutine expand

end program test_expand

This seems to compile and run without warnings/errors with both NagWare's f95
and g95.

Now my question:

What is the preferable way? The pointer solution of this post, or the solution with the allocatable array?

My guess would be the solution with pointers, since that solution doesn't
require 'a standardized extension to f95', so all compilers should be able to
support that solution.

For the record: i prefer my code to compile with as many compilers as possible, so i prefer not to make use of extensions that might not be supported by all
compilers. This would mean my best choice is the pointer-solution, right?

Regards,
Bart

--
"Share what you know. Learn what you don't."

Michael Metcalf

2005-02-24, 3:59 pm

"Bart Vandewoestyne" <MyFirstName.MyLastName@telenet.be> wrote in message
news:1109246757.70825@seven.kulnet.kuleuven.ac.be...
>
> For the record: i prefer my code to compile with as many compilers as
possible, so i >prefer not to make use of extensions that might not be
supported by all
> compilers. This would mean my best choice is the pointer-solution, right?

Right, by your definition of best. Of, course, it is an abuse of pointers,
but there is no other alternative if you want to keep highly portable (a
laudable aim and one that the standard is supposed to foster). Note that
your array 'x' is better defined as allocatable rather than pointer, as it
local to expand and you are not using any other pointer properties.

Regards,

Mike Metcalf

David Frank

2005-02-24, 3:59 pm

"Michael Metcalf" <metcalfmetcalf@compuserve.com> wrote in message
news:cvkhi9$pdb$00$1@news.t-online.com...
>
> "Bart Vandewoestyne" <MyFirstName.MyLastName@telenet.be> wrote in message
> news:1109246757.70825@seven.kulnet.kuleuven.ac.be...
> possible, so i >prefer not to make use of extensions that might not be
> supported by all
>
> Right, by your definition of best. Of, course, it is an abuse of pointers,
> but there is no other alternative if you want to keep highly portable (a
> laudable aim and one that the standard is supposed to foster). Note that
> your array 'x' is better defined as allocatable rather than pointer, as it
> local to expand and you are not using any other pointer properties.
>
> Regards,
>
> Mike Metcalf
>
>

Whats the shortest syntax to do this F2003?

Bart Vandewoestyne

2005-02-24, 3:59 pm

In article <cvkhi9$pdb$00$1@news.t-online.com>, Michael Metcalf wrote:
>
> [...]
> Note that
> your array 'x' is better defined as allocatable rather than pointer, as it
> local to expand and you are not using any other pointer properties.

I've changed x to an allocatable instead of a pointer, and it also
works, but can you explain a bit more why it is better to use an
allocatable for x instead of a pointer?

I do understand that x is only a local for the subroutine, and that I am
not using any pointer properties for x, but I don't understand why an
allocatable is preferable to a pointer. Is it a performance issue?

Regards,
Bart

--
"Share what you know. Learn what you don't."

Richard E Maine

2005-02-24, 3:59 pm

In article <1109253026.269514@seven.kulnet.kuleuven.ac.be>,
Bart Vandewoestyne <MyFirstName.MyLastName@telenet.be> wrote:

> I do understand that x is only a local for the subroutine, and that I am
> not using any pointer properties for x, but I don't understand why an
> allocatable is preferable to a pointer. Is it a performance issue?

Two issues.

One is performance, yes. Pointers are more flexible than allocatables in
many ways. Compilers have to allow for that flexibility, which often
causes slowdowns. The possibility of different pointers aliasing
(pointing to the same or overlapping storage) can be a big performance
hit. The possibility of noncontiguity can be another. Both those can
happen in pointers, but not in allocatables.

The other issue is user errors. Pointers are just plain error-prone.
Fortran pointers are designed to avoid some of the kinds of errors that
you get in languages such as C, but there still plenty of ways to go
wrong. A simple, localized use like this isn't likely to run into much
in the way of that kind of problem, but as a general statement, the use
of pointers adds substantially to the odds of programming errors. We
quite regularly see such errors here on clf - sometimes caused by
misunderstanding of how pointers work, and sometimes caused by mistakes
that even experts make.

Allocatables, on the other hand, tend to work pretty much like many
people intuitively expect. This greatly lowers the error rate. On the
other hand, some things about allocatables are hard for compilers to get
right. They compiler bug frequency for allocatable stuff is
significantly higher than for other more mature features. That statement
is based on general observation rather than collecting hard statistics,
I admit. I expect compilers to mature in this area eventually, but it is
a slow process.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Dick Hendrickson

2005-02-24, 3:59 pm

David Frank wrote:
> "Michael Metcalf" <metcalfmetcalf@compuserve.com> wrote in message
> news:cvkhi9$pdb$00$1@news.t-online.com...
>
>
>
> Whats the shortest syntax to do this F2003?
>
>
I probably should think about this more, but if "p" is
an allocatable array, then I believe
p = [p,4,5]
will do what the OP wanted.

F2003 calls for deallocation of the LHS and reallocation to
the right stuff if it differs from the size, etc., of the
RHS. "p" will become an array of size 5 with the original
3 values and 4 and 5 as new values.

Be careful though. The statement
p(:) = [p,4,5]
doesn't do the magic de/re allocation.

Dick Hendrickson

Michael Metcalf

2005-02-24, 4:00 pm

"David Frank" <dave_frank@hotmail.com> wrote in message
news:421dd680$0$60068$ec3e2dad@news.usenetmonster.com...
>
>
> Whats the shortest syntax to do this F2003?
>
>

Here's an example of expanding an array from Section 17.5.3 of "Fortran
95/200 Explained":

real, allocatable :: a(:,:), temp(:,:)
:
! Increase size of a to (n,m)
allocate(temp(n,m))
temp(1:size(a,1),1:size(a,2)) = a
call move_alloc(temp,a)
! a is now (n,m), and temp is unallocated

Only one copy is required, rather than two at the monent.

Regards,

Mike Metcalf

Michael Metcalf

2005-02-24, 4:00 pm

"Bart Vandewoestyne" <MyFirstName.MyLastName@telenet.be> wrote in message
news:1109253026.269514@seven.kulnet.kuleuven.ac.be...
>
> I do understand that x is only a local for the subroutine, and that I am
> not using any pointer properties for x, but I don't understand why an
> allocatable is preferable to a pointer. Is it a performance issue?
>
In general, yes. Defining a pointer sets up a descriptor that includes the
address of its target, if any. When you reference a pointer, there is an
automatic dereference 'behind the scenes'. All that starts to become an
overhead in storage and execution if you don't need any properties of a
pointer other than the fact that it can be dynamically allocated.

Regards,

Mike Metcalf

Sponsored Links