Fortran - Array with unknown number of dimensions

Author

Array with unknown number of dimensions

Dave Taylor

2005-04-18, 8:57 am

Hi there; I'm struggling with a problem in which F95 seems to demand
more information than I have available before runtime.

The situation is this - I am writing a routine that uses pre-existing
library routines (written in C) to extract data from our archive of
experimental and analysed data. This data may have any number of
dimensions (although in practice it almost always will be in the range
0-3). I am returning, in the function result, the relevant info as a
structure, and I want to return the data as a field of this. The
problem is that I don't know how many dimensions to declare this field
with until the function is run and it has performed tests on the data.
I have two possible bodge solutions, neither of which have universal
support...

i) Return it as a 1-D array, regardless of its true dimensionality. The
user can then reshape it.
ii) Return it as an array of the maximum dimensionality permitted by
Fortran (7?), with trivial trailing dimensions (e.g. Return a 6x6x6 3-D
array as a 7-D array of size (6 6 6 1 1 1 1) ).

Ideally, though, I'd like to return it as an array with the correct
number of dimensions. Is this possible?

Many thanks,
Dave Taylor,
UKAEA Fusion

Michael Metcalf

2005-04-18, 8:57 am

"Dave Taylor" <dtaylor@ukaea.org.uk> wrote in message
news:1113813835.778586.120570@l41g2000cwc.googlegroups.com...
> i) Return it as a 1-D array, regardless of its true dimensionality. The
> user can then reshape it.
> ii) Return it as an array of the maximum dimensionality permitted by
> Fortran (7?), with trivial trailing dimensions (e.g. Return a 6x6x6 3-D
> array as a 7-D array of size (6 6 6 1 1 1 1) ).
>
Both these are certainly valid (and 7 is, indeed, the maximum allowed rank).
In addition, I cobbled together the snippet below, where you have a
structure that defines all possible ranks but activates just one of them.
Just which one can be tested, either by the value of 'rank' or by which
pointer is associated. In any case, dynamic rank as such is not available.

Regards,

Mike Metcalf

module arrays
type t
integer :: rank
integer, pointer :: scalar => null()
integer, pointer, dimension(:) :: vector => null()
integer, pointer, dimension(:,:) :: rank2 => null()
integer, pointer, dimension(:,:,:) :: rank3 => null()
end type t
end module arrays

use arrays
type(t) j
j = myfunc()
if(j%rank == 2) print *, shape(j%rank2), size(j%rank2), j%rank2

contains
type(t) function myfunc()
myfunc%rank = 2
allocate(myfunc%rank2(2, 3))
myfunc%rank2 = 6
end function myfunc

END Program

Dave Taylor

2005-04-18, 3:58 pm

Thanks for the reply; it's not every day a question is answered
personally by the author of your bible on the topic...

Michael Metcalf wrote:
> Both these are certainly valid (and 7 is, indeed, the maximum allowed
rank).
> In addition, I cobbled together the snippet below, where you have a
> structure that defines all possible ranks but activates just one of
them.
> Just which one can be tested, either by the value of 'rank' or by
which
> pointer is associated.

In the interests of trying to prevent the users from using their brains
more than absolutely necessary, I think the best solution from my
perspective is to declare it with the maximum possible number of
dimensions. This should be transparent to the users, shouldn't it,
barring calls to SHAPE? I suspect a call to MATMUL using one of these
would not compile, either. Hmm; I shall s

a few opinions around
here.

> In any case, dynamic rank as such is not available.

is the bottom line here... So my "ideal" solution is not possible. Are
there future plans to modify this in Fortran, or would this be a bad
idea? I'm still getting to grips with the F2003 standard, but I don't
recall allocatable rank being discussed.

Cheers,
Dave

Michael Metcalf

2005-04-18, 3:58 pm

"Dave Taylor" <dtaylor@ukaea.org.uk> wrote in message
news:1113837291.554790.243440@o13g2000cwo.googlegroups.com...
> I suspect a call to MATMUL using one of these
> would not compile, either.
>

Depends. Remember that a(:, :, 1, 1, 1, 1, 1) is a rank-2 array.

