Add element to a polymorphic arrar

[\"Vladimír Fuka <"name.surnameat]

Hi,
it is easy to reallocate an add an item to a non-polymorphic array

subroutine add_element_int(a,e)
integer,allocatable,intent(inout) :: a(:)
integer,intent(in) :: e
integer,allocatable :: tmp(:)

if (.not.allocated(a)) then
a = [e]
else
allocate(tmp(size(a)),source=a)
deallocate(a)
allocate(a(size(tmp)+1))
a(1:size(tmp)) = tmp
a(size(tmp)+1) = e
end if
end subroutine

(It is supposed to contain just couple of items and most of the work is
getting them, otherwise I know how to do a stack that grows in powers of
two.)

With a polymorphic array I ended up with this ugly thing (does not
compile in gfortran because of compiler error). Is there some other
possibility that does not use type(*) (which is not yet supported by
ifort)? Among problems I encountered is that c_loc() and c_sizeof() do not
accept non-interoperable arguments and that I cannot do any assignments
even if the dynamic types are the same (check by same_type_as()).

subroutine add_element_poly(a,e)
use iso_c_binding
class(*),allocatable,intent(inout),target :: a(:)
class(*),intent(in),target :: e
class(*),allocatable,target :: tmp(:)
type(c_ptr) :: dummy
integer i

interface
function memcpy(dest,src,n) bind(C,name="memcpy") result(res)
import
type(c_ptr) :: res
integer(c_intptr_t),value :: dest
integer(c_intptr_t),value :: src
integer(c_size_t),value :: n
end function
end interface

if (.not.allocated(a)) then
allocate(a(1), source=e)
else
allocate(tmp(size(a)),source=a)
deallocate(a)
allocate(a(size(tmp)+1),mold=e)
dummy = memcpy(loc(a(1)),loc(tmp),sizeof(tmp))
dummy = memcpy(loc(a(size(tmp)+1)),loc(e),sizeof(e))
end if
end subroutine

[Wolfgang Kilian]

Isn't that what MOVE_ALLOC is for?

-- Wolfgang

--
E-mail: firstnameini...@domain.de
Domain: yahoo

[\"Vladimír Fuka <"name.surnameat]

Yes, it is better
call move_alloc(a,tmp)
instead of
allocate(tmp(size(a)),source=a)
deallocate(a)


but my question is more about the assignment (or moving bit patterns) to a
polymorphic sub-array.

Tato zpráva byla vytvořena převratným poštovním klientem Opery:
http://www.opera.com/mail/

[Wolfgang Kilian]

On 05/17/2013 02:01 PM, \"Vladimír Fuka <\"name.surname at >> wrote:
>

> Hi,
> it is easy to reallocate an add an item to a non-polymorphic array

> [,,,]

>
> With a polymorphic array I ended up with this ugly thing (does not
> compile in gfortran because of compiler error). Is there some other
> possibility that does not use type(*) (which is not yet supported by
> ifort)? Among problems I encountered is that c_loc() and c_sizeof() do
> not accept non-interoperable arguments and that I cannot do any
> assignments even if the dynamic types are the same (check by
> same_type_as()).

> [...]

Yes, the Fortran2003 standard forbids intrinsic assignment with a
polymorphic l.h.s. It is allowed in F2008.

Without F2008 support, the problem is how to ALLOCATE and copy in a
single step.

Here is code that works with NAG Fortran. It doubles the size of a
polymorphic CLASS(*) array and copies the values to a temporary. The
copy back to the array is then avoided by using MOVE_ALLOC, but there is
probably still a redundant copy involved.

program main
implicit none

class(*), dimension(:), allocatable :: a, tmp
integer :: newsize

allocate (integer :: a (5))
select type (a)
type is (integer); a = [1, 2, 3, 4, 5]
end select

! Expand array:
newsize = 2 * size (a)
allocate (tmp (newsize), &
source = reshape (spread (a, 2, 2), [newsize]))
call move_alloc (tmp, a)
! Now a has doubled its size.

select type (a)
type is (integer); print *, a
end select

end program main


(TRANSFER or UNPACK are other candidates besides SPREAD.)

The experts may answer whether the RESHAPE and SPREAD intrinsics are
allowed by the Standard within the source expression of an ALLOCATE
statement. We recently got a response from Intel support that they are
not (because they allocate an anonymous array behind the scenes), but
I'm not convinced. Maybe the above example is standard-conforming.

[\"Vladimír Fuka <"name.surnameat]

> Yes, the Fortran2003 standard forbids intrinsic assignment with a
> polymorphic l.h.s. It is allowed in F2008.
>
> Without F2008 support, the problem is how to ALLOCATE and copy in a
> single step.

I wouldn't mind F2008 but I do mind that my compilers do not support this
feature.

> [snip]

Interesting, but the problem, hot get the new value there remains.

> (TRANSFER or UNPACK are other candidates besides SPREAD.)
>
> The experts may answer whether the RESHAPE and SPREAD intrinsics are
> allowed by the Standard within the source expression of an ALLOCATE
> statement. We recently got a response from Intel support that they are
> not (because they allocate an anonymous array behind the scenes), but
> I'm not convinced. Maybe the above example is standard-conforming.

