Pointer-valued function to access inner components
Salvatore
Is the following code legal? Specifically, I am referring to the
pointer assignment in the main program. I am getting three different
behaviours from three different compilers....
Thanks
Salvatore Filippone
module class_mesh
type mesh
integer :: n
real(kind(1.d0)), allocatable :: area(:)
end type mesh
contains
subroutine create_mesh(msh)
type(mesh), intent(out) :: msh
msh%n=10
allocate(msh%area(msh%n))
return
end subroutine create_mesh
end module class_mesh
module class_field
use class_mesh
implicit none
private ! Default
public :: create_field, field
public :: msh_
type field
private
character(len=32) :: name
type(mesh), pointer :: msh => null()
integer :: isize(2)
end type field
interface msh_
module procedure msh_
end interface
interface create_field
module procedure create_field
end interface
contains
subroutine create_field(fld,msh)
type(field), intent(out) :: fld
type(mesh), intent(in), target :: msh
fld%msh => msh
fld%name = 'Temp'
fld%isize = 1
end subroutine create_field
function msh_(fld)
type(mesh), pointer :: msh_
type(field), intent(in) :: fld
msh_ => fld%msh
end function msh_
end module class_field
module class_scalar_field
use class_field
implicit none
private ! Default
public :: create_field, scalar_field
public :: msh_
type scalar_field
private
type(field) :: base
real(kind(1.d0)), allocatable :: x(:)
end type scalar_field
interface create_field
module procedure create_scalar_field
end interface
interface msh_
module procedure get_scalar_field_msh
end interface
contains
subroutine create_scalar_field(fld,msh)
use class_mesh
! Mandatory arguments
type(scalar_field), intent(out) :: fld
type(mesh), intent(in), target :: msh
call create_field(fld%base,msh)
allocate(fld%x(10))
end subroutine create_scalar_field
function get_scalar_field_msh(fld)
use class_mesh
!
type(mesh), pointer :: get_scalar_field_msh
type(scalar_field), intent(in), target :: fld
get_scalar_field_msh => msh_(fld%base)
end function get_scalar_field_msh
end module class_scalar_field
program test_pnt
use class_mesh
use class_scalar_field
implicit none
!
type(mesh) :: msh
type(mesh), pointer :: mshp
type(scalar_field) :: quality
call create_mesh(msh)
call create_field(quality,msh)
mshp => msh_(quality)
end program test_pnt
Richard Maine
> Is the following code legal? Specifically, I am referring to the
> pointer assignment in the main program. I am getting three different
> behaviours from three different compilers....
The pointer assignment is ok. I need to leave pretty soon and don't have
time to adequately peruse the rest of the code.
However, just reading the subject line rings alarm bells. I truly wish
that the standard did not allow pointer-values functions. It does allow
them, but I consider their problems to be far more severe than their
benefits. Pointer-valued functions are extremely error prone, even for
expert users.
If you adopt a style that keeps away from the most problematic
error-prone areas, then you pretty much give up what functions are for.
The thing that most interestingly distinguishes a function from a
subroutne is that you can evaluate a function as part of an expression.
You aren't restricted to just things like y=f(x); you can also have
y=f(x)+42, or whatever. But do that with a pointer-valued function and
you probably introduced a memory leaks or other bug.
Yes, there are isolated cases where you can avoid the problems. But to
me, those cases are too isolated to counter the error-proneness of teh
whole thing.
I do note that you've got generics that have specifics of the same name
as the generic. The standard does allow that, but it wouldn't shock me
if some compilers get as confused by it as I do. (I didn't see at first
that was what was going on).
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
Salvatore
Well, the whole idea of this function, as proposed by a coworker, was
to get acces (and only read access at that) via a pointer to an inner
vector which is otherwise hidden (avoiding a copy). Indeed the whole
thing works if I change the function+assignment to a subroutine call.
However I got a whole range of behaviour from different compilers,
from working code to segfault at runtime to compile time error; even
if the code is dangerous, and I can agree it should be changed, I
wanted to know whether I had a right to complain against the compiler.
Thanks
Salvatore
mm148881
If you use a different name for the specific function msh_ in the
class_field module, for instance mshb_, your code compiles on g95 (v
0.90), intel (v 9.1) and pgf95 (v 6.2) and gives the expected results
for the three compilers. I tried g95 and ifort on powerpc and intel
based iMacs, then pgf90 on an opteron.
g95 seems to be getting really confused by the specific function name
being identical to the generic one.
Sincerely
Massimo
glen herrmannsfeldt
> On Feb 14, 11:00 pm, nos...@see.signature (Richard Maine) wrote:
(snip on pointer value functions)
>>If you adopt a style that keeps away from the most problematic
>>error-prone areas, then you pretty much give up what functions are for.
>>The thing that most interestingly distinguishes a function from a
>>subroutne is that you can evaluate a function as part of an expression.
>>You aren't restricted to just things like y=f(x); you can also have
>>y=f(x)+42, or whatever. But do that with a pointer-valued function and
>>you probably introduced a memory leaks or other bug.
Pointer valued functions are a problem in C, but one must live
with them, anyway. There are a few possibilities in C:
1) Allocating memory and returning a pointer to it. That is the case
that can cause a memory leak, though even there allocation is
normally done by a function, malloc().
2) Returning a pointer to static data. This works well as long
as you copy the result somewhere else before calling the function
again. Many C library routines work this way.
3) Returning a pointer to one of the arguments, or some offset
into one of the arguments. The former is for convenience, the
latter probably doesn't work in Fortran.
>>Yes, there are isolated cases where you can avoid the problems. But to
>>me, those cases are too isolated to counter the error-proneness of teh
>>whole thing.
(snip)
> Well, the whole idea of this function, as proposed by a coworker, was
> to get acces (and only read access at that) via a pointer to an inner
> vector which is otherwise hidden (avoiding a copy). Indeed the whole
> thing works if I change the function+assignment to a subroutine call.
> However I got a whole range of behaviour from different compilers,
> from working code to segfault at runtime to compile time error; even
> if the code is dangerous, and I can agree it should be changed, I
> wanted to know whether I had a right to complain against the compiler.
That sounds like one of the cases that could work. Since you aren't
allocating new memory, there is no memory leak by not releasing it.
Not knowing Fortran pointers as well as C pointers, can you have a
pointer to an array element or structure element? As far as I know,
you can't have a pointer to part of an array, as is common in C.
-- glen
Jan Vorbrüggen
> pointer to an array element or structure element? As far as I know,
> you can't have a pointer to part of an array, as is common in C.
Huh? I thought you can have a pointer refer to any legal slice of an array -
that would include a single element among the possibilities...
Jan
robert....@sun.com
wrote:
> Not knowing Fortran pointers as well as C pointers, can you have a
> pointer to an array element or structure element? As far as I know,
> you can't have a pointer to part of an array, as is common in C.
If you can have a pointer to an array, you can have a pointer
to an element of that array.
Bob Corbett
Salvatore
1. the code compiles cleanly with XLF 10.1 on an AIX, but I've got
reports that older versions of XLF spew out error messages.
2. Intel Fortran 9.1 on LInux digests it OK.
3. GNU Fortran 4.2: a bug has been acknowledged by the developers, and
is being fixed.
Regards
Salvatore FIlippone
glen herrmannsfeldt
> On Feb 15, 12:20 pm, glen herrmannsfeldt wrote:
>>Not knowing Fortran pointers as well as C pointers, can you have a
>>pointer to an array element or structure element? As far as I know,
>>you can't have a pointer to part of an array, as is common in C.
> If you can have a pointer to an array, you can have a pointer
> to an element of that array.
I know Fortran doesn't let you do all the things that C does, I
believe intentionally. I especially believe that Fortran doesn't
have anything like C's p++ to point to the next element of an array.
So, if I have, say:
REAL, TARGET:: A(100)
REAL, POINTER:: B
B => A(I)
Will point to the I'th element of A. Given a pointer to
the I'th element, I can't point it to any other elements.
Maybe not so far off, though:
REAL, TARGET:: A(100)
REAL, POINTER:: B(:)
B => A(I:UBOUND(A))
B points to the subarray starting at I, and
B => B(2:UBOUND(B))
will point to a subarray, similar to B++ in C.
To be more C like, maybe it should be
B(0:) => B(LBOUND(B)+1,UBOUND(B))
A function could then return a pointer to an element or
subarray of one of its arguments without worries about a memory
leak, and that result could be used in an expression.
I will guess that C's B--,
B(0:) => B(LBOUND(B)-1,UBOUND(B))
doesn't work, though.
-- glen
David Frank
"glen herrmannsfeldt" <g...@ugcs.caltech.edu> wrote in message
news:rOedndd57LX4CUjY...@comcast.com...
> robert....@sun.com wrote:
>
>> On Feb 15, 12:20 pm, glen herrmannsfeldt wrote:
>
>>>Not knowing Fortran pointers as well as C pointers, can you have a
>>>pointer to an array element or structure element? As far as I know,
>>>you can't have a pointer to part of an array, as is common in C.
>
>> If you can have a pointer to an array, you can have a pointer
>> to an element of that array.
>
> I know Fortran doesn't let you do all the things that C does, I
> believe intentionally. I especially believe that Fortran doesn't
> have anything like C's p++ to point to the next element of an array.
>
> So, if I have, say:
>
> REAL, TARGET:: A(100)
> REAL, POINTER:: B
>
> B => A(I)
>
> Will point to the I'th element of A. Given a pointer to
> the I'th element, I can't point it to any other elements.
>
> -- glen
Not true, below is acceptable to CVF...
program test
REAL, TARGET:: A(100) = [1:100]
REAL, POINTER:: B
integer :: i = 50
B => A(I)
write (*,*) b ! = 50
b = b+1
write (*,*) b ! = 51 therefore the increment was "intelligent"
end program
David Frank
"David Frank" <dave_...@hotmail.com> wrote in message
news:v4hBh.2691$Jl....@newsread3.news.pas.earthlink.net...
>
> program test
> REAL, TARGET:: A(100) = [1:100]
> REAL, POINTER:: B
> integer :: i = 50
>
> B => A(I)
> write (*,*) b ! = 50
> b = b+1
> write (*,*) b ! = 51 therefore the increment was "intelligent"
>
> end program
>
I misled myself, the increment is NOT "intelligent" just straight
arithmetic.
>
Walter Spector
> ...
> Not knowing Fortran pointers as well as C pointers, can you have a
> pointer to an array element or structure element?
Yes, and yes. As long as both of the latter have the "target" attribute.
> As far as I know,
> you can't have a pointer to part of an array, as is common in C.
Sure you can. A F90 pointer can even point to a non-unit stride
"slice" of a matrix:
real, target :: a(100,100)
real, pointer :: aslice(:)
:
aslice => a(i,:) ! Row-wise slice of A
W.
P.S., Bonus question: How can one get 'aslice' to point to the diagonal?
(And yes, it can be done.)
Michael Metcalf
"Walter Spector" <w6ws_xt...@earthlink.net> wrote in message
news:45D5BEFA...@earthlink.net...
> P.S., Bonus question: How can one get 'aslice' to point to the diagonal?
> (And yes, it can be done.)
Like this?
real, target :: a(100,100)
real, pointer :: aslice(:)
:
aslice => f(a)
contains
function f(a)
real, intent(in) :: a(:, :)
real, dimension(:), pointer :: f
allocate(f(ubound(a, 1)))
forall(i = 1:10) f(i) = a(i, i)
end function f
end
Regards,
Mike Metcalf
Walter Spector
>
> "Walter Spector" <w6ws_xt...@earthlink.net> wrote in message
> news:45D5BEFA...@earthlink.net...
> > P.S., Bonus question: How can one get 'aslice' to point to the diagonal?
> > (And yes, it can be done.)
>
> Like this?
No, that points to a *copy* of the diagonal.
W.
Richard Maine
> robert....@sun.com wrote:
>
> > On Feb 15, 12:20 pm, glen herrmannsfeldt wrote:
>
> >>Not knowing Fortran pointers as well as C pointers, can you have a
> >>pointer to an array element or structure element? As far as I know,
> >>you can't have a pointer to part of an array, as is common in C.
>
> > If you can have a pointer to an array, you can have a pointer
> > to an element of that array.
>
> I know Fortran doesn't let you do all the things that C does, I
> believe intentionally. I especially believe that Fortran doesn't
> have anything like C's p++ to point to the next element of an array.
It's true that Fortran doesn't allow such pointer arithmetic, but that's
a separate question. Fortran *DOES*, as Robert mentioned, allow you to
point at an element to an array or a part (slice) of an array. Period.
The fact that you can't subsequently increment that pointer doesn't
change that. Yes, it does things differently than C does, but that again
is a separate question. Let's not say that Fortran doesn't allow a
pointer to an array element just because it isn't exactly the same as a
C pointer.
glen herrmannsfeldt
> Michael Metcalf wrote:
>
>>"Walter Spector" <w6ws_xt...@earthlink.net> wrote in message
>>news:45D5BEFA...@earthlink.net...
>>
>>>P.S., Bonus question: How can one get 'aslice' to point to the diagonal?
>>>(And yes, it can be done.)
>>Like this?
>> function f(a)
>> real, intent(in) :: a(:,
>> real, dimension(:), pointer :: f
>> allocate(f(ubound(a, 1)))
>> forall(i = 1:10) f(i) = a(i, i)
>> end function f
> No, that points to a *copy* of the diagonal.
I just noticed that.
My first thought was RESHAPE, but now I believe that makes
a copy, too. Maybe RESHAPE and TRANSFER?
For C, you make arrays of pointers, not exactly the same, but
close enough much of the time. For PL/I, the DEFINED attribute
does it, at least without pointers.
DCL A(10,10) FLOAT, B(10) DEFINED A(1SUB,1SUB);
-- glen
glen herrmannsfeldt
> glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
(snip)
>>I know Fortran doesn't let you do all the things that C does, I
>>believe intentionally. I especially believe that Fortran doesn't
>>have anything like C's p++ to point to the next element of an array.
> It's true that Fortran doesn't allow such pointer arithmetic, but that's
> a separate question. Fortran *DOES*, as Robert mentioned, allow you to
> point at an element to an array or a part (slice) of an array. Period.
> The fact that you can't subsequently increment that pointer doesn't
> change that. Yes, it does things differently than C does, but that again
> is a separate question. Let's not say that Fortran doesn't allow a
> pointer to an array element just because it isn't exactly the same as a
> C pointer.
As I wrote in the one you replied to,
B => B(LBOUND(B)+1,UBOUND(B))
should be fairly similar to C's B++, but, as far as
I know, you can't do B--. That doesn't seem too
unreasonable to me.
-- glen
Richard Maine
> Walter Spector wrote:
>
> > No, that points to a *copy* of the diagonal.
>
> I just noticed that.
>
> My first thought was RESHAPE, but now I believe that makes
> a copy, too. Maybe RESHAPE and TRANSFER?
No. You are focussing on the details and missing the bigger picture.
RESHAPE and TRANSFER are intrinsic functions, which is to say that they
are functions. That's all you need to know. You don't need to ask what
particular functions they are or look into fine details of their
definitions. A function returns a value. That is the most fundamental
point of the definition of a function. It isn't some obscure detail; it
is the central defining point. The other stuff about functions counts as
detail.
A function does not somehow return its argument(s). It typically returns
a value based on its arguments, but that is a fundamentally different
thing.
One might describe this as saying that the meaning is always that of a
copy. An optimizer might occasionally elide the copy if it concludes
that doing so does not change the meaning. But the meaning is still that
of a copy.
Now a function that returns a pointer would be a different matter. The
pointer itself would indeed be a "copy", but it might point to somethng
that isn't.
There is no point in wading through the list of intrinsic functions to
see which one might do the trick. They won't.
I'm pretty sure I recall the trick that Walter alludes to. It has been
posted here before. I'll avoid completely spoiling the fun by giving the
fiull answer. (This has the secondary benefit of avoiding the
embarassment of messing it up. :-)) But I'll offer a hint in the term
"sequence association".
glen herrmannsfeldt
> glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
(snip)
>>My first thought was RESHAPE, but now I believe that makes
>>a copy, too. Maybe RESHAPE and TRANSFER?
> No. You are focussing on the details and missing the bigger picture.
> RESHAPE and TRANSFER are intrinsic functions, which is to say that they
> are functions. That's all you need to know. You don't need to ask what
> particular functions they are or look into fine details of their
> definitions. A function returns a value. That is the most fundamental
> point of the definition of a function. It isn't some obscure detail; it
> is the central defining point. The other stuff about functions counts as
> detail.
I think I agree with you, thought I am not completely sure in
the case of pointer valued functions and pointer assignment.
(Though I don't want to get into the value of a pointer
discussion again.) At the implementation level, you want an
array descriptor describing the discontiguous array that is
the diagonal of the original array. Consider the pointer
assignment
B => C(I,:)
if I understand pointers, B is now a pointer to the non-contiguous
array of the Ith row of C. It should also be possible to write
a function, F(C,I) such that
B => F(C,I)
results in the same B as the above assignment.
So, can one write
B => DIAG(C)
where DIAG is a pointer valued function given a TARGET square
matrix C? I thought that was what we were trying to do.
-- glen
James Giles
...
> B => DIAG(C)
>
> where DIAG is a pointer valued function given a TARGET square
> matrix C? I thought that was what we were trying to do.
Well, there's another (possibly clearer) way using F2003.
Let C_BASE be a rank-one array with N^2 elements (where
N is the size of the square array you want. Now let C be a
rank-two pointer assigned as follows:
C(1:N, 1:N) => C_BASE
C is now your square array. Finally let C_DIAG be a rank-one
array assigned as follows:
C_DIAG(1:N) => C_BASE(1 : N*N : N+1)
C_DIAG is now the name of the diagonal elements of C.
So, Fortran has now reinvented about a third of the linear subscript
transformation feature that Ken Kennedy recommended in 1982.
I guess that's progress.
Note:
Except for different names, the above is the same as the example
in the F2003 document given in note 7.44 (bottom of page 145).
Their example didn't involve a square array. I guess it still qualifies
as one of the diagonals of the array.
--
J. Giles
"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare
Richard Maine
> Richard Maine wrote:
> > glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
> >>My first thought was RESHAPE, but now I believe that makes
> >>a copy, too. Maybe RESHAPE and TRANSFER?
>
> I think I agree with you, thought I am not completely sure in
> the case of pointer valued functions and pointer assignment.
Thus my comment, which you elided,
>> Now a function that returns a pointer would be a different matter.
>> The pointer itself would indeed be a "copy", but it might point to
>> somethng that isn't.
> At the implementation level,...
You don't need to get into implementation-level issues. Everything
important here is well defined by the language standard, without
reference to the implementation level.
> It should also be possible to write
> a function, F(C,I) such that
>
> B => F(C,I)
Yes. But it better be a function that returns a pointer. You cannot have
a non-pointer expression as a target. And yes, it is possible to write a
function to do the trick. I wouldn't recommend it, as I never recommend
writing functions that return pointers. You can also do the same thing
here with a subroutine.
But you won't find any intrinsic functions that do anything close. For a
start, you won't find many intrinsic functions that return pointer
results. You can forget about all the non-pointer ones, which was my
point. Null returns a pointer, as do some of the new C interop things
(ok; they are not intrinsic functions, but instead are functions in an
intrinsic module. That's a distinction that matters for picky points of
definition, but not here). Null won't help. The C interop ones can, as
James showed. But the trick can also be done without them; see my
previous hint.
Michael Metcalf
common/c/a
real, pointer :: a(:, :)
real, pointer :: aslice(:)
allocate(a(10, 10))
aslice => f()
print *, aslice
contains
function f()
real, pointer :: f(:)
common/c/a
real, target :: a(100)
a = (/ ( i, i = 1, 100) /)
f => a(1:100:11)
Walter Spector
> > No, that points to a *copy* of the diagonal.
> >
> Agreed. Is this then what you meant (following Richard's hint)?
>
Well, placing things in COMMON doesn't make them very dynamic.
But that is certainly one approach.
I was thinking more of the following:
program testdiag
implicit none
real, allocatable, target :: a(:,:)
real, pointer :: adiag(:)
integer :: i, j
integer :: n
write (*,'(a)', advance='no') 'size of matrix? '
read *, n
allocate (a(n,n))
a = 0
print *, 'A before:'
do, j=1, size (a,1)
write (*,'(100f5.2)') (a(j,i), i=1, size (a,2))
end do
adiag => getdiag (a, size (a,1))
adiag = 42.42
print *
print *, 'A after:'
do, j=1, size (a,1)
write (*,'(100f6.2)') (a(j,i), i=1, size (a,2))
end do
contains
function getdiag (t, tn)
! Return a F90 pointer which points to the diagonal of T.
integer, intent(in) :: tn
real, target :: t(tn+1,*)
real, pointer :: getdiag(:)
getdiag => t(1,:tn)
end function
end program
The above works fine, in practice, with a lot of different
compilers. I've used it with Intel, SGI IRIX, g95, gfortran,
and Salford compilers. However, it might be argued that it violates
the letter of the Standard... So one must be careful not to use it
in a situation where the compiler will potentially perform
a copyin/copyout on array A during the procedure call.
When copyin/copyout does get in the way, Cray pointers can come in
handy. Or, I presume, some of the new F2003 intrinsics. Time
will tell.
W.
Richard Maine
> "Walter Spector" <w6ws_xt...@earthlink.net> wrote in message
> news:45D5C9C7...@earthlink.net...
>> "Walter Spector" <w6ws_xt...@earthlink.net> wrote in message
>> news:45D5BEFA...@earthlink.net...
>> P.S., Bonus question: How can one get 'aslice' to point to the
>> diagonal?
>> (And yes, it can be done.)
> Agreed. Is this then what you meant (following Richard's hint)?
[code using common]
That variant hadn't occured to me :-) It would probably work (I didn't
check every detail - just noticed the general idea), but is a bit
restrictive in some ways.
I had in mind the same general principle using argument sequence
association instead of common. But yes, the general principle is the
same; the difference is just in what method you use to set up the
sequence association.
highegg
wrote:
> Michael Metcalf wrote:
>
> > "Walter Spector" <w6ws_xthiso...@earthlink.net> wrote in message
using some more non-standard extensions and dirty tricks, one can even
implement a function
that can take the diagonal of an arbitrary array section (i.e., takes
assumed shape array as an argument). The following program works with
gfortran, g95 (with -Wno-globals), Intel and Pathscale:
module diagm
implicit none
contains
function diag(mat)
real,dimension(:,:),intent(in),target:: mat
real,dimension(:),pointer:: diag
integer:: n,inc
integer,parameter:: sizeof_real = 4
interface
subroutine dirty_pointer_1d(n,x,inc,ptr)
integer,intent(in):: n,inc
real:: x ! an intentional lie ...
real,dimension(:),pointer:: ptr
end subroutine
end interface
n = minval(shape(mat))
if (n == 0) then
nullify(diag)
else if (n == 1) then
diag => mat(1:1,1)
else
inc = (loc(mat(2,2)) - loc(mat(1,1))) / sizeof_real
call dirty_pointer_1d(n,mat(1,1),inc,diag)
end if
end function
end module
! this one has to be external
subroutine dirty_pointer_1d(n,x,inc,ptr)
integer,intent(in):: n,inc
real,target:: x(0:*)
real,dimension(:),pointer:: ptr
ptr => x(0:n*inc:inc)
end subroutine
program testp
use diagm
real,dimension(:,:),allocatable:: a
real,dimension(:),pointer:: dg
allocate(a(15,10))
a = 0
write(*,'(10F4.1)') transpose(a)
write(*,'(40("-"))')
read *
dg => diag(a)
dg = 1.1
write(*,'(10F4.1)') transpose(a)
write(*,'(40("-"))')
read *
dg => diag(a(4:12:2,1:6))
dg = 2.2
write(*,'(10F4.1)') transpose(a)
write(*,'(40("-"))')
read *
dg => diag(a(15:1:-1,:))
dg = 3.3
write(*,'(10F4.1)') transpose(a)
end program
robert....@sun.com
wrote:
> "Walter Spector" <w6ws_xthiso...@earthlink.net> wrote in message
>
> news:45D5C9C7...@earthlink.net...> Michael Metcalf wrote:
>
> >>news:45D5BEFA...@earthlink.net...
> >> > P.S., Bonus question: How can one get 'aslice' to point to the
> >> > diagonal?
> >> > (And yes, it can be done.)
>
> >> Like this?
>
> > No, that points to a *copy* of the diagonal.
>
> Agreed. Is this then what you meant (following Richard's hint)?
>
> common/c/a
> real, pointer :: a(:, :)
> real, pointer :: aslice(:)
>
> allocate(a(10, 10))
> aslice => f()
> print *, aslice
>
> contains
>
> function f()
> real, pointer :: f(:)
> common/c/a
> real, target :: a(100)
> a = (/ ( i, i = 1, 100) /)
> f => a(1:100:11)
> end function f
>
> end
Your post might confuse some people. The pointer aslice
is not pointing to the diagonal of anything. Note that
the variable a in the main program and the a in the
internal function are different variables referencing
different storage. If you made the variable a in the main
program an explicit-shape 10x10 array instead of an array
pointer, you would get the effect you want by storage
association. EQUIVALENCE would also work, also because
of storage association.
Bob Corbett
Ron Shepard
Walter Spector <w6ws_xt...@earthlink.net> wrote:
> When copyin/copyout does get in the way, Cray pointers can come in
> handy. Or, I presume, some of the new F2003 intrinsics. Time
> will tell.
[1D array pointing to the diagonal elements of a 2D array...]
F2003 allows pointer assignment of a multimensional array to a
1-dimensional array, so this problem is solved in a fairly
straightforward way. F90 and F95 do not allow that kind of pointer
assignment.
I remember that this surprised me when I first tried it in f90
because I thought that was the main point of introducing pointers
into the language in the first place. I had in mind something more
like a runtime equivalence using pointer assignment syntax. That is
apparently not what the language designers had in mind. Now the
language seems to have finally caught up (at least when f2003
compilers actually become available).
$.02 -Ron Shepard
Richard Maine
> [1D array pointing to the diagonal elements of a 2D array...]
>
> F2003 allows pointer assignment of a multimensional array to a
> 1-dimensional array, so this problem is solved in a fairly
> straightforward way.
Well, no, because tthe rank mismatch is allowed in only one direction -
and it is the wrong direction for this problem. That's a significant
limitation of the feature as it is in f2003. Yes, you can do it if you
realize ahead of time that you'll need this and thus you make the
original array 1-D. But that has this particular requirement forcing you
to change the original code. You can't just take a random 2-D array and
apply this directly.
The reason, by the way, is that a 2-d array might have strides that
would be non-uniform if you tried to map a 1-d array onto it. The
standard doesn't want to force that kind of capability on array dope
vectors; that would have lots of bad consequences. (It would be very
flexible, but the flexibility would come at a higher cost than anyone
wants to pay). If you start with a 1-D array, you never have such
problems; higher rank ones can map on it it just fine in all cases.
There is a proposal for f2008 which improves on this situation by
allowing you to specify that an array is contiguous. If an array is
contiguous, then mappings of any rank should work. And after all, most
arrays are contiguous, so that's not as big a limitation as the target
having to be rank 1. I think that proposal is still "go". It is one I
favor. The contiguous declaration has other benefits as well; this one
is almost a side benefit, but a nice one.
glen herrmannsfeldt
(snip)
> The reason, by the way, is that a 2-d array might have strides that
> would be non-uniform if you tried to map a 1-d array onto it. The
> standard doesn't want to force that kind of capability on array dope
> vectors; that would have lots of bad consequences. (It would be very
> flexible, but the flexibility would come at a higher cost than anyone
> wants to pay). If you start with a 1-D array, you never have such
> problems; higher rank ones can map on it it just fine in all cases.
I hadn't thought of that one before, so I looked up how PL/I
handles this problem. It seems, at least in the current versions,
that it passes a copy of the array if it is iSUB defined.
(iSUB allows non-uniform mapping even if the base array is
contiguous.) In the usual case where PL/I creates a copy
it does not copy the result back, as in the general case
of an expression.
-- glen
James Giles
...
> The reason, by the way, is that a 2-d array might have strides that
In what cases? A 2-d array with stride has a constant interval
between rows and a constant interval between columns. The
only rank-1 linear reductions of that have a constant stride that's
the sum (or difference) of constant multiples of those two existing
strides. It's still uniform. For example, the diagonal of a square
array (even if that array is not contiguous) is the first element
followed by N-1 additional elements whose stride is the sum of
the column stride and the row stride of the original square array.
This works even if one or both of the strides of the original square
array are negative. It even works, if using some mechanism not
permitted in current Fortran, your square array's dope vector
denotes the transpose of some memory array. (Note that trans-
position is also a linear subscript transformation.) Any linear
subscript transformation can be applied to the "dope vector" that
results any previous such transformation with the same kind of
operations that were applied to get the first "dope vector" in the
first place.
Ken Kennedy's recomendation for linear subscript transformations
not only included a description of the feature from the user's point
of view, but also included the mathematics needed by the vendor
to implement the feature. There's some complexity to be sure. But
it's not all *that* hard.
Now, maybe you're worried about non-linear transformations? If
so, starting with a flat array (as Fortran 2003 now permits) won't
let you do those. The current (2003) feature has no more (or less)
transforming capability than Ken Kennedy's proposal - it's just a lot
less convenient. It also unnecessarily involves pointers, which are
only needed when aliasing is involved. Most slicing of arrays ought
to be done without aliasing.
Ron Shepard
nos...@see.signature (Richard Maine) wrote:
> Ron Shepard <ron-s...@NOSPAM.comcast.net> wrote:
>
> > [1D array pointing to the diagonal elements of a 2D array...]
> >
> > F2003 allows pointer assignment of a multimensional array to a
> > 1-dimensional array, so this problem is solved in a fairly
> > straightforward way.
>
> Well, no, because tthe rank mismatch is allowed in only one direction -
> and it is the wrong direction for this problem.
Yes, I knew that. :-) It would be straightforward at least in those
situations in which the 2D array was already based on a 1D array,
which is all of the cases in which this is needed to replace the f77
situations (i.e. through equivalence or storage association).
Going beyond the f77 replacement situation, the spacing between the
elements of a diagonal of any 2D array in fortran is constant, even
if the 2D array has strange strides between elements of the rows and
columns. Therefore, it should not place any additional burden on
the dope vector of a 1D array to point to those diagonal elements.
I think that the fact that f90, f95, and now f2003 did not allow
this kind of pointer assignment or stride specification in a simple
and straightforward way was a mistake in the specification, a design
flaw. After all, this is a very common requirement, and it was a
common programming technique with f77 code. For example, think of
all the times you use a level-1 BLAS call to operate on a diagonal
(or a codiagonal sequence) of a 2D array. There should have been an
array syntax equivalent to that common practice from the beginning.
At least in f2003 you can work around the problem in a legal way by
allocating the original array as 1D, although, as you point out, it
may require changing old working code in order to achieve this.
Another downside of this approach is that it requires pointer
assignment, which in many cases is unnecessary and suppresses
optimization possibilities whereas the equivalent stride notation
(if it were allowed) would be both simpler and more efficient.
These operations are easy to write with standard mathematical
notation, and a language that is based on formula translation should
have accommodated that common need.
$.02 -Ron Shepard <--IMO, of course
Brooks Moses
> Richard Maine wrote:
>> The reason, by the way, is that a 2-d array might have strides that
>> would be non-uniform if you tried to map a 1-d array onto it. [...]
>
> In what cases? A 2-d array with stride has a constant interval
> between rows and a constant interval between columns. The
> only rank-1 linear reductions of that have a constant stride that's
> the sum (or difference) of constant multiples of those two existing
> strides. It's still uniform. For example, the diagonal of a square
> array (even if that array is not contiguous) is the first element
> followed by N-1 additional elements whose stride is the sum of
> the column stride and the row stride of the original square array.
Suppose we have ALLOCATE(A(4,4)). Consider the array A(2:3,2:3). If we
map a 1-D array onto that, we have the first and second elements
adjacent, then a skip over two elements between the second and third,
and then the third and fourth are adjacent again.
Thus, for some A1 => A(2:3,2:3), the elements of A1 have a non-uniform
stride.
- Brooks
--
The "bmoses-nospam" address is valid; no unmunging needed.
glen herrmannsfeldt
> Richard Maine wrote:
>>The reason, by the way, is that a 2-d array might have strides that
>>would be non-uniform if you tried to map a 1-d array onto it. [...]
> In what cases? A 2-d array with stride has a constant interval
> between rows and a constant interval between columns. The
> only rank-1 linear reductions of that have a constant stride that's
> the sum (or difference) of constant multiples of those two existing
> strides.
I hadn't thought of it before, either. If the source is a subarray
of the original array.
REAL A(10,10)
CALL X(A(3:4,5:6))
... (assume appropriate interface)
SUBROUTINE X(B)
REAL B(:,:)
> It's still uniform. For example, the diagonal of a square
> array (even if that array is not contiguous) is the first element
> followed by N-1 additional elements whose stride is the sum of
> the column stride and the row stride of the original square array.
As a 2D array, it is uniform in each dimension separately,
but not, in general, as a 1D array. The diagonal is a special
case, which conveniently is uniform.
The PL/I case is worse, you can make non-linear transformations.
DCL A(10,10) FLOAT, B(10) DEFINED A(1SUB,1SUB);
is a linear map to the diagonal, but
DCL C(1000) FLOAT, D(31) DEFINED C(1SUB**2);
maps elements of D onto C(1), C(4), C(9), ... C(961).
-- glen
James Giles
...
> Suppose we have ALLOCATE(A(4,4)). Consider the array A(2:3,2:3). If
> we map a 1-D array onto that, we have the first and second elements
> adjacent, then a skip over two elements between the second and third,
> and then the third and fourth are adjacent again.
>
> Thus, for some A1 => A(2:3,2:3), the elements of A1 have a non-uniform
> stride.
That's not a linear transformation. You're attempting a map
in array element order. I didn't say that would work. No doubt
if that's what you want you'll have to make a contiguous copy.
It's similar to another case. Suppose that you've got a rank-1
array of size 60 and you transform that separately to a rank-2
array that's 15×4 and another rank-2 array that's 20×3. They
can't be mapped to each other by any simple linear mapping.
You're relying on array element order.
Still, linear mappings accomplish more of the kinds of transforms
that people really need than Fortran currently allows and more
conveniently. Any that linear transforms don't do, current (and
proposed) Fortran can't do either.
Richard Maine
> That's not a linear transformation. You're attempting a map
> in array element order.
Well, yes. But that's exactly how the remapping in question is defined.
It is defined as being in array element order. It isn't defined as a
linear transformation. Yes, this is exactly the problem. Maybe you had
some other proposed feature in mind, but the feature that is in the
standard is based on array element order. Yes, I guarantee that this is
exactly the reason why it isn't allowed for other than rank-1 targets.
That isn't speculation; I was there and remember.
(I also remember wishing we could allow it the other way and perhaps
proposing that it be allowed in the cases where it works, but that
required more concepts to be added than would fly at the time. The
contiguous attribute stuff proposed for f2008 ought to fill in at least
the most important cases.)
James Giles
> James Giles <james...@worldnet.att.net> wrote:
>
>> That's not a linear transformation. You're attempting a map
>> in array element order.
>
> Well, yes. But that's exactly how the remapping in question is
The remapping in question is the diagonal of an array. I had
no idea any other remapping was "in question". I suspect
I'm not much interested in most non-linear remappings since
a direct appeal to erray element order seems to me to automatically
be better handled by copying.
Ron Shepard
Brooks Moses <bmoses...@cits1.stanford.edu> wrote:
> Suppose we have ALLOCATE(A(4,4)). Consider the array A(2:3,2:3). If we
> map a 1-D array onto that, we have the first and second elements
> adjacent, then a skip over two elements between the second and third,
> and then the third and fourth are adjacent again.
>
> Thus, for some A1 => A(2:3,2:3), the elements of A1 have a non-uniform
> stride.
No this isn't correct. Consider the pointer assignment or the
argument association with B=>A(2:3,2:3). The spacing between between
elements in a column, B(:,k) is uniform (=1), the spacing between
elements in a row B(k,:) is uniform (=4), and the spacing between
elements in a diagonal or codiagonal sequence is uniform (=5). Of
course, you need more than a 2x2 B(:,:) to see this, but it is true.
Even if the original array A(:,:) is itself some complicated slice
with strange increments, the B(:,:) array still have uniform spacing
between elements in the above cases, and it should be possible to
map those sequences to a 1D fortran array. As I said before, it was
very common in f77 code to take advantage of these relations.
However, F90 (and now F95 and F2003) array stride notation does not
allow them to be exploited, and the F2003 pointer assignment only
allows it to be exploited in a very special situation.
That's not to say that the F2003 generalization is a bad one. It
will be very useful once compilers are generally available
(hopefully within the next couple of years). However, it does not
solve all the problem that could be solved with the current
dope-vector approach, and it does not do everything with array
notation and pointer assignment that could be done with standard f77
code 30 years ago using storage association tricks.
One area where this is significant is with tensor representations
and transformations of those representations from one basis to
another. These can all be done with minimal data movement using
subscript indexing and f77 storage association. But if you try it
with array syntax notation, it requires a large amount of copying
and memory transfers to accommodate all of the TRANSPOSE() and
RESHAPE() operations that are required. Here is a typical example
from one of my codes:
do i = nmat, 1, -1
! The transpose in the following statement is equivalent
! to a cyclic permutation of the indices.
w(:) = reshape( transpose( matmul( a_pt(i)%a, &
& reshape( w, shape=ishape2 ) ) ), shape=ishape1 )
enddo ! after nmat cycles, the indices
! are returned back to the original order.
w(:) holds the tensor and a_pt(:)%a are the transformation arrays,
one for each rank. All of those transpose and reshape operations
involve making copies of the original data. Those copies are
unnecessary in the f77 version of this same code, and they would be
unnecessary in the f90 version of the code if there were a more
flexible array syntax.
Actually, they would be unnecessary if there were some kind of
pointer assignment that did the same thing as a reshape() and a
transpose(), but operated only on the dope vectors. My mental
picture of this is like a runtime equivalence statement where the
data stays in place, but the strides, dimensions, and bounds are
specified with runtime expressions. I think the f2003 pointer
assignment will handle the reshape(), but I don't think it handles
the transpose().
The above code is about 10 years old, and this was one of my first
programs that used extensive f90 array syntax throughout. I was
very dissapointed when I realized how limiting the array notation
was compared to the old f77 style of coding. Especially considering
that tensor manipulations should be one of the strong points for a
language with array notation. The above code is short and concise,
it just isn't very efficient. The f77 code is longer and somewhat
complicated, but it is efficient. I wanted to have the best of both.
$.02 -Ron Shepard
Richard Maine
> Richard Maine wrote:
> > James Giles <james...@worldnet.att.net> wrote:
> >
> >> That's not a linear transformation. You're attempting a map
> >> in array element order.
> >
> > Well, yes. But that's exactly how the remapping in question is
> > defined. [...]
>
> The remapping in question is the diagonal of an array. I had
> no idea any other remapping was "in question".
The remapping I was referring to was the f2003 feature of allowing
rank-n pointers to point to remapped rank-1 targets. In particular, I
was explaining why the standard does not allow the other way around -
rank-1 pointers to rank-n targets. My comment was about the standard -
not about anything else.
Yes, in the particular case of a diagonal, it always "works", but I was
explaining the reason for the lack of the feature in the standard. The
explanation given is the reason for that. The standard does not have a
special case for the diagonal. If someone wanted a reason for not having
that, I would think it fairly obvious - that it is a bit much of a
special case to have a seperate feature for. I don't recall anyone ever
proposing such a thing as a special-case feature.
James Giles
...
> The remapping I was referring to was the f2003 feature of allowing
> rank-n pointers to point to remapped rank-1 targets. In particular, I
> was explaining why the standard does not allow the other way around -
> rank-1 pointers to rank-n targets. My comment was about the standard -
> not about anything else.
>
> Yes, in the particular case of a diagonal, it always "works", but I
> was explaining the reason for the lack of the feature in the
> standard. The explanation given is the reason for that. The standard
> does not have a special case for the diagonal. If someone wanted a
> reason for not having that, I would think it fairly obvious - that it
> is a bit much of a special case to have a seperate feature for. I
> don't recall anyone ever proposing such a thing as a special-case
> feature.
It was in 1982. And it's not a special case. Rather, the capabilities
of Fortran's pointer remappings are special cases of it. It was a
completely general linear remapping of subscripts.
Richard Maine
> Richard Maine wrote:
> ...
> > The remapping I was referring to was the f2003 feature of allowing
> > rank-n pointers to point to remapped rank-1 targets. In particular, I
> > was explaining why the standard does not allow the other way around -
> > rank-1 pointers to rank-n targets. My comment was about the standard -
> > not about anything else.
> >
> > Yes, in the particular case of a diagonal, it always "works", but I
> > was explaining the reason for the lack of the feature in the
> > standard. The explanation given is the reason for that. The standard
> > does not have a special case for the diagonal. If someone wanted a
> > reason for not having that, I would think it fairly obvious - that it
> > is a bit much of a special case to have a seperate feature for. I
> > don't recall anyone ever proposing such a thing as a special-case
> > feature.
>
> It was in 1982. And it's not a special case. Rather, the capabilities
> of Fortran's pointer remappings are special cases of it. It was a
> completely general linear remapping of subscripts.
Oh. You are talking about the IDENTIFY (or some such; I think that was
the name) proposal in early f8x drafts.
Whatever. I will no longer continue this thread. I cited the reason for
the lack of a particular feature in the f2003 standard. That was, in
fact, the reason for that particular decision at that particular time. I
was there when the decision was made (in the late 90's - not in 1982)
and I know for certain that this was the reason.
I am not at the moment interested in discussing some other feature or
some other decision for some other standard. I am also not interested in
discussing the wisdom or lack thereof of the decision. I seem to be
having quite enough trouble in communicating a simple matter of fact
about the single item that I was trying to say. I'm quite willing to
accept that the failure in communication is mine. In any case, it
doesn't seem like a good basis for a broader discussion. Feel free to
continue with other discussions without me should you be so inclined.
James Giles
...
> Oh. You are talking about the IDENTIFY (or some such; I think that was
> the name) proposal in early f8x drafts.
No. I'm not. I've posted it before. You've read it (and responded
about it before). It was Ken Kennedy's alternative to ALIAS/IDENTIFY
and the other feature SET RANGE/RANGE.
It's of little interest now (to Fortran - it's something designers
of other languages should consider before following Fortran's
mistakes). As I said, with the F2003 pointer stuff, and the
F2008 CONTIGUOUS stuff, they've reinvented (in less convenient
form) about a third of the linear subscript transformation feature.
"All" that's missing is the ability to transform from other than
rank-1 and the removal of the necessity to use pointers to do
the job. It'll probably never happen and if it does I see no
way for it to avoid lots of complexity in order to retain backward
compatibility.
Too bad. Among other things it would have handled all of the
stuff Ron Shepard was asking about.
Brooks Moses
> In article <45D77B1B...@cits1.stanford.edu>,
> Brooks Moses <bmoses...@cits1.stanford.edu> wrote:
>> Suppose we have ALLOCATE(A(4,4)). Consider the array A(2:3,2:3). If we
>> map a 1-D array onto that, we have the first and second elements
>> adjacent, then a skip over two elements between the second and third,
>> and then the third and fourth are adjacent again.
>>
>> Thus, for some A1 => A(2:3,2:3), the elements of A1 have a non-uniform
>> stride.
>
> No this isn't correct. Consider the pointer assignment or the
> argument association with B=>A(2:3,2:3). The spacing between between
> elements in a column, B(:,k) is uniform (=1), the spacing between
> elements in a row B(k,:) is uniform (=4), and the spacing between
> elements in a diagonal or codiagonal sequence is uniform (=5).
Yes, it is correct. :) At least, it is for the case I was talking
about, which is mapping a 1-D array onto A(2:3,2:3), such that A1 is
one-dimensional.
I agree that it works fine when both sides have the same rank, which is
what you're talking about.
Brooks Moses
> Ron Shepard wrote:
>> In article <45D77B1B...@cits1.stanford.edu>,
>> Brooks Moses <bmoses...@cits1.stanford.edu> wrote:
>>> Suppose we have ALLOCATE(A(4,4)). Consider the array A(2:3,2:3). If we
>>> map a 1-D array onto that, we have the first and second elements
>>> adjacent, then a skip over two elements between the second and third,
>>> and then the third and fourth are adjacent again.
>>>
>>> Thus, for some A1 => A(2:3,2:3), the elements of A1 have a non-uniform
>>> stride.
>> No this isn't correct. Consider the pointer assignment or the
>> argument association with B=>A(2:3,2:3). The spacing between between
>> elements in a column, B(:,k) is uniform (=1), the spacing between
>> elements in a row B(k,:) is uniform (=4), and the spacing between
>> elements in a diagonal or codiagonal sequence is uniform (=5).
>
> Yes, it is correct. :) At least, it is for the case I was talking
> about, which is mapping a 1-D array onto A(2:3,2:3), such that A1 is
> one-dimensional.
Perhaps I should have been more clear: Such that A1 is one-dimensional
and refers to all of the elements of A(2:3,2:3) in array-element order.
James Giles
It's plain text. As usual, an errors in transcribing the
paper onto the computer are mine. The link is:
http://home.att.net/~jamesgiles/K_Kennedy.txt
Note that he mainly addresses the ALIAS/IDENTIFY feature
of the F8x proposal. All he says about it applies equally to
the current POINTER feature of Fortran. Indeed, the only
difference is that pointers may also do dynamic allocation
which ALIAS/IDENTIFY could not do. (That makes it worse:
as the paper points out, the main flaw is having a single feature
do too much.)
A small number of simple features that work orthogonally
is almost always better than a big complex feature that does
a whole lot.