>
> is the bottom line here... So my "ideal" solution is not possible. Are
> there future plans to modify this in Fortran, ...

Not that I'm aware of.

Regards,

Mike Metcalf

Richard E Maine

2005-04-18, 3:58 pm

In article <1113813835.778586.120570@l41g2000cwc.googlegroups.com>,
"Dave Taylor" <dtaylor@ukaea.org.uk> wrote:

> Ideally, though, I'd like to return it as an array with the correct
> number of dimensions. Is this possible?

Mike gave you the basic answer (namely, "no"), but I'd like to elaborate
on one or two things.

1. There is one context in which a function can return arrays of varying
rank. That's an elemental function. However, that doesn't seem to fit
your application, because an elemental function does get the actual rank
of a particular invocation determined at compile time. Elementalness
allows you to use a single piece of source code to support multiple
ranks, but it still does require that the determination be done at
compile time.

2. Fortran requires that the rank of *EVERYTHING* be determined at
compile time. Changing this would be a pretty big deal - not impossible,
but a big deal.

3. Related to point 2, exactly how would you think that you could do
anything useful with an array of unknown rank? You asked just about
whether you could create one, but I want to know what you'd do with it.
You couldn't reference an element or slice of it, because the syntax of
referencing an element or slice requires you to know the rank in order
to write down the source code. About the only kinds of things you would
be able to do would be whole array elemental operations... which would
give you results with the same problem.

Of course, you could pass the array to C code. Or you could use sequence
association to pass it to a Fortran subroutine written for a fixed rank
(such as 1). These are basically the same thing. But if all the uses are
going to end up being places where the rank is fixed, then the variable
rank didn't do you any good in the first place - might as well use one
of the two approaches that you mentioned (namely, fix the rank either at
1 or 7 - those being two "obvious" choices, or maybe at the max rank
expected for your application if you can define such a thing).

--
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

glen herrmannsfeldt

2005-04-18, 8:58 pm

Dave Taylor wrote:
> Hi there; I'm struggling with a problem in which F95 seems to demand
> more information than I have available before runtime.

> The situation is this - I am writing a routine that uses pre-existing
> library routines (written in C) to extract data from our archive of
> experimental and analysed data. This data may have any number of
> dimensions (although in practice it almost always will be in the range
> 0-3). I am returning, in the function result, the relevant info as a
> structure, and I want to return the data as a field of this. The
> problem is that I don't know how many dimensions to declare this field
> with until the function is run and it has performed tests on the data.
> I have two possible bodge solutions, neither of which have universal
> support...

> i) Return it as a 1-D array, regardless of its true dimensionality. The
> user can then reshape it.

> ii) Return it as an array of the maximum dimensionality permitted by
> Fortran (7?), with trivial trailing dimensions (e.g. Return a 6x6x6 3-D
> array as a 7-D array of size (6 6 6 1 1 1 1) ).

> Ideally, though, I'd like to return it as an array with the correct
> number of dimensions. Is this possible?
>

What does the C program expect? If you use Fortran assumed-size arrays
the only requirement is that the total number of elements in the caller
not be less than the callee uses.

You can call a Fortran subroutine with a 2D array at one point, and
later in the same program call it with a 3D array, within the standard.
The called routine then must know where the elements are. (Usually
using a 1D array and calculating the position of the element using the
method required by the standard.)

For C programs doing dynamic allocation of 2D arrays, it is often done
using an array of pointers to arrays. In C terms, (int**). Otherwise,
at least before C99, the dimensions must be compile time constants.
This form is, as far as I know, incompatible with Fortran arrays.

-- glen

Dave Taylor

2005-04-21, 8:58 am

Richard E Maine wrote:
> 3. Related to point 2, exactly how would you think that you could do
> anything useful with an array of unknown rank? You asked just about
> whether you could create one, but I want to know what you'd do with
it.

The generic data-reading function I am writing is intended for the use
of many users. These users will know what rank of data they are
expecting (*), so they can declare their local arrays appropriately.
The function, on the other hand, does not know what rank of data to
expect until it has tested the data item in question. This is simply an
effort to remain as generic as possible, in the spirit of the
underlying routines and the equivalent available wrappers in other
languages. There are various possible solutions involving either asking
the user to specify more information or to make two function calls for
each data-reading, but these are rather contrary to the aforementioned
philosophy (spirit) of this.

