Most elegant syntax for inverting a permutation?
Oskar Enoksson
want to replace the content of P with the inverse permutation.
One way is to allocate a new array IP(:) and then
ALLOCATE(IP(N))
DO I=1,N
IP(P(I)) = I
END DO
P(:) = IP(:)
DEALLOCATE(IP)
I was pretty satisfied when realized I could avoid the need for an extra
array using a forall loop as follows:
FORALL(I=1:N)
P(P(I)) = I
END FORALL
But it made me wonder - perhaps there is an even more elegant way that
e.g. even avoids the loop variable I? Maybe using some fancy intrinsic?
Or maybe some other obvious way that I have failed to see?
/Oskar
Dick Hendrickson
P(P) = [(I, I=1,N)]
Dick Hendrickson
Richard Maine
> I was pretty satisfied when realized I could avoid the need for an extra
> array using a forall loop as follows:
>
> FORALL(I=1:N)
> P(P(I)) = I
> END FORALL
Don't be so sure that this aviods the extra array. FORALL is notorious
for often involving temporary arrays in its implementation. There may
well be an extra array here, even though you don't "see" it in the code.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
glen herrmannsfeldt
> I have an integer array P(:) containing a permutation of 1,2,...,N and
> want to replace the content of P with the inverse permutation.
(snip)
> I was pretty satisfied when realized I could avoid the
> need for an extra array using a forall loop as follows:
> FORALL(I=1:N)
> P(P(I)) = I
> END FORALL
> But it made me wonder - perhaps there is an even more elegant way that
> e.g. even avoids the loop variable I? Maybe using some fancy intrinsic?
> Or maybe some other obvious way that I have failed to see?
You don't avoid the extra array with FORALL, you just don't
see it.
It can be done in place with the assumption that all values
are positive and fit in a signed integer variable with a few
scalar variables. It is very unlikely that FORALL will
do something like this.
(It takes one extra bit per value.
Using the sign bit allows one to avoid a temporary array.)
integer x(10)
data x /1,3,2,5,6,7,4,10,8,9/
print *,x
n=10
j=1
k=x(j)
do i=1,n
if(x(k).lt.0) then
do j=1,n
k=x(j)
if(k.gt.0) exit
enddo
endif
l=x(k)
x(k)=-j
j=k
k=l
enddo
do i=1,n
x(i)=-x(i)
enddo
print *,x
stop
end
-- glen
Oskar Enoksson
> Oskar Enoksson <nob...@nowhere.org> wrote:
>
>> I was pretty satisfied when realized I could avoid the need for an extra
>> array using a forall loop as follows:
>>
>> FORALL(I=1:N)
>> P(P(I)) = I
>> END FORALL
>
> Don't be so sure that this aviods the extra array. FORALL is notorious
> for often involving temporary arrays in its implementation. There may
> well be an extra array here, even though you don't "see" it in the code.
Yes that I'm aware of. I'm just looking for an clear, compact syntax
that does the job.
I even guess it's impossible to avoid O(N) extra storage when inverting
a permutation, except for special classes of permutations.
glen herrmannsfeldt
(snip)
> I even guess it's impossible to avoid O(N) extra storage when inverting
> a permutation, except for special classes of permutations.
It requires one bit per element. If all the values are positive
(assuming array indexing from one) the sign bit is available.
-- glen
Richard Harter
That depends upon the price you are willing to pay. You can
invert a permutation with no extra storage provided that you are
willing to use an O(n*log(n)) algorithm.
Richard Harter, c...@tiac.net
http://home.tiac.net/~cri, http://www.varinoma.com
But the rhetoric of holistic harmony can generate into a kind of
dotty, Prince Charles-style mysticism. -- Richard Dawkins
highegg
You can write the forall as one-liner, or use
P(P) = [(i,i=1,N)] ! (or even [1:N] with Intel's extension)
but of course, if you need to do it repeatedly, it's best to do
call invert(P)
paul.rich...@gmail.com
The temporary is a necessary price to pay to get the correct result.
The FORALL and index tricks above have a dependency that is not
resolved. Try:
integer :: p(4), q(4)
p = (/2,4,1,3/)
FORALL(I=1:4)
q(P(I)) = I
END FORALL
p = q
print *, " Correct inverse perm = ", p
p = (/2,4,1,3/)
FORALL(I=1:4)
P(P(I)) = I
END FORALL
PRINT *, " Using FORALL = ", P
p = (/2,4,1,3/)
p(p) = (/(I, I=1,N)/)
print *, " Using vector index = ", p
END
I get:
Correct inverse perm = 3 1
4 2
Using FORALL = 3 1
4 3
Using vector index = 3 1
4 3
"You don't get owt for nowt", I'm afraid.
Paul
paul.rich...@gmail.com
> p(p) = (/(I, I=1,N)/)
I should either have retained the parameter or replaced 'N' with 4.
Paul
Charles Russell
perhaps there is an even more elegant way
What does "elegant" mean: clear, brief, tricky?
C programmers seem especially fond of idioms that are readily
comprehended only by an expert.
David Frank
"Oskar Enoksson" <nob...@nowhere.org> wrote in message
news:SyPNi.10566$ZA....@newsb.telia.net...
Maybe I dont understand the question, but isnt below the obvious solution?
integer,parameter :: N = 4
integer :: p(N) = [22,44,11,-33]
p = p(N:1:-1)
Dick Hendrickson
> Dear All,
>
> The temporary is a necessary price to pay to get the correct result.
> The FORALL and index tricks above have a dependency that is not
> resolved. Try:
I think you are wrong. There are no "unresolved" dependencies.
Indeed, from the standard's point of view, there are no dependencies.
The form of an assignment statement is
variable = expr
In 7.4.4.2.3, execution of FORALL, it says
"Execution of a forall-assignment-stmt that is an assignment-stmt causes
the evaluation of expr and all expressions within variable for all
active combinations of index-name values. These evaluations may be
done in any order. After all these evaluations have been performed, each
expr value is assigned to the corresponding variable. The assignments
may occur in any order."
In 7.4.1.3 interpretation of intrinsic assignment, it says
"Execution of an intrinsic assignment causes, in effect, the evaluation
of the expression expr and all expressions within variable (7.1.8), the
possible conversion of expr to the type and type parameters
of variable (Table 7.9), and the definition of variable with the
resulting value. The execution of the assignment shall have the same
effect as if the evaluation of all operations in expr and variable
occurred before any portion of variable is defined by the assignment."
In both cases, the compiler is clearly required to completely evaluate
both the expression and any subscripts in the variable before doing any
assignments. This often requires temporaries.
What compiler did you use? Will you submit a bug report?
Dick Hendrickson
In both cases,
Oskar Enoksson
No that's something else. I wanted an inverse permutation PP i.e. such
that PP(P(I)) = I for all I=1..N
Oskar Enoksson
Something like that is what I wanted yes, thanks. I just forgot to
mention that I must use Fortran 95 and the above is Fortran 2003, right?
Oskar Enoksson
4.2.1 and latest gfortran (trunk) gives wrong result. A bug in gfortran
it seems.
/Oskar
Oskar Enoksson
... but "P(P) = (/(I, I=1,N)/)" is Fortran 95 and nicer than the forall
loop imo. Thanks!
/Oskar
Oskar Enoksson
/Oskar
paul.rich...@gmail.com
> I think you are wrong. There are no "unresolved" dependencies.
> Indeed, from the standard's point of view, there are no dependencies.
You are quite right but there is a dependency, I think; the vector
index has to be stored before the assignments does it not? Thus a
temporary is necessary.
> > "You don't get owt for nowt", I'm afraid.
..well, I was right about that bit anyway!
Oskar raised the necessary bug report.
Cheers
Paul
paul.rich...@gmail.com
>
> Oskar raised the necessary bug report.
I have also raised a separate PR33749 for the assignment expression.
Thanks for pointing this out.
Paul