Home > Archive > Fortran > June 2004 > whole array manipulation
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
whole array manipulation
|
|
|
| I've recently started using whole array manipulations. There is one
problem - I don't dynamically size my arrays, since I usually read in
data from files where I previously don't know how much data is in the
file (I use open(...,end=blah) to set the size -- incidentally, it
would be nice if open(...,end=blah) could take names instead of just
label numbers - same as do loops, and "EXIT outer" etc)
If want a statement to act on an array as following:
v2vr = SUM(vr3d*vr3d*mass)
where all the quantities are arrays, but the arrays with valid data
are somewhat shorter than the full length of the arrays.
I thought I might be able to restrict sum to just the parts of the
arrays I am interested in, as follows, but it turns out this is not
valid syntax (I feel it ought to be - it is in scripting languages
that deal with whole array operations):
v2vr = SUM((vr3d*vr3d*mass)(1:numparts))
Otherwise I have to go through every mention of an array, and put the
(1:numparts) there. This seems to go against the simplifying nature of
whole array operations:
v2vr(1:numparts) = SUM(vr3d(1:numparts)*vr3d(1:numparts)*ma
ss(1:numparts))
Is there a way I can transfer the whole array into an array of the
required length, and then use sum on that, without involving the
manual use of an intermediate dynamically allocated array?
What would you do?
--
TimC -- http://astronomy.swin.edu.au/staff/tconnors/
Happiness is a warm hard disk.
| |
| beliavsky@aol.com 2004-06-25, 7:43 pm |
|
TimC <tconnors@no.astro.spam.swin.accepted.edu.here.au> wrote:
>I've recently started using whole array manipulations. There is one
>problem - I don't dynamically size my arrays, since I usually read in
>data from files where I previously don't know how much data is in the
>file (I use open(...,end=blah) to set the size -- incidentally, it
>would be nice if open(...,end=blah) could take names instead of just
>label numbers - same as do loops, and "EXIT outer" etc)
>
>If want a statement to act on an array as following:
>
> v2vr = SUM(vr3d*vr3d*mass)
>
>where all the quantities are arrays, but the arrays with valid data
>are somewhat shorter than the full length of the arrays.
>
>I thought I might be able to restrict sum to just the parts of the
>arrays I am interested in, as follows, but it turns out this is not
>valid syntax (I feel it ought to be - it is in scripting languages
>that deal with whole array operations):
>
> v2vr = SUM((vr3d*vr3d*mass)(1:numparts))
>
>Otherwise I have to go through every mention of an array, and put the
>(1:numparts) there. This seems to go against the simplifying nature of
>whole array operations:
> v2vr(1:numparts) = SUM(vr3d(1:numparts)*vr3d(1:numparts)*ma
ss(1:numparts))
Why don't you write x**2 instead of x*x. The compiler will generate the same
code.
>
>Is there a way I can transfer the whole array into an array of the
>required length, and then use sum on that, without involving the
>manual use of an intermediate dynamically allocated array?
>
>What would you do?
Btw if your array starts with 1, as is the default in Fortran, x(1:n) is
the same as x(:n).
If you find it a nuisance to write x(1:n) repeatedly in your code, you can
define a POINTER (say xp) with points to x(1:n) and then use xp in your code.
You must declare array x to have the TARGET attribute. I wonder if doing
can reduce the speed a program.
Sometimes I read a data file to see how many lines and columns there are,
ALLOCATE arrays properly, and then read the data into the file. It's not
elegant, but it's easy to code and may not affect program speed if much more
time is spent computing than reading data. Alternatively, you could dimension
a large array of character variables, read the whole data file into this
array, and then ALLOCATE the arrays used in computations.
----== Posted via Newsfeed.Com - Unlimited-Uncensored-Secure Usenet News==----
http://www.newsfeed.com The #1 Newsgroup Service in the World! >100,000 Newsgroups
---= 19 East/West-Coast Specialized Servers - Total Privacy via Encryption =---
| |
| James Van Buskirk 2004-06-25, 7:43 pm |
| "TimC" <tconnors@no.astro.spam.swin.accepted.edu.here.au> wrote in message
news:slrncd340p.m0n.tconnors@tellurium.ssi.swin.edu.au...
> v2vr = SUM((vr3d*vr3d*mass)(1:numparts))
SUM has a MASK optional argument as in:
v2vr = SUM(vr3d**2*mass, MASK = (/(.TRUE.,i=1,numparts), &
(.FALSE.,i=numparts+1,SIZE(vr3d))/)
> Otherwise I have to go through every mention of an array, and put the
> (1:numparts) there.
Oh, I can see now that you would want to define a helper array:
LOGICAL, allocatable :: helper(:)
! ...
ALLOCATE(helper(size(vr3d)))
! ...
helper(:numparts) = .TRUE.
helper(numparts+1:) = .FALSE.
! ...
v2vr = SUM(vr3d**2*mass, MASK = helper)
> Is there a way I can transfer the whole array into an array of the
> required length, and then use sum on that, without involving the
> manual use of an intermediate dynamically allocated array?
You could have a pointer point at the applicable section of the
array.
> What would you do?
I would normally create a data structure to read the data into,
such as a linked list, and only when the data is read would I
allocate the array and copy the data into it.
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
| |
| Richard Maine 2004-06-25, 7:43 pm |
| "James Van Buskirk" <not_valid@comcast.net> writes:
>
> I would normally create a data structure to read the data into,
> such as a linked list, and only when the data is read would I
> allocate the array and copy the data into it.
Me too. In fact, that's not only what I would do, it's what I do do
in comparable situations. Life is just much simpler if the arrays are
sized just right. Then you don't have to worry about unused parts of
the array and generally keep track of the array size and the used size
as two distinct things. You only have the array size...and the
compiler keeps track of that for you. It ends up being a pretty big
savings in code legibility.
--
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
| |
| Eric K. 2004-06-25, 7:43 pm |
| TimC <tconnors@no.astro.spam.swin.accepted.edu.here.au> wrote in message news:<slrncd340p.m0n.tconnors@tellurium.ssi.swin.edu.au>...
>v2vr(1:numparts) = SUM(vr3d(1:numparts)*vr3d(1:numparts)*ma
ss(1:numparts))
>
> Is there a way I can transfer the whole array into an array of the
> required length, and then use sum on that, without involving the
> manual use of an intermediate dynamically allocated array?
>
> What would you do?
I assume that you want is a scalar-valued result like
v2vr = sum(vr3d(1:numparts)*vr3d(1:numparts)*ma
ss(1:numparts))
Unfortunately, Fortran doesn't allow your preferred syntax. If you have
several instances of the above construct, then you might define a little
helper function like this:
v2vr = f(vr3d,mass,numparts)
...rest of your code...
return
contains
real function f(x,y,n)
integer,intent(in):: n
real,intent(in):: x(n),y(n)
f=sum(y*x**2)
return
end function f
end subroutine
Or, if all your arrays have size NUMPARTS, then you might do all of the work
in another subroutine:
... compute NUMPARTS
call work(vr3d,mass,position,numparts)
where you have something like
subroutine work(vr3d,mass,position,n)
integer,intent(in):: n
real,intent(inout):: vr3d(n), mass(n), position(n)
... now all arrays are "exactly" sized
--Eric
| |
|
| On Thu, 17 Jun 2004 at 13:20 GMT, beliavsky@aol.com (aka Bruce)
was almost, but not quite, entirely unlike tea:
> TimC <tconnors@no.astro.spam.swin.accepted.edu.here.au> wrote:
>
> Why don't you write x**2 instead of x*x. The compiler will generate the same
> code.
Does it?
I thought this came up before, and it was decided that if x is
negative, x**2 is not defined (some implementations might get it right
though)
I don't know how this would apply to arrays; I would guess it would
still be undefined...
--
TimC -- http://astronomy.swin.edu.au/staff/tconnors/
"Legacy (adj): an uncomplimentary computer-industry epithet that
means 'it works'." -- Anthony DeBoer
| |
| James Giles 2004-06-25, 7:43 pm |
| TimC wrote:
> On Thu, 17 Jun 2004 at 13:20 GMT, beliavsky@aol.com (aka Bruce)
> was almost, but not quite, entirely unlike tea:
>
> Does it?
>
> I thought this came up before, and it was decided that if x is
> negative, x**2 is not defined (some implementations might get it right
> though)
Nope, in most implementations, x**2 is perfectly well defined (unless
it overflows, or x is NaN). The standard doesn't actually say much about
exponentiation except that ** is the power operator (and the X1**X2 is
the principal value for X1 complex). I suspect that any implementation
that refused to accept X**2 for X negative would be laughed out of existence.
The only operator that's clearly defined by the standard is integer
divide.
The usual rule of thumb is that X**N (where N is any integer) is
(mathematically) the same as X*X*X...<total of N>...*X (though
not necessarily computationally), whereas X**Y (Y is floatingpoint)
is the same as EXP(Y*LOG(X)). The standard doesn't say so though.
An implementation is actually free to do X**2.0 as X*X if it wants.
I guess it's also free to do X**2 as EXP(2*LOG(X)) if it wants as
well - though, as I say, that would earn it a laughable reputation.
--
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
| |
| Ron Shepard 2004-06-25, 7:43 pm |
| In article
<slrn-0.9.7.4-23964-8100-200406181251-tc@hexane.ssi.swin.edu.au>,
TimC <tconnors@no.astro.spam.swin.accepted.edu.here.au> wrote:
> I thought this came up before, and it was decided that if x is
> negative, x**2 is not defined (some implementations might get it right
> though)
If the exponent is an integer 2, then it is alright for X to be
negative. In fact, this is a fairly common situation.
If the exponent is a floating point number, then I think it is
illegal for X to be negative, even if the exponent has an integer
value (2.0 or otherwise). However, I think X can be complex with a
negative real part and a zero imaginary part and the exponent can be
a floating point number of any value. This is from memory, so it
could be wrong, but I think this may be the situation you were
thinking of.
$.02 -Ron Shepard
| |
| Gerry Thomas 2004-06-25, 7:43 pm |
|
"TimC" <tconnors@no.astro.spam.swin.accepted.edu.here.au> wrote in message
news:slrn-0.9.7.4-23964-8100-200406181251-tc@hexane.ssi.swin.edu.au...
[...]
It amazes moi how many Fortran flops believe they can outwit a compiler
which in all probability isn't even written in Fortran.
> --
> TimC -- http://astronomy.swin.edu.au/staff/tconnors/
> "Legacy (adj): an uncomplimentary computer-industry epithet that
> means 'it works'." -- Anthony DeBoer
Love that legacy bit, a big wet one to your pooch, at least. As applied to
Fortran, it works for now, courtesy of Sgt Maine (genuflect, genuflect,
genuflect) and his self-serving superannuated dozen or so disciples, but in
the 10-15 yr plan, just where's the sustainable legacy when you need it?
Not with a Fortran vendor and the dirty baker's dozen will be as dead as
the unresuscitatable standard that they now so vainly try to maintain.
Better to just pull the plug (or gun, in the case of Americans) and be done
with it.
BTW, someone let a black piglet (they called it a hog but it looked more
like a wog to the lads and all) loose at the local Tranzac and it ran round
blithely supping whatever it was offered until it started to spread eagle
on all fours at which point the fuzz arrived due to a complaint from some
community nosey parker, Ms Piggy left the building and was last sighted
orbiting the CN tower, although this psychotic claim hasn't been
independently confirmed by a reliable source (oxymoronic sober diggers are
not automatically disqualified). Sounds allegorically like the destiny of
Fortran.
--
You're Welcome,
Gerry T.
______
"We are all in this together son!" -- Harry Tuttle (Robert DeNero), in
Brazil (Terry Gilliam).
| |
| glen herrmannsfeldt 2004-06-25, 7:43 pm |
| James Giles wrote:
(snip)
> The usual rule of thumb is that X**N (where N is any integer) is
> (mathematically) the same as X*X*X...<total of N>...*X (though
> not necessarily computationally), whereas X**Y (Y is floatingpoint)
> is the same as EXP(Y*LOG(X)). The standard doesn't say so though.
> An implementation is actually free to do X**2.0 as X*X if it wants.
> I guess it's also free to do X**2 as EXP(2*LOG(X)) if it wants as
> well - though, as I say, that would earn it a laughable reputation.
I have complained a number of times on comp.lang.c that C doesn't
have an integer exponentiation operation, as Fortran (and some
other languages) have.
As above, while the actual algorithm isn't defined in the
standard, it is explicitly allowed to apply an integer power
to all numeric types. There is a convenient algorithm for
computing X**N using O(log N) multiplications, Google for
http://www2.toki.or.id/book/AlgDesignManual/BOOK/BOOK2/NODE54.HTM#SECTION02361000000000000000
seems to have an explanation of the algorithm, though normally it
is done with a loop instead of recursion.
Also, for small constant N Fortran compilers are pretty good
about expanding it inline, with no subroutine call.
-- glen
| |
| David Ham 2004-06-25, 7:43 pm |
| On Fri, 18 Jun 2004 00:07:30 -0500
Ron Shepard <ron-shepard@NOSPAM.comcast.net> wrote:
> In article
> <slrn-0.9.7.4-23964-8100-200406181251-tc@hexane.ssi.swin.edu.au>,
> TimC <tconnors@no.astro.spam.swin.accepted.edu.here.au> wrote:
>
>
> If the exponent is an integer 2, then it is alright for X to be
> negative. In fact, this is a fairly common situation.
>
> If the exponent is a floating point number, then I think it is
> illegal for X to be negative, even if the exponent has an integer
> value (2.0 or otherwise).
Do you mean illegal in the sense of prohibited (which would surprise me
in the Fortran standard) or is it simply the case that the standard
doesn't require that the exponent operator be defined?
David
| |
| Tim Prince 2004-06-25, 7:43 pm |
|
"glen herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
news:NuxAc.67376$HG.34113@attbi_s53...
> James Giles wrote:
> (snip)
>
>
>
> I have complained a number of times on comp.lang.c that C doesn't
> have an integer exponentiation operation, as Fortran (and some
> other languages) have.
>
> As above, while the actual algorithm isn't defined in the
> standard, it is explicitly allowed to apply an integer power
> to all numeric types. There is a convenient algorithm for
> computing X**N using O(log N) multiplications, Google for
>
>
http://www2.toki.or.id/book/AlgDesignManual/BOOK/BOOK2/NODE54. HTM#SECTION02361000000000000000
> seems to have an explanation of the algorithm, though normally it
> is done with a loop instead of recursion.
>
> Also, for small constant N Fortran compilers are pretty good
> about expanding it inline, with no subroutine call.
That's where C misses out, unless you have a compiler which in-lines pow()
for small integral powers, which is allowable, but not normally done. C
pow() must have built in a detection of exact integral powers, as pow() is
required to "Do The Right Thing" for negative operands raised to integral
power. That has been known to spill over into Fortran run-time library
implementations of the X**Y operator (e.g. g77), encouraging people to
produce non-portable code.
| |
| Richard Maine 2004-06-25, 7:43 pm |
| TimC <tconnors@no.astro.spam.swin.accepted.edu.here.au> writes:
> I thought this came up before, and it was decided that if x is
> negative, x**2 is not defined
That is not true. If x is negative, then x**2.0 is not defined,
but that is *VERY* different from x**2, which is perfectly fine.
That difference was the whole point.
--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain
| |
| Dick Hendrickson 2004-06-25, 7:43 pm |
|
Richard Maine wrote:
> TimC <tconnors@no.astro.spam.swin.accepted.edu.here.au> writes:
>
>
>
>
> That is not true. If x is negative, then x**2.0 is not defined,
> but that is *VERY* different from x**2, which is perfectly fine.
> That difference was the whole point.
>
Richard,
Are you sure it isn't defined? I was looking through
section 7 and couldn't see anything that said x**2.0
was undefined for negative x. (Of course, I also couldn't
see anything that said it was defined ;} ) Basically, the
standard says "** is the exponentiation operator, period".
The common understanding of x**y is EXP(y*LOG(x)); but
that's not specified in the standard. I'm not sure that
x**2.0 is anymore "undefined" that is x*2.0 or even
the ever popular 2+2. Furthermore, as James said, a
compiler is allowed to transform x**2 into x**2.0, most
of them don't.
x**2.0 is sort of an edge case, but I think it's a poor
implememtation that doesn't special case integral values
for the exponent. In that sense, it's undefined and
unportable. But the standard doesn't single this out
as being especially undefined (unless I didn't read
enough ;{ )
Dick Hendrickson
| |
| Richard Maine 2004-06-25, 7:43 pm |
| Dick Hendrickson <dick.hendrickson@att.net> writes:
> Richard Maine wrote:
[color=darkred]
> Are you sure it isn't defined? I was looking through
> section 7 and couldn't see anything that said x**2.0
> was undefined for negative x....
Yes, I'm sure. (Which doesn't always mean I'm right. :-))
At first, I couldn't find the prohibition either, and was going
to explain how I thought that it almost followed anyway (in that
I suppose a processor's arithmetic could define it, but that this
definition couldn't depend on whether the 2.0 was a constant or
the value of a variable, and I bet that the processors that
handle this do it only for constants).
But I thought I recalled it being more explicit than that and
I eventually found it. In the f3003 DIS, 2nd sentence of the
2nd para of 7.1.8 (on pg 128).
"Raising a negative-valued primary of type real to a real
power is prohibited."
--
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
| |
| Richard Maine 2004-06-25, 7:43 pm |
| David Ham <d.a.ham@citg.tudelft.nl> writes:
> On Fri, 18 Jun 2004 00:07:30 -0500
> Ron Shepard <ron-shepard@NOSPAM.comcast.net> wrote:
>
> Do you mean illegal in the sense of prohibited (which would surprise me
> in the Fortran standard) or is it simply the case that the standard
> doesn't require that the exponent operator be defined?
Be surprised, I guess. Illegal as in prohibited. See the explicit
quote from the standard in another post I made in this thread.
But I don't know why it would surprise you, so I speculate that
you don't appreciate what "illegal as in prohibited" means in the
standard.
The standard's requirements and prohibitions are, unless stated
otherwise (and this one doesn't state otherwise) prohibitions and
requirements on the program, not on the compiler. Thus it is
prohibited for the program to have x**2.0 for negative x; to put it
another way, any program that does that is nonconforming.
However, it is *NOT* illegal for a compiler to accept the code
and make it work as one might hope (and some compilers do). That
would be a perfectly fine extension.
--
Richard Maine
email: my last name at domain
domain: summertriangle dot net
| |
| David Ham 2004-06-25, 7:43 pm |
| On Fri, 18 Jun 2004 18:02:43 -0700
Richard Maine <nospam@see.signature> wrote:
> David Ham <d.a.ham@citg.tudelft.nl> writes:
>
>
>
> Be surprised, I guess. Illegal as in prohibited. See the explicit
> quote from the standard in another post I made in this thread.
>
> But I don't know why it would surprise you, so I speculate that
> you don't appreciate what "illegal as in prohibited" means in the
> standard.
>
> The standard's requirements and prohibitions are, unless stated
> otherwise (and this one doesn't state otherwise) prohibitions and
> requirements on the program, not on the compiler.
That would be my mistake. Despite the fact that I've seen you say it on
multiple occasions, I keep forgetting that most of the rules in the
standard are addressed to programs not compilers. I think I was
therefore misinterpreting the OP as saying it would be illegal for the
compiler to "do the right thing" with x**2.0 (which would surprise me).
In fact, it is illegal for the program to "ask" the compiler to deal
with that situation which is not surprising at all.
Thanks for straightening that out for me,
David
> Thus it is
> prohibited for the program to have x**2.0 for negative x; to put it
> another way, any program that does that is nonconforming.
>
> However, it is *NOT* illegal for a compiler to accept the code
> and make it work as one might hope (and some compilers do). That
> would be a perfectly fine extension.
>
> --
> Richard Maine
> email: my last name at domain
> domain: summertriangle dot net
|
|
|
|
|