Cheers,
Dave

(*) or at least they will know once they have examined it, having read
it in using this routine.

Duane Bozarth

2005-04-21, 3:59 pm

Dave Taylor wrote:
>
> Richard E Maine wrote:
> it.
>
....
> ...There are various possible solutions involving either asking
> the user to specify more information or to make two function calls for
> each data-reading, but these are rather contrary to the aforementioned
> philosophy (spirit) of this.

Although I also fail to see the "how" of using it, since you don't seem
to have a problem there I'll skip that to simply comment that while two
function calls to determine the size may see "anti-philosophical", it
would seem that from a users' standpoint that could/would be a hidden
detail...

Jan Vorbrüggen

2005-04-21, 3:59 pm

If you are just doing I/O, you might as well use a buffer, i.e., a
one-dimensional array. You can always pass such a buffer to a routine
expecting an array of any rank, provided the dummy argument has the
same size (or smaller) than the actual. The routine being called,
however, must be declaring the dummy argument with the "correct" rank.

Jan

glen herrmannsfeldt

2005-04-21, 3:59 pm

Dave Taylor wrote:

> Richard E Maine wrote:
>
[color=darkred]
> The generic data-reading function I am writing is intended for the use
> of many users. These users will know what rank of data they are
> expecting (*), so they can declare their local arrays appropriately.
> The function, on the other hand, does not know what rank of data to
> expect until it has tested the data item in question. This is simply an
> effort to remain as generic as possible, in the spirit of the
> underlying routines and the equivalent available wrappers in other
> languages. There are various possible solutions involving either asking
> the user to specify more information or to make two function calls for
> each data-reading, but these are rather contrary to the aforementioned
> philosophy (spirit) of this.

This seems to me what assumed size arrays were designed to do.

You can write something like

SUBROUTINE ADD(A,B,C,N)
REAL A(N),B(N),C(N)
DO 1 I=1,N
1 A(I)=B(I)+C(I)
RETURN
END

This will now add the first N elements of any rank
(assumed size) array passed to it. This was not at all unusual
in the Fortran 66 days, though most would dimension the arrays (1)
instead of (N) as the compilers didn't check them.

Note though that it can't do one common case, using only
part of a multidimensional array.

REAL A(100,100),B(100,100),C(100,100)

(fill 10 by 10 subarrays of B and C)
DO 1 J=1,10
DO 1 I=1,10
B(I,J)=I+J
C(I,J)=I-J
1 CONTINUE

CALL ADD(A,B,C,100)

won't add the right elements.

For this case, more complicated routines were written such as

SUBROUTINE ADD(A,B,C,L,M,N)
REAL A(1),B(1),C(1)
DO 1 I=1,L
DO 1 J=1,M
K=(I-1)*N+J
1 A(K)=B(K)+C(K)
RETURN
END

With this one you can add the subarray with

CALL ADD(A,B,C,10,10,100)

If you find some of the matrix libraries from the Fortran 66 days
all these tricks are there. Sometimes a subroutine was used to do
the subscript calculation so that the complicated code only needed
to be done once.

For the arbitrary rank case, one could be written that would take
arrays for the actual dimensions and desired element and would
compute the offset, all legal Fortran even in newer standards.

-- glen

Richard E Maine

2005-04-21, 8:57 pm

In article <1114070703.213451.295490@l41g2000cwc.googlegroups.com>,
"Dave Taylor" <dtaylor@ukaea.org.uk> wrote:

> Richard E Maine wrote:
> it.
>
> The generic data-reading function I am writing is intended for the use
> of many users. These users will know what rank of data they are
> expecting (*), so they can declare their local arrays appropriately.
> The function, on the other hand, does not know what rank of data to
> expect until it has tested the data item in question.

Then this sounds like a perfectly fine application for one of the
solutions you originally mentioned - declaring the array to be or rank 1.

