fortranparameter-passingslice

Sequence association vs. modern Fortran interfaces


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:

  1. Remove the interface block. We go back to the days of the wild west, where you hope and pray you got all the 10+ arguments to BLAS/LAPACK functions correctly. Not a nice debugging experience if thing do go wrong, and it does not solve the invalid element problem.
  2. Pass a slice of the array:
        ...
        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.

  3. Fortran 2008 pointer magic:
        ...
        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?


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.