Array indexing question
Paul van Delst
I have some definitions that look like so:
INTEGER, PARAMETER :: N_NLTE_CHANNELS = 3
INTEGER, PARAMETER :: N_CHANNELS = 8
! Define some pretend channel numbers and their coefficient indices
INTEGER, PARAMETER :: SENSOR_CHANNEL(N_CHANNELS) = (/10, 12, 17, 20, 22, 30, 33, 34/)
INTEGER, PARAMETER :: C_INDEX(N_CHANNELS) = (/ 0, 0, 0, 1, 2, 0, 0, 3/)
INTEGER, PARAMETER :: NLTE_CHANNEL(N_NLTE_CHANNELS) = (/20, 22, 34/)
(they're parameters since they are in my unit test code).
As you can see the C_INDEX array is simply a count of those channels in SENSOR_CHANNEL that are listed in the
NLTE_CHANNEL array. The C_INDEX array is used to reference a separate array containing data for the "NLTE-affected"
channels (in this case 20, 22, and 34).
Anyway, I was wanting to build the C_INDEX array dynamically if I was only given the SENSOR_CHANNEL and NLTE_CHANNEL
arrays. Doing it via loops and counters, I can do this:
integer :: i, n
integer :: cidx(N_CHANNELS)
cidx = 0
i=1
do n=1,N_CHANNELS
if ( SENSOR_CHANNEL(n) == NLTE_CHANNEL(i) ) then
cidx(n) = i
i = i + 1
end if
end do
print *, cidx
which gives me
0 0 0 1 2 0 0 3
the same as C_INDEX above.
No worries. But does anyone have a method to do this more simply? E.g. maybe using UNPACK? Or MERGE? Or...?
cheers,
paulv
James Van Buskirk
news:ijmho8$d8b$1...@speranza.aioe.org...
> 0 0 0 1 2 0 0 3
C:\gfortran\clf\Aindex>type Aindex.f90
program Aindex
implicit none
integer,parameter :: SENSOR_CHANNEL(*) = &
[10,12,17,20,22,30,33,34]
integer,parameter :: NLTE_CHANNEL(*) = [20,22,34]
integer,parameter :: N_NLTE_CHANNELS = size(NLTE_CHANNEL)
integer,parameter :: N_CHANNELS = size(SENSOR_CHANNEL)
integer i
integer,parameter :: C_INDEX(*) = unpack( &
vector = [(i,i=1,size(SENSOR_CHANNEL))], &
mask = any(spread(SENSOR_CHANNEL,2,size(NLTE_CHANNEL))== &
spread(NLTE_CHANNEL,1,SIZE(SENSOR_CHANNEL)),2), &
field = 0)
character(20) fmt
write(fmt,'(a,i0,a)') '(a,t19,',size(SENSOR_CHANNEL),'(i3:","))'
write(*,fmt) 'SENSOR_CHANNEL = ',SENSOR_CHANNEL
write(fmt,'(a,i0,a)') '(a,t19,',size(NLTE_CHANNEL),'(i3:","))'
write(*,fmt) 'NLTE_CHANNEL = ',NLTE_CHANNEL
write(*,'(a,t19,i3)') 'N_NLTE_CHANNELS = ',N_NLTE_CHANNELS
write(*,'(a,t19,i3)') 'N_CHANNELS = ',N_CHANNELS
write(fmt,'(a,i0,a)') '(a,t19,',size(C_INDEX),'(i3:","))'
write(*,fmt) 'C_INDEX = ',C_INDEX
end program Aindex
C:\gfortran\clf\Aindex>gfortran Aindex.f90 -oAindex
C:\gfortran\clf\Aindex>Aindex
SENSOR_CHANNEL = 10, 12, 17, 20, 22, 30, 33, 34
NLTE_CHANNEL = 20, 22, 34
N_NLTE_CHANNELS = 3
N_CHANNELS = 8
C_INDEX = 0, 0, 0, 1, 2, 0, 0, 3
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
Gib Bogle
I suppose I'm revealing my ignorance, but to me Paul's code is much
easier to read and understand.
glen herrmannsfeldt
(snip)
>> integer,parameter :: C_INDEX(*) = unpack(&
>> vector = [(i,i=1,size(SENSOR_CHANNEL))],&
>> mask = any(spread(SENSOR_CHANNEL,2,size(NLTE_CHANNEL))==&
>> spread(NLTE_CHANNEL,1,SIZE(SENSOR_CHANNEL)),2),&
>> field = 0)
(snip)
> I suppose I'm revealing my ignorance, but to me Paul's code is much
> easier to read and understand.
That is usual for this type of question.
Note though, in this case, the OP wanted a PARAMETER, which
restricts which functions can be used in its evaluation.
Not that it matters so much in this case, but these solutions
often use more memory and likely run slower. Specfically, for:
spread(SENSOR_CHANNEL,2,size(NLTE_CHANNEL))==
spread(NLTE_CHANNEL,1,SIZE(SENSOR_CHANNEL))
As to the original question, it isn't so obvious that wasting a
few array element by not using the indirect C_INDEX would be so bad.
Though that depends much on the actual problem, assuming that
this is a simplified example.
-- glen
sturlamolden
> I suppose I'm revealing my ignorance, but to me Paul's code is much
> easier to read and understand.
It's better to be explicit than concise. There is no price for winning
the obfuscated array syntax contest, only unmaintainable code.
Sturla
James Van Buskirk
news:534dc89c-8c68-48f8...@n1g2000yqm.googlegroups.com...
Nice to see that I've accreted a couple more cheering fans!
Obviously you are cheering because I was able to make C_INDEX into
a named constant so that it can be protected from alteration by
the program. Some compilers might even be smart enough to put it
in read-only memory so that code like:
call sub(C_INDEX)
...
subroutine sub(x)
integer x(*)
x(5) = 42
end subroutine sub
Could be terminated with extreme prejudice at runtime. Unfortunately
I don't own a compiler that's smart enough to take this action.
Or was the adulation because one could look at the expression
defining C_INDEX and know, without looking through all the INCLUDE
files and modules that make up the program of which this expression
was only an infinitesimal part, that it was the same length as
SENSOR_CHANNEL, even if SENSOR_CHANNEL got its bounds changed from
1:N_CHANNELS during some cycle of code maintainance?
Or that similarly one could tell from the expression in isolation
that it compared with all elements of NLTE_CHANNEL?
Or because I was able to use UNPACK just as the O.P. requested?
Whether for any of these reasons or otherwise, I accept your
spontaneous applause.
<bows>
sturlamolden
> Whether for any of these reasons or otherwise, I accept your
> spontaneous applause.
Here in Scandinavia, bragging about your own achievements is
considered shameful.
I have proven my skills in vectorization before (not just Fortran, but
also Matlab and Python), I don't have to repeat -- particularly when
it's not useful.
Obviously you did not understand the clarity part. Please read the Zen
of Python (yes I know this is comp.lang.fortran), line 2, 3, 7 and 17:
- Explicit is better than implicit.
- Simple is better than complex.
- Readability counts.
- If the implementation is hard to explain, it's a bad idea.
http://www.python.org/dev/peps/pep-0020/
Sturla
glen herrmannsfeldt
> Obviously you did not understand the clarity part. Please read the Zen
> of Python (yes I know this is comp.lang.fortran), line 2, 3, 7 and 17:
> - Explicit is better than implicit.
> - Simple is better than complex.
> - Readability counts.
> - If the implementation is hard to explain, it's a bad idea.
> http://www.python.org/dev/peps/pep-0020/
Much of the time, this is true. It seems to be especially
true when trying to do something that is easy to do with DO
loops, but not so easy with array intrinsic functions.
But there are also cases, such as sort algorithms or the FFT,
where the simple one is O(n**2), and the more complex and
hard to explain one is O(n log n). (Though for small n, the
simpler one is often good enough or faster.)
-- glen
sturlamolden
> Much of the time, this is true. It seems to be especially
> true when trying to do something that is easy to do with DO
> loops, but not so easy with array intrinsic functions.
Many with experience from Matlab or Python falsely get the idea that
DO loops are inherently evil. Compiled DO loops in Fortran are not,
which is an important destinction. Thus, Fortran's array intrinsics
are for writing equations easily, not for the speed of compiled
programs.
Sturla
Paul van Delst
> "sturlamolden" <sturla...@yahoo.no> wrote in message
> news:534dc89c-8c68-48f8...@n1g2000yqm.googlegroups.com...
>
>> On 18 Feb, 22:13, Gib Bogle <g.bo...@auckland.ac.nz> wrote:
>
>>> I suppose I'm revealing my ignorance, but to me Paul's code is much
>>> easier to read and understand.
>
>> It's better to be explicit than concise. There is no price for winning
>> the obfuscated array syntax contest, only unmaintainable code.
>
> Nice to see that I've accreted a couple more cheering fans!
Good one. :o)
> Obviously you are cheering because I was able to make C_INDEX into
> a named constant so that it can be protected from alteration by
> the program. Some compilers might even be smart enough to put it
> in read-only memory so that code like:
>
> call sub(C_INDEX)
> ...
> subroutine sub(x)
> integer x(*)
> x(5) = 42
> end subroutine sub
>
> Could be terminated with extreme prejudice at runtime. Unfortunately
> I don't own a compiler that's smart enough to take this action.
>
> Or was the adulation because one could look at the expression
> defining C_INDEX and know, without looking through all the INCLUDE
> files and modules that make up the program of which this expression
> was only an infinitesimal part, that it was the same length as
> SENSOR_CHANNEL, even if SENSOR_CHANNEL got its bounds changed from
> 1:N_CHANNELS during some cycle of code maintainance?
This is pretty much the crux of the issue. I don't want to have to change source (or test) codes simply because I add or
remove a particular satellite sensor from the list to be used (or tested). I would probably remember what needs to be
changed, but the other 5-10 or so people that also develop the software might not (we don't have a continuous
integration tool set up yet that would highlight these sorts of things).
> Or that similarly one could tell from the expression in isolation
> that it compared with all elements of NLTE_CHANNEL?
>
> Or because I was able to use UNPACK just as the O.P. requested?
>
> Whether for any of these reasons or otherwise, I accept your
> spontaneous applause.
I dunno about anyone else, but:
<claps>
:o)
Thanks for the code.
cheers,
paulv
Paul van Delst
I agree, but....
It's not an either-or matter. It's a continuum[*]. Being too explicit can produce code that may be easy to understand,
but requires constant, even if easy, maintenance. Being too concise can produce code that may be infinitely flexible,
but unintelligible to all but the cognoscenti. Or vice-versa(?)
I like having the ability to consider problems from both sides of the spectrum[*].
cheers,
paulv
[*] given the field I am in, the sort-of puns are intended :o)
James Van Buskirk
> James Van Buskirk wrote:
>> Or was the adulation because one could look at the expression
>> defining C_INDEX and know, without looking through all the INCLUDE
>> files and modules that make up the program of which this expression
>> was only an infinitesimal part, that it was the same length as
>> SENSOR_CHANNEL, even if SENSOR_CHANNEL got its bounds changed from
>> 1:N_CHANNELS during some cycle of code maintainance?
> This is pretty much the crux of the issue. I don't want to have to change
> source (or test) codes simply because I add or
> remove a particular satellite sensor from the list to be used (or tested).
> I would probably remember what needs to be
> changed, but the other 5-10 or so people that also develop the software
> might not (we don't have a continuous
> integration tool set up yet that would highlight these sorts of things).
Even given these constraints we have freedom to choose between
different implementations. Here is a more 'loopy' version, you
decide whether it's more readable:
C:\gfortran\clf\Aindex>type Cindex.f90
program Cindex
implicit none
integer,parameter :: SENSOR_CHANNEL(*) = &
[10,12,17,20,22,30,33,34]
integer,parameter :: NLTE_CHANNEL(*) = [20,22,34]
integer,parameter :: N_NLTE_CHANNELS = size(NLTE_CHANNEL)
integer,parameter :: N_CHANNELS = size(SENSOR_CHANNEL)
integer i
integer,parameter :: C_INDEX(*) = unpack( &
vector = [(i,i=1,size(SENSOR_CHANNEL))], &
mask = [(any(SENSOR_CHANNEL(i) == NLTE_CHANNEL), &
i=lbound(SENSOR_CHANNEL,1),ubound(SENSOR_CHANNEL,1))], &
field = 0)
character(20) fmt
write(fmt,'(a,i0,a)') '(a,t19,',size(SENSOR_CHANNEL),'(i3:","))'
write(*,fmt) 'SENSOR_CHANNEL = ',SENSOR_CHANNEL
write(fmt,'(a,i0,a)') '(a,t19,',size(NLTE_CHANNEL),'(i3:","))'
write(*,fmt) 'NLTE_CHANNEL = ',NLTE_CHANNEL
write(*,'(a,t19,i3)') 'N_NLTE_CHANNELS = ',N_NLTE_CHANNELS
write(*,'(a,t19,i3)') 'N_CHANNELS = ',N_CHANNELS
write(fmt,'(a,i0,a)') '(a,t19,',size(C_INDEX),'(i3:","))'
write(*,fmt) 'C_INDEX = ',C_INDEX
end program Cindex
C:\gfortran\clf\Aindex>gfortran Cindex.f90 -oCindex
f951.exe: internal compiler error: in add_init_expr_to_sym, at
fortran/decl.c:14
18
Please submit a full bug report,
with preprocessed source if appropriate.
See <http://gcc.gnu.org/bugs.html> for instructions.
steve
> sturlamolden wrote:
> > On 18 Feb, 22:13, Gib Bogle <g.bo...@auckland.ac.nz> wrote:
>
> >> I suppose I'm revealing my ignorance, but to me Paul's code is much
> >> easier to read and understand.
>
> > It's better to be explicit than concise. There is no price for winning
> > the obfuscated array syntax contest, only unmaintainable code.
>
> I agree, but....
>
> It's not an either-or matter. It's a continuum[*]. Being too explicit can produce code that may be easy to understand,
> but requires constant, even if easy, maintenance. Being too concise can produce code that may be infinitely flexible,
> but unintelligible to all but the cognoscenti. Or vice-versa(?)
>
I would think that the large comment above the code would
explain what the code is doing. :-)
--
steve