Sequence association vs. modern Fortran interfaces

Question

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:

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.
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.
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?

Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Meta Stack Overflow, or in Stack Overflow Chat. Comments continuing discussion may be removed. — mickmackusa
– mickmackusa ♦, Commented Aug 19 at 9:13

Peter Klausler · Accepted Answer · 2025-08-17 17:05:52Z

2

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.

answered Aug 17 at 17:05

Peter Klausler

543 bronze badges

Sign up to request clarification or add additional context in comments.

Collectives™ on Stack Overflow

Sequence association vs. modern Fortran interfaces

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related