I suspect that you aren't aware of how array passing worked in f77 and
earlier (and still in the f77-compatible cases of f90). In particular,
it is the feature known these days as sequence association. If it is
only in the user routines that the actual rank is used, then I think
this is the feature for you.

In your routine, just declare the array to be or rank 1, doing any index
arithmetic that you might need. When you pas the array as an actual
argument to the user routine, the ranks do *NOT* have to match as long
as the user routine declares the dummy array to be either explicit shape
or assumed size (i.e. all the possibilities that existed in f77; just
avoid assumed shape). All that has to match is the total number of
elements and the layout in memory.

This works fine as long as you avoid assumed-shape dummies; they do
require that the rank match (otherwise the shape information wouldn't
"fit"). Oh yes, and it also doesn't work if the subroutine is generic
(as in the technical Fortran term - not the broader English usage),
because generic resolution depends on rank, or elemental; but you won't
be in those cases.

--
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

Dave Taylor

2005-04-22, 3:59 pm

Richard E Maine wrote:
> In your routine, just declare the array to be or rank 1, doing any
index
> arithmetic that you might need. When you pas the array as an actual
> argument to the user routine, the ranks do *NOT* have to match as
long
> as the user routine declares the dummy array to be either explicit
shape
> or assumed size (i.e. all the possibilities that existed in f77; just

> avoid assumed shape). All that has to match is the total number of
> elements and the layout in memory.

Were I passing the array back as an argument (or a function result?), I
can see how this might suffice, but the situation is more complex. As I
wrote above, this array is being returned as a field of a structure.
The derived type definition for this structure is given in the module
which contains my reading function. Thus, users can simply preface
their declarations with "USE DATA_READ", and then manipulate the data
in the returned structure as they wish.

I am loath to return each item of information in the structure as a
separate output argument to the function, as this would require the
user to declare a large number of extra variables, many of which they
would not be interested in.

Richard E Maine

2005-04-22, 3:59 pm

In article <1114180468.591516.89300@o13g2000cwo.googlegroups.com>,
"Dave Taylor" <dtaylor@ukaea.org.uk> wrote:

[on the question of arrays with arbitrary rank]

> Richard E Maine wrote:
[color=darkred]
> Were I passing the array back as an argument (or a function result?), I
> can see how this might suffice, but the situation is more complex. As I
> wrote above, this array is being returned as a field of a structure.

Ah. I'm beginning to get the picture. Hmm... Possibilities that occur
to me are.

1&2. The two you mentioned. Declare the array with rank 1 or as the max
rank. The user just has to deal with it. Not very "nice" to the poor
user, I agree.

