(j3.2006) If SIZE(B) > SIZE(A), what part of B does RESHAPE(B, SHAPE(A)) contain, if both A, B are rank N arrays ?
Toon Moene
toon
Fri Dec 9 13:59:23 EST 2011
It might sound academic, but this question results from a real
(standards' violation) problem I solved yesterday.
Somewhere in the 3.5 million lines of code of our weather forecasting
code, someone had written:
A(:,:,:,:) = B(:,:,:,:)
Looked perfectly OK, has worked for years on various platforms, until I
tried it at home using gfortran 4.6.2.
BWAM - segmentation fault.
Due to the way the code was set up, I knew both arrays were
allocatables, so my first thought was: one of them is not allocated.
The second thought was: B is larger than A. So I added some debug
prints that showed that both were allocated, and B was larger than A
(it's shape differed in one position).
So why did this "work" with the other compilers ? Fortunately, I had
seen *that* one before; the answer is very simple: The other compilers
use the shape of the LHS of the assignment to construct their
scalarization loops, and gfortran uses the shape of the RHS; after all,
if you have to scalarize expressions, in case of
PRINT*, A+B+C
why wouldn't you re-use that code for assignments ?
So now my task was to change this non-standard code (that was working
correctly under the assumption that the shape of the LHS determined the
copying) into something that would work independent of this assumption.
I came up with:
ISHAPE = SHAPE(A) ! ISHAPE is a local 4 element integer array
A(:,:,:,:)=B(LBOUND(B,1):ISHAPE(1)-LBOUND(B,1)+1,&
& LBOUND(B,2):ISHAPE(2)-LBOUND(B,2)+1,&
& LBOUND(B,3):ISHAPE(3)-LBOUND(B,3)+1,&
& LBOUND(B,4):ISHAPE(4)-LBOUND(B,4)+1)
However, today I thought: I'm just doing RESHAPE's job here - by hand.
Certainly, one must be able to replace the above with the following:
A = RESHAPE(B, SHAPE(A))
Unfortunately, perusing paragraph 13.7.140 of the 2008 Standard
(actually the 09-007r2 pdf), I can't convince myself that this is the
case. It says remarkably little about RESHAPEing a larger array to a
smaller: "The elements of the result ... are those of SOURCE in normal
array element order ...".
This doesn't sound like the rank N analogon of the "upper-left hand
corner" of a rank 2 array, which the above hand-coded assignment
certainly is.
Can RESHAPE replace the above convolution ?
Thanks in advance,
--
Toon Moene - e-mail: toon at moene.org - phone: +31 346 214290
Saturnushof 14, 3738 XG Maartensdijk, The Netherlands
At home: http://moene.org/~toon/; weather: http://moene.org/~hirlam/
Progress of GNU Fortran: http://gcc.gnu.org/wiki/GFortran#news
More information about the J3
mailing list