Consider the following (traditional) Fortran code, which performs matrix-vector multiplication z = X y
, but instead of the full matrix, we exclude the top two rows:
subroutine matrix_vector_product(m, n, x, y, z)
integer :: m, n
real :: x(m, n), y(n), z(m-2)
interface
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
character :: TRANS
integer :: INCX, INCY, LDA, M, N
real :: ALPHA, BETA, A(LDA,*), X(*), Y(*)
end
end interface
call SGEMV("N", m-2, n, 1.0, x(3, 1), m, y, 1, 0.0, z, 1)
! ^~~~~~~ we pass a scalar here.
end
Note that we pass a scalar as argument where actually the matrix is required. This works because of a concept called sequence association (Fortran scalars are passed by pointer and arrays are passed by pointer-to-first-element, so we can simply pass the start position of the array).
However, if the interface block is present, gfortran emits a diagnostic, because this form of sequence association is technically illegal (you are not allowed to pass a scalar in place of an array). The other problem is that x(1,3)
is not a valid element for m = 2
, even though the matrix-vector product stays well-defined. -fcheck=bounds
then complains and aborts the program.
There are a couple of workarounds:
...
call SGEMV("N", m-2, n, 1.0, x(3:, 1:), m-2, y, 1, 0.0, z, 1)
Annoyingly, this creates a temporary, because the interface requires a contiguous array, which the slice is not.
...
real, contiguous, pointer :: px_flat(:), px(:, :)
px_flat(1:m*n) => x
px(1:m, 1:n) => px_flat(3:)
call SGEMV("N", m-2, n, 1.0, px, m, y, 1, 0.0, z, 1)
This is the most complicated solution of them all: we associate a pointer with the sequence of elements and then reshape that pointer into the desired matrix, which we then pass into the subroutine. This is verbose and also somewhat hacky: px
does not really have the "true" shape of the array ... but at least it works and does not cause gfortran to complain.
Is there a simpler/better solution?
If the desired array section is discontiguous, sequence association isn't going to work easily when there are multiple columns. I don't emit an error or warning for this example from flang-new, since it's a valid usage of sequence association, but please be careful.