views:

482

answers:

3

I'm going through a Fortran code, and one bit has me a little puzzled.

There is a subroutine, say

SUBROUTINE SSUB(X,...)
REAL*8 X(0:N1,1:N2,0:N3-1),...
...
RETURN 
END

Which is called in another subroutine by:

CALL SSUB(W(0,1,0,1),...)

where W is a 'working array'. It appears that a specific value from W is passed to the X, however, X is dimensioned as an array. What's going on?

+2  A: 

in fortran variables are passed by address. So W(0,1,0,1) is value and address. so basically you pass subarray starting at W(0,1,0,1).

aaa
Thank you for the feedback. I suspected this, but, could you explain how the subarray is shaped? That is where it starts, and how it maps to the new sizing. Does it read/map sequentially into the new sizing?
ccook
+4  A: 

This is non-uncommon idiom for getting the subroutine to work on a (rectangular in N-dimensions) subset of the original array.

All parameters in Fortran (at least before Fortran 90) are passed by reference, so the actual argument is resolved as a location in memory. A location inside the space for the, and the subroutine only manipulates part of the array.

Biggest issue: you have to be aware of how the array is laid out in memory and how Fortran's array indexing scheme works. Fortran uses column major array ordering which is the opposite convention from c. Consider an array that is 5x5 in size (and index both directions from 0 to make the comparison with c easier). In both languages 0,0 is the first element in memory. In c the next element in memory is [0][1] but in Fortran it is (1,0). This affects which indexes you drop when choosing a subspace: if the original array is A(i,j,k,l), and the subroutine works on a three dimensional subspace (as in your example), in c it works on Aprime[i=constant][j][k][l], but in Fortran in works on Aprime(i,j,k,l=constant).

The other risk is wrap around. The dimensions of the (sub)array in the subroutine have to match those in the calling routine, or strange, strange things will happen (think about it). So if A is declared of size (0:4,0:5,0:6,0:7), and we call with element A(0,1,0,1), the receiving routine is free to start the index of each dimension where ever it likes, but must make the sizes (4,5,6) or else; but that means that the last element in the j direction actually wraps around! The thing to do about this is not use the last element. Making sure that that happens is the programmers job, and is a pain in the butt. Take care. Lots of care.

dmckee
Thank you, the example and link are particularly useful. It makes good sense for passing say 0,0 or 0,5, but in this case where it is 0,1,0,1, how does that second 1 come into play? Would the next elements be 1,1,0,1 and then after the first dimension 0,2,0,1?
ccook
Your comment crossed with my edit. The last index given in the calling routine is effectively constant in the called routine.
dmckee
Ah, sorry about that. The edits certainly help, ty again. I suppose I am just having trouble with that second 1 now. Bringing that to 2d, is that to say start from the index of (0,1) the first column, second row, and sequentially 'map' into the array?
ccook
Okay, I'm a fool. The second dimension starts from one and (0,1,0,1) is (0,0,0,const). Thanks!!!
ccook
Actually, even if it does not apply this time, that question *does* come up.
dmckee
My computational methods professor would be *so* proud of me. He was the very master of this nonsense. Thankfully he only made us prove we could wrap out little heads around it once.
dmckee
Thanks again :) I love the tangents that a code can send us on.
ccook
+2  A: 

This is called "sequence association". In this case, what appears to be a scaler, an element of an array (actual argument in caller) is associated with an array (implicitly the first element), the dummy argument in the subroutine . Thereafter the elements of the arrays are associated by storage order, known as "sequence". This was done in Fortran 77 and earlier for various reasons, here apparently for a workspace array -- perhaps the programmer was doing their own memory management. This is retained in Fortran >=90 for backwards compatibility, but IMO, doesn't belong in new code.

M. S. B.
What is the modern approach?
ccook