3. Since the user better know at compile time what rank he expects (even
if your library code is generic and doesn't know), you could have
multiple versions and have the user just call the appropriate one.
Internally, your multiple versions might even call a common core, with
only the top level needing to be different to provide the right
interface.

I was about to comment that you could use a generic procedure
disambiguated based on rank, but then I realized that the rank
difference is inside of the structure, so that doesn't work. You'll need
a different type name for each different rank. (Then you can make the
procedure generic based on that type name, but the user still has to
know and specify the different type names, so the genericity doesn't
help much).

4. You could have multiple arrays in the same structure - one for each
rank. Only "fill" the appropriate ones. To avoid wasting space, make
them pointer or (assuming your compiler supports the allocatable TR),
allocatable. Actually, you probably are making them pointer or
allocatable anyway; in fact, you probably said, but I don't have the
earlier posts handily in front of me to look.

5. If you had f2003, I might be tempted to suggest using inheritance.
Inheritance allows you to take a base type and essentially "add
components". There is (quite a bit) more to it than that and I'm not
sure it would be the best solution here, but it isn't completely
implausible. This does have a little bit the ring of the kind of things
that inheritance can help with. You can't really do it in Fortran yet
anyway until f2003 compilers get out. (The current version of the NAG
compiler does support f2003 inheritance, but I doubt you want to be tied
to a specific compiler).

6. Speaking of f2003, f2003 allows (in multiple ways actually - one
directly via pointer assignment, and one indirectly via C pointer
interop features) a target of rank 1 to be pointed to by a pointer of
higher rank. This could make the "always pass it as rank 1" approach
less painful for the user.

--
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

James Van Buskirk

2005-04-23, 3:57 am

"Richard E Maine" <nospam@see.signature> wrote in message
news:nospam-CFF323.08193422042005@news.supernews.com...

> 6. Speaking of f2003, f2003 allows (in multiple ways actually - one
> directly via pointer assignment, and one indirectly via C pointer
> interop features) a target of rank 1 to be pointed to by a pointer of
> higher rank. This could make the "always pass it as rank 1" approach
> less painful for the user.

Of course, you can try to do this in f95 as well.

http://groups-beta.google.com/group...e69049fe8d1a66e

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

Richard E Maine

2005-04-23, 3:57 am

In article <FfSdnQadudS88fTfRVn-gA@comcast.com>,
"James Van Buskirk" <not_valid@comcast.net> wrote:

> "Richard E Maine" <nospam@see.signature> wrote in message
> news:nospam-CFF323.08193422042005@news.supernews.com...
>
>
> Of course, you can try to do this in f95 as well.
>
> http://groups-beta.google.com/group...m/thread/3e6904
> 9fe8d1a66e

Oh yeah. I'd forgotten about that trick... and the questions about
whether or not it necessarily conforms to the standard. (No, I'm not
prepared to answer that question... and I'm not sure I ever will be -
it's a mess. :-)).

But in f2003 there's a simple one-liner with no complicated guru
contemplation needed. If I recall the syntax correctly,

ptr(lower_1:upper_1,lower_2:upper_2) => target

ought to do the trick for rank2.

--
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

Ron Shepard

2005-04-23, 3:57 am

In article <nospam-AF7A38.14452622042005@news.supernews.com>,
Richard E Maine <nospam@see.signature> wrote:

> But in f2003 there's a simple one-liner with no complicated guru
> contemplation needed. If I recall the syntax correctly,
>
> ptr(lower_1:upper_1,lower_2:upper_2) => target
>
> ought to do the trick for rank2.

There are two things in that statement that make me wonder. Why
wasn't the ability to alias 2D arrays to 1D arrays included in f90?
This is a direct extension of what is done through dummy argument
declarations in f77, so it seemed like an obvious feature to include
in f90. There are some nonstandard (but usually portable)
workarounds in f90 involving subroutines and dummy arguments, but
doing it the straightforward way as above is simpler and clearer.

Second, why wasn't the ability to specify lower and upper bounds in
pointer assignments included in f90? One of the major advantages of
f77 over C has always been the flexibility of specifying lower and
upper bounds in array declarations. But the f90 pointer assignments
removed this flexibility and made that aspect of fortran look like
the less flexible C (but with 1-based indexes rather than 0-based
indexes).

I'm glad that f2003 has added this flexibility, but I wonder why it
wasn't there 15 years ago.

$.02 -Ron Shepard

NuclearWizard

2005-04-23, 3:57 pm

I hope this solution isn't already posted, but if you use a subroutine
instead of a function then you can have multiple data reading
subroutines, for rank 0-3 for instance, and overload the functions to
one generic subroutine,

CALL GET_DATA( data , [arguments])

instead of

data = GET_DATA( [arguments] )

Then use an interface statement in the "preamble" of the module.

INTERFACE GET_DATA
MODULE PROCEDURE GET_DATA_scalar
MODULE PROCEDURE GET_DATA_rank1
MODULE PROCEDURE GET_DATA_rank2
MODULE PROCEDURE GET_DATA_rank3
ENDINTERFACE

Then after the CONTAINS, have all the specific procedures for getting
various sets of data. Of course this method probably requires some
duplicating of code, which is the only downside. I know a function
illuminates the desired use of this procedure much better, but often
times I'm forced not to use functions because you can't overload based
on the result of the function---it's annoying but I understand why it
has to be this way. I like the 7-D array idea, too. It will be easier
to write, the only downside is that there will need to exist a 7-D
dummy array to use at all times. However, this could be avoided by
users if they wrapped an appropriate RESHAPE function around the
GET_DATA function.

data = RESHAPE( GET_DATA([arguments]) , desired_shape )

as opposed to

dummy_data = GET_DATA([arguments])
data = dummy_data( : , : , 1 , 1 , 1 , 1 , 1 ) !for rank-2
example

In summary, I think these are the best 2 ways to solve this problem:
1. return a 7-D array and let the user deal with it (hopefully with the
appropriate call to RESHAPE so we don't need to have a dummy 7-D
array),
2. change the data getting procedure from a function to a subroutine so
it can be overloaded with specific procedures for data of various ranks.

NuclearWizard

2005-04-23, 3:57 pm

I hope this solution isn't already posted, but if you use a subroutine
instead of a function:

CALL GET_DATA( data , [arguments])

instead of

data = GET_DATA( [arguments] )

Then you can have multiple data reading subroutines, for rank 0-3 for
instance, and overload the subroutines to one generic subroutine. Use
an interface statement in the "preamble" of the module which contains
the specific subroutines.

INTERFACE GET_DATA
MODULE PROCEDURE GET_DATA_scalar
MODULE PROCEDURE GET_DATA_rank1
MODULE PROCEDURE GET_DATA_rank2
MODULE PROCEDURE GET_DATA_rank3
ENDINTERFACE

Then after the CONTAINS, have all the specific subroutines for getting
various sets of data. Of course this method probably requires some
duplicating of code, which is the only downside. I know a function
illuminates the desired use of this procedure much better, but often
times I'm forced not to use functions because you can't overload based
on the result of the function---it's annoying but I understand why it
has to be this way. I like the 7-D array idea, too. It will be easier
to write than the overloaded subroutine option, the only downside is
that there will need to exist a 7-D dummy array at all times. However,
this could be avoided by
users if they wrapped an appropriate RESHAPE function around the
GET_DATA function.

data = RESHAPE( GET_DATA([arguments]) , desired_shape )

as opposed to

dummy_data = GET_DATA([arguments])
data = dummy_data( : , : , 1 , 1 , 1 , 1 , 1 ) !for rank-2
example

In summary, I think these are the best 2 ways to solve this problem:
1. return a 7-D array and let the user deal with it (hopefully with the
appropriate call to RESHAPE so we don't need to have a dummy 7-D
array),
2. change the data getting procedure from a function to a subroutine so
it can be overloaded with specific procedures for the data of various
ranks.

Dick Hendrickson

2005-04-23, 3:57 pm

Ron Shepard wrote:
> In article <nospam-AF7A38.14452622042005@news.supernews.com>,
> Richard E Maine <nospam@see.signature> wrote:
>
>
>
>
> There are two things in that statement that make me wonder. Why
> wasn't the ability to alias 2D arrays to 1D arrays included in f90?
Probably no really good reason (other than the usual "that's
the way the vote came out"). There were believed to be
problems when the right hand side was complicated. Just how
do you, in general, map

ptr(lower_1:upper_1,lower_2:upper_2) =>
target(lower:upper:stride)

Or, suppose target was itself the end point in a string
of rank changing remappings? Or was an assumed shape
dummy argument that happened to be discontinuous on the
calling side.

I think the question wasn't "can it be done", but more
along the lines of "Is it important enough to do in
what is an already a big change to F77 when it's not
clear if it can be done well."

> This is a direct extension of what is done through dummy argument
> declarations in f77, so it seemed like an obvious feature to include
> in f90. There are some nonstandard (but usually portable)
> workarounds in f90 involving subroutines and dummy arguments, but
> doing it the straightforward way as above is simpler and clearer.
>
> Second, why wasn't the ability to specify lower and upper bounds in
> pointer assignments included in f90?
I think sort of the same answer. The demand/utility didn't
appear to outweigh the perceived complexity.

Besides, you could specify lower bounds on the right hand side

ptr => target(lower_1:upper_1,lower_2:upper_2)

and (maybe) that would/could be good enough for most uses.

Also, many people viewed pointers as having an assumed shape
dummy argument like quality and dummy arguments couldn't
have their actual lower bounds sensibly imported into a
subroutine.

Those are all just guesses and recollections.

Dick Hendrickson

> One of the major advantages of
> f77 over C has always been the flexibility of specifying lower and
> upper bounds in array declarations. But the f90 pointer assignments
> removed this flexibility and made that aspect of fortran look like
> the less flexible C (but with 1-based indexes rather than 0-based
> indexes).
>
> I'm glad that f2003 has added this flexibility, but I wonder why it
> wasn't there 15 years ago.

>
> $.02 -Ron Shepard

Richard E Maine

2005-04-23, 3:57 pm

In article
<ron-shepard-05E5E8.00251023042005@comcast.dca.giganews.com>,
Ron Shepard <ron-shepard@NOSPAM.comcast.net> wrote:

> Why
> wasn't the ability to alias 2D arrays to 1D arrays included in f90?
> This is a direct extension of what is done through dummy argument
> declarations in f77, so it seemed like an obvious feature to include
> in f90.

Add the usual caveat about this kind of question. I have no direct
knowledge of this at all. The only thing I can contribute to it are
3rd-hand observations just like yours.

The main observation that I would add is that the f77 feature you are
talking about is sequence association. Sequence association pretty much
had to be in f90 because it not only was in f77, but it was very widely
used. However...

Sequence association has lots of complications, limitations, and strange
special cases. It doesn't really "fit" very well with the newer parts of
the language. I rather doubt that there was a lot of enthusiasm about
extending sequence association into new areas.

And I'd guess it would have been a non-trivial extension. Also, you'd
have still gotten questions about why it had such "strange" limitations,
just like I'm quite sure there will be questions about the "strange"
limitations of the new f2003 feature. In particular, the target *MUST*
be rank 1. You can point a rank 2 pointer at a rank 2 target, but then
you can't go the other way. And you can't point a rank 3 pointer at a
rank 2 target, so it isn't just relative rank. It is only for rank 1 as
the target. I'm quite sure that a lot of people are going to find this
limitation strange, but it introduces all kinds of complexity to allow
anything else.

Also note that pointers were a quite late addition to f90.

Oh, and if you want to contemplate an oddity in pointers and bounds, I'd
draw your attention more to the difference between assumed shape and
pointer dummy arguments. This is without any rank funniness. You can
specify the lower bounds of an assumed shape dummy. In fact, you
basically always specify them, explicitly or implicitly. There isn't an
option to assume the lower bounds. Yes, there are good reasons for this
as at least the default. It would be messier to write most subroutines
if this weren't the case, and a lot of people would screw it up.
Basically, you would have to think about lower bounds even if you never
used nondefault ones. But occasionally it is inconvenient. Contrast that
with pointer dummies where you always assume the lower bound and cannot
specify it. Yes, I see how it came about, but it still seems odd to me.

> Second, why wasn't the ability to specify lower and upper bounds in
> pointer assignments included in f90?

Well upper bounds is the same question. Note that you still can't do
that unless the target is of rank 1. If the target isn't rank 1, then
the memory layout might not be "right" unless you throw in a lot of
complicated and arbitrary-sounding restrictions.

For lower bounds, I have no special insight other than the note about
pointers being a late addition (and a pretty complicated one for
something as late as they were), and the usual comment about hindsight
not being very much the same as being there.

--
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

Dave Taylor

2005-04-26, 8:59 am

NuclearWizard wrote:
> you can have multiple data reading subroutines, for rank 0-3 for
> instance, and overload the subroutines to one generic subroutine.

Unfortunately, the rank cannot be inferred from the desired input
arguments without opening up the data item in question and testing it.
If one could have the testing part of the code contained in the
interface block prior to (and informing) the overloading statements,
this would work, but I don't think that this is allowed.

Dave

NuclearWizard

2005-04-26, 4:00 pm

Oh so the USER doesn't even know what rank of data they are expecting?
With the average user knowing so little about the data they will be
getting, it makes me wonder what they will do with it. Maybe just send
it to another routine that actually performs some calculations with the
data? If this is the case, then why have the user be the mailman,
delivering data to the calculational routine? I would say it would be
better to add some optional argument(s) to the calculational routine
which specify what data to use. Then the GET_DATA procedure would be
internal to the calculational routine and hidden from the user.

If I was designing a data-getting procedure to allow users to
potentially manipulate data and use this data in their own routines, I
would write GET_DATA to allow usage like this:

!user declares a variable
REAL(prec),POINTER :: data2(:,:)

!user gets wind speed data
CALL GET_DATA( data2 , 'wind_speeds' , date='01-01-01' )
!data2 is now associated with wind speed data

!user calls some function that needs wind speeds
mesh = GET_MESH( wind_speeds=data2 )

!user calls some subroutine that needs data derived from wind speeds
[mesh]
CALL CALC( mesh , init , soln )

!eventually deallocates wind speed data
DEALLOCATE( data2 )

where for this case the user has looked up the data key, 'wind_speeds',
in some table and sees that wind speed data is a returned as an array
of latitude and longitude (rank 2) so they need to send GET_DATA some
rank 2 structure or their code will get a runtime error and come to a
screeching halt. (I rather have it yield a compile-time error but I
feel like it is impossible in this situation with a data-lookup via
some key.)

The specific subroutine for rank 2 data might look something like this.

SUBROUTINE GET_DATA_rank2( data , key , date )
REAL(prec),INTENT(OUT) :: data(:,:)
CHARACTER(*),INTENT(IN) :: key
CHARACTER(8),INTENT(IN),OPTIONAL :: date
....
SELECT CASE(key)
CASE('wind_speeds')
!access database to get 2D wind speed data
CASE('something_else_2D')
!etc.
CASE DEFAULT
WRITE(*,*)'**RUNTIME ERROR**'
WRITE(*,*)' key='//key//' does not exist or is not rank 2!'
STOP
ENDSELECT
....
ENDSUBROUTINE

Keep in mind these problems with the 7-D array:
1. every calculational routine will need to accept 7-D data for all its
data-type arguments---but then only use the dimensions it needs
2. OR the user DOES know what rank the data is and can send it to a
calculational routine as data(:,:,:,:,1,1,1), for a 4D array for
instance.

I think 1. is prone to error and 2. indicates that the option I
described or Mike Metcalf's option (which would require some fairly
easy manipulation of calculational routines to accept the new data
structure) is viable. Also with 7D arrays required for everything, you
will lose some speed with all the unnecessary extra array index
mappings (7 to 1 vs. the actual number of dimensions to 1.)

glen herrmannsfeldt

2005-04-27, 3:59 pm

Richard E Maine wrote:

(snip)

> Oh, and if you want to contemplate an oddity in pointers and bounds, I'd
> draw your attention more to the difference between assumed shape and
> pointer dummy arguments. This is without any rank funniness. You can
> specify the lower bounds of an assumed shape dummy. In fact, you
> basically always specify them, explicitly or implicitly. There isn't an
> option to assume the lower bounds. Yes, there are good reasons for this
> as at least the default.

> It would be messier to write most subroutines
> if this weren't the case, and a lot of people would screw it up.
> Basically, you would have to think about lower bounds even if you never
> used nondefault ones. But occasionally it is inconvenient. Contrast that
> with pointer dummies where you always assume the lower bound and cannot
> specify it. Yes, I see how it came about, but it still seems odd to me.

For comparison purposes only, that is pretty much how PL/I does it.
Loops over the array should use LBOUND and HBOUND for the limits.
Not having legacy assumed size arrays might have helped, though.
(There is also a DIM function, similar to Fortran's SIZE.)

In the implementations I know, the array descriptor contains the address
of the element with all subscripts zero, even if it isn't within the
array.

[color=darkred]
> Well upper bounds is the same question. Note that you still can't do
> that unless the target is of rank 1. If the target isn't rank 1, then
> the memory layout might not be "right" unless you throw in a lot of
> complicated and arbitrary-sounding restrictions.

> For lower bounds, I have no special insight other than the note about
> pointers being a late addition (and a pretty complicated one for
> something as late as they were), and the usual comment about hindsight
> not being very much the same as being there.

On the other hand, traditional PL/I pointers have no dimension
information at all. This may have changed over the years, though.

-- glen

Sponsored Links