Would be nice to know. Maybe this is why i had problems with `allocate(a,
source = [e])` (crashes gfortran 4.8).

[Richard Maine]

Wolfgang Kilian <kil...@invalid.com> wrote:

> The experts may answer whether the RESHAPE and SPREAD intrinsics are
> allowed by the Standard within the source expression of an ALLOCATE
> statement. We recently got a response from Intel support that they are
> not (because they allocate an anonymous array behind the scenes), but
> I'm not convinced. Maybe the above example is standard-conforming.

I don't understand why they wouldn't be allowed. Hmm. On checking the
standard, I suppose I see how one could read it that way, but I think it
would be a misreading. I do see that the words are a bit ambiguous, but
there just would be no reason for a restriction quite like that. The
standard (well f2003) does say

"Neither stat-variable, source-expr, nor errmsg-variable shall be
allocated within the ALLOCATE statement within which it appears..."

but that's really talking about restricting them from being
allocate-objects. This restriction is to avoid some order-of-execution
problems and circular dependencies. Just having a temporary array
involved in one of the expressions isn't the same sort of thing at all.
Yes, that temporary array will have to be allocated and deallocated by
the compiler, but that's just like what happens with temporary arrays in
any expression. It would just seem arbitrary to disallow thing like that
here.

I suppose I can see misreading the words that way, but I can't see them
as being intended that way.

--
Richard Maine
email: last name at domain . net
domain: summer-triangle

[Ian Harvey]

On 2013-05-17 10:01 PM, \"Vladimír Fuka <\"name.surname at >> wrote:
>
> Hi,
> it is easy to reallocate an add an item to a non-polymorphic array
>
> subroutine add_element_int(a,e)
> integer,allocatable,intent(inout) :: a(:)
> integer,intent(in) :: e
> integer,allocatable :: tmp(:)
>
> if (.not.allocated(a)) then
> a = [e]
> else

If you are lazy like me, then:

a = [ a, e ]

is less typing and clearer from a programmer intent point of view (but
likely to have much more execution overhead with current compilers).

Note you are not just talking a polymorphic array - you've got an
unlimited polymorphic array. Consequently:

Perhaps I'm off base here, but I think your intent is to have a generic
procedure (in the sense of "applicable to any type") that appends to an
array, where the thing being appended has the same dynamic type as the
array being appended to. This requirement is potentially compile time
checkable in the non-polymorphic case, but is runtime implied with your
memcpy approach. At the time the call sites of you procedure are
compiled, really all you want to do is "replace" INTEGER with whatever
the relevant dynamic type is for the arguments at that call site.

Basically - you are trying to do generic programming.

Fortran's processor based support for generic programming is poor -
though as often posted on c.l.f. there are source code "hacks" using
INCLUDE etc that help you around this. While these hacks are messy from
a source code appearance point of view, they will probably give you a
much better outcome from the point of view of the code that the compiler
executes, and potentially things like compiler error checking, etc.

(And then you can join the chorus of people begging for better
in-language support of generic programming.)

This is probably further off-base, but a similar situation could be
where you want to have an array of things that could be of any dynamic
type, and you want to append a new object (again of any dynamic type) to
it. In that case, the array is of container type:

TYPE :: Container
CLASS(*), ALLOCATABLE :: item
END TYPE Container

and then append operation is then on a *non-polymorphic* array dummy of
that type:

SUBROUTINE add_element_(a, e)
TYPE(Container), INTENT(INOUT), ALLOCATABLE :: a(:)
CLASS(*), INTENT(IN) :: e
IF (.NOT. ALLOCATED(a)) THEN
a = [ Container(e) ]
ELSE
a = [ a, Container(e) ] ! Being lazy again.
END IF
END SUBROUTINE add_element_

(The definition of a unlimited polymorphic component via a structure
constructor might be erroneous in F2003 (not sure, ridiculous
restriction if so), but there are a couple of ways of working around
this if that is the case and I think the use is ok in F2008.)

If all elements are of the same dynamic type, then you have equivalent
storage (but not exactly the same "use") as your polymorphic array
growing thing. Client code would need to be doing lots of SELECT TYPE
on the array elements in order to do anything useful with the stored
objects, but presumably client code knows what it has been stuffing
inside the array.

(If you are using ifort then this thread

http://software.intel.com/en-us/forums/topic/326077

identifies some potential compiler issues with this approach (see in
particular the .f90 that appear as attachments later on), that may have
been fixed by now. In the source examples in that thread the items are
just polymorphic, not unlimited polymorphic, but many aspects are common.)

[Ian Harvey]

On 2013-05-17 11:51 PM, Wolfgang Kilian wrote:
...

> The experts may answer whether the RESHAPE and SPREAD intrinsics are
> allowed by the Standard within the source expression of an ALLOCATE
> statement. We recently got a response from Intel support that they are
> not (because they allocate an anonymous array behind the scenes), but
> I'm not convinced. Maybe the above example is standard-conforming.

