Are Fortran 90 arrrays contiguous in memory?
Carl Christian Kjelgaard Mikkelsen
Suppose I have
subroutine test1(m,n)
integer :: m,n
double precision, dimension(:,:), allocatable :: a
allocate(a(m,n))
!... omitted do-loops to initialize a
a(i,j)=(j-1)*n+i
!...
end subroutine test1
then I am reasonably certain that the elements of a occupy contiguous
locations in memory. This is true, as long as the allocation did not
fail, right?
On the other hand, suppose I have
subroutine test2(m,n)
integer :: m,n
double precision, dimension(:,:), pointer :: a
allocate(a(m,n))
!... omitted do-loops to initialize a
a(i,j)=(j-1)*n+i
!...
end subroutine test2
I assume that a succesful call to allocate creates not one but two
arrays: one array contains enough memory to store mn adresses, while the
other array can hold mn double precision numbers. Is this correct?
It seems that the first array must be contiguous in memory or address
calculation would be hard, but can I be certain that the second array is
contiguous?
I am using the compiler ifort version 8.1.
My thanks
Carl Christian K. Mikkelsen
Gordon Sande
<cmik...@cs.purdue.edu> said:
Fortran does not have any way to DIRECTLY describe an array of pointers so
your declaration is just a pointer to a 2-d array. But what a pointer! It has
bounds and all sorts of thinngs that a C programmer would kill for.
(You question
seems to be based on knowing C or something like it.)
Fortran does allow for arrays of thingys, more formally know as derived
data types.
And the thingy can have a component which is a pointer. Most of the
time the Fortran
pointer just does what you want as a dynamic equivalence, to use older
fortran notions.
However it does byte when you try to do explicitly what you ask for.
What is the point
of what you explicitly asked for? Hopefully it was a simplification of
what you really
want. For realistic problems the awkwardness is about the same as other
languages.
Richard Maine
> subroutine test1(m,n)
>
> integer :: m,n
>
> double precision, dimension(:,:), allocatable :: a
> allocate(a(m,n))
>
> !... omitted do-loops to initialize a
> a(i,j)=(j-1)*n+i
> !...
>
> end subroutine test1
>
> then I am reasonably certain that the elements of a occupy contiguous
> locations in memory. This is true, as long as the allocation did not
> fail, right?
Yes, The standard doesn't guarantee memory locations per se, but
realistically that's going to be the case. It might help if you
mentioned why you cared. For the most part, a basic idea of many modern
constructs is that you aren't supposed to have to worry about exact
memory addresses; the language takes care of that. There are definitely
exceptions, where memory addresses can matter to you for one reason or
other, but if you want the best answer, you should specify why you care.
If you are concerned about performance issues, far better to investigate
the performance directly than to try to make guesses inferred from
storage layout. If you are concerned about C interoperability, I
recommend looking at the f2003 C interop features specifically aimed at
that.
> subroutine test2(m,n)
>
> integer :: m,n
>
> double precision, dimension(:,:), pointer :: a
> allocate(a(m,n))
>
> !... omitted do-loops to initialize a
> a(i,j)=(j-1)*n+i
> !...
>
> end subroutine test2
>
> I assume that a succesful call to allocate creates not one but two
> arrays: one array contains enough memory to store mn adresses, while the
> other array can hold mn double precision numbers. Is this correct?
No, no, no. Not even close. You seem to be thinking this would be an
array of pointers. It is not. It is *ONE* pointer to an array. The
double precision array will be continguous, just like the allocatable
case. Pointers can point to non-contiguous arrays, but that won't in
practice happen with an allocate statement; only with other ways of
assigning a pointer to a target.
Then there will be the pointer part, but that part will not be just an
address. It will have an address, plus information about bounds and
stride. It will have a fixed size, independent of the sie of the target.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
Carl Christian Kjelgaard Mikkelsen
I have had problems in the past which I could not precisely nail down,
because I was writting an MPI code calling LAPACK. I thought something
went wrong when I passed submatrices to, say, the Gaussian elimination
routine, but no it seems to be working perfectly, as in
program pt
implicit none
integer :: m, i, j, info
integer :: ipiv(2)
double precision, dimension(:,:), pointer :: a
character(len=9) :: myformat='(20F10.4)'
m=4
allocate(a(m,m))
do i=1,m
do j=1,m
a(i,j)=j**i
end do
end do
call dgetrf(2,2,a(1:3:2,1:3:2),2,ipiv,info)
do i=1,m
print myformat, (a(i,j), j=1,m)
end do
end program pt
which generates a 4 by 4 Vandermonde matrix, LU factors a particular 2
by 2 submatrix and returns
1.0000 2.0000 3.0000 4.0000
1.0000 4.0000 9.0000 16.0000
1.0000 8.0000 24.0000 64.0000
1.0000 16.0000 81.0000 256.0000
i.e. the correct result
1 3
1 24
embedded in the original matrix.
However, I have plenty of other possibilities for bugs.
Thanks
CCKM
Ron Shepard
Carl Christian Kjelgaard Mikkelsen <cmik...@cs.purdue.edu> wrote:
> call dgetrf(2,2,a(1:3:2,1:3:2),2,ipiv,info)
You probably could use a loc() function, or look at the assembler
code, to see what is happening here, but it is almost certainly true
that the compiler is generating local copies of this array and doing
copy-in/copy-out. This is because the subroutine has only an
implicit interface, so it must assume that adjacent array elements
must be contiguous in memory. On the other hand, if the subroutine
had an explicit interface, and that interface used implicit shape
dummy arrays, then those copies would probably not occur (they
might, but more likely not).
$.02 -Ron Shepard