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

[Robert Corbett]

On 12/09/11 10:59, Toon Moene wrote:
> 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 ?

The intrinsic RESHAPE does not do what you want.  The code you
presented can be simplified.  Try

       A = B( : SIZE(A, 1) - LBOUND(B, 1) + 1, &
              : SIZE(A, 2) - LBOUND(B, 2) + 1, &
              : SIZE(A, 3) - LBOUND(B, 3) + 1, &
              : SIZE(A, 4) - LBOUND(B, 4) + 1)

OSS Fortran will have an easier time generating code for the
revised version.  I suspect other compilers will as well.

Robert Corbett