If that was their response then that is absurd - to paraphrase: 'the
standard doesn't allow it because our implementation doesn't support it'.

I could imagine a response 'we don't currently support that F2003
feature because of our current implementation'. That's completely
reasonable. There are several F2003 features that fall into that sort
of category for that compiler. Given they've stated that they'll get to
full F2003 eventually, presumably eventually they'll adapt their
implementation.

Another reasonable response would be "the standard doesn't allow it -
see this paragraph or this interpretation'.

[\"Vladimír Fuka <"name.surnameat]

Dne Mon, 20 May 2013 03:40:41 +0200 Ian Harvey <ian_h...@bigpond.com>
napsal(a):

Thanks, I especially like

> If you are lazy like me, then:
>
> a = [ a, e ]

which I was not sure it is allowed (I tried an implied do instead of the
array), but as Fortran does not support arrays of arrays it can work.
Unfortunately I have some problems at runtime with this. When I tried that
in gfortran-4.8 the array is full of garbage
and valgrind gives warnings like:

Invalid read of size 1
==9845== at 0x4C2DE21: memcpy@@GLIBC_2.14 (in
/usr/lib64/valgrind/vgpreload_memcheck-amd64-linux.so)
==9845== by 0x4A3E2E: __geometricshapes_MOD_arrayfromobst
(geometricshapes.f90:1776)
==9845== by 0x831BEA: customsolidbodies_ (in
/home/lada/f/CLMM-work/bodies/GAUK/CLMM)
==9845== by 0x4BB0B9: __solidbodies_MOD_initsolidbodies
(solidbodies.f90:182)
==9845== by 0x693D5E: __initial_MOD_initboundaryconditions
(initial.f90:1869)
==9845== by 0x8311C2: MAIN__ (main.f90:28)
==9845== by 0x831B20: main (main.f90:3)
==9845== Address 0x91aa47b is 363 bytes inside a block of size 364 free'd
==9845== at 0x4C2AF6C: free (in
/usr/lib64/valgrind/vgpreload_memcheck-amd64-linux.so)
==9845== by 0x4A3B02: __geometricshapes_MOD_arrayfromobst
(geometricshapes.f90:1776)
==9845== by 0x831BEA: customsolidbodies_ (in
/home/lada/f/CLMM-work/bodies/GAUK/CLMM)
==9845== by 0x4BB0B9: __solidbodies_MOD_initsolidbodies
(solidbodies.f90:182)
==9845== by 0x693D5E: __initial_MOD_initboundaryconditions
(initial.f90:1869)
==9845== by 0x8311C2: MAIN__ (main.f90:28)
==9845== by 0x831B20: main (main.f90:3)


1776 is the line above. I do not know if it is a known limitation or not.


Vladimir

[Ian Harvey]

On 2013-05-21 6:26 PM, \"Vladimír Fuka <\"name.surname at >> wrote:
> Dne Mon, 20 May 2013 03:40:41 +0200 Ian Harvey <ian_h...@bigpond.com>
> napsal(a):
>
> Thanks, I especially like
>
>> If you are lazy like me, then:
>>
>> a = [ a, e ]
>
>
> which I was not sure it is allowed (I tried an implied do instead of the
> array), but as Fortran does not support arrays of arrays it can work.
> Unfortunately I have some problems at runtime with this. When I tried
> that in gfortran-4.8 the array is full of garbage

...

I tried a little test program here with 4.6.3 (supplied with distro) and
4.8.0 (built by myself, might have the odd hack) and saw no issues.

My little test program was:

PROGRAM array_concat
IMPLICIT NONE
INTEGER, ALLOCATABLE :: array(:)
!***
array = [ 1, 2 ]
CALL add_element(array, 3)
PRINT *, array
CONTAINS
SUBROUTINE add_element(a, e)
INTEGER, INTENT(INOUT), ALLOCATABLE :: a(:)
INTEGER, INTENT(IN) :: e
!***
a = [ a, e ]
END SUBROUTINE add_element
END PROGRAM array_concat

(The subroutine assumes the a argument is always allocated.)

[Richard Maine]

Ian Harvey <ian_h...@bigpond.com> wrote:

> SUBROUTINE add_element(a, e)
> INTEGER, INTENT(INOUT), ALLOCATABLE :: a(:)
> INTEGER, INTENT(IN) :: e
> !***
> a = [ a, e ]
> END SUBROUTINE add_element
> END PROGRAM array_concat
>
> (The subroutine assumes the a argument is always allocated.)

Yes. That is important, as the a = [a, e] is invalid if "a" is not
allocated on entry. Note that the OP's code did not have that
assumption, treating the unallocated case separately. In fact, it wasn't
entirely clear to me whether that might be related to the error that he
reported in the revised version.

[Wolfgang Kilian]

To be fair, their response was of that type. But see Richard's comment
above.

[\"Vladimír Fuka <"name.surnameat]

> I tried a little test program here with 4.6.3 (supplied with distro) and
> 4.8.0 (built by myself, might have the odd hack) and saw no issues.

See http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57354 , I hope it is a
valid code.

Vladimir