Multidimensional arrays are stored in column-major order. This means the
left-most (inner-most) index addresses elements contiguously.
From a practical point this means that the array slice
V(:, 1) is
contiguous, while the stride between elements in the slice
is the dimension of the columns. This is important when passing array
slices to procedures which expect to work on contiguous data.
The locality of the memory is important to consider depending on your application, usually when performing operations on a multidimensional the sequential access should always advance in unity strides.
In the following example the inverse distance between two sets of points
is evaluated. Note that the points are stored contiguously in the arrays
xyz2, while the inner-most loop is advancing the left-most
index of the matrix
subroutine coulomb_matrix(xyz1, xyz2, a) real(dp), intent(in) :: xyz1(:, :) real(dp), intent(in) :: xyz2(:, :) real(dp), intent(out) :: a(:, :) integer :: i, j do i = 1, size(a, 2) do j = 1, size(a, 1) a(j, i) = 1.0_dp/norm2(xyz1(:, j) - xyz2(:, i)) end do end do end subroutine coulomb_matrix
Another example would be the contraction of the third dimension of a rank three array:
do i = 1, size(amat, 3) do j = 1, size(amat, 2) do k = 1, size(amat, 1) cmat(k, j) = cmat(k, j) + amat(k, j, i) * bvec(i) end do end do end do
Contiguous array slices can be used in array-bound remapping to allow usage of higher rank arrays as lower rank arrays without requiring to reshape and potentially create a temporary copy of the array.
For example this can be used to contract the third dimension of a rank three array using a matrix-vector operation:
subroutine matmul312(amat, bvec, cmat) real(dp), contiguous, intent(in), target :: amat(:, :, :) real(dp), intent(in) :: bvec(:) real(dp), contiguous, intent(out), target :: cmat(:, :) real(dp), pointer :: aptr(:, :) real(dp), pointer :: cptr(:) aptr(1:size(amat, 1)*size(amat, 2), 1:size(amat, 3)) => amat cptr(1:size(cmat)) => cmat cptr = matmul(aptr, bvec) end subroutine matmul312