序列关联与现代Fortran接口：矩阵向量乘的更优方案？

考虑以下传统Fortran代码，它执行矩阵向量乘法z = X y，但排除了矩阵的前两行：

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)
    !                            ^~~~~~~ 此处传入的是标量
end

注意，我们在需要传入矩阵的位置传入了一个标量。这之所以可行，是因为一种名为**序列关联（sequence association）**的机制：Fortran中标量按指针传递，数组按指向首个元素的指针传递，因此可以直接传递数组的起始位置。

然而，如果存在接口块，gfortran会发出诊断信息，因为这种形式的序列关联在技术上是非法的（不允许用标量替代数组传递）。另一个问题是，当m = 2时，x(3, 1)不是有效元素，尽管此时矩阵向量乘法的定义仍然成立，启用-fcheck=bounds会报错并终止程序。

现有几种解决办法：

• 移除接口块：退回到无类型检查的阶段，需手动确保BLAS/LAPACK函数的10余个参数全部正确，调试体验不佳，且无法解决无效元素的问题。
• 传递数组切片：

  ...
    call SGEMV("N", m-2, n, 1.0, x(3:, 1:), m-2, y, 1, 0.0, z, 1)
  

  但这会创建临时数组，因为接口要求传入连续数组，而该切片并非连续存储。
• Fortran 2008指针技巧：

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

  这是最复杂的解决方案：将指针与元素序列关联，再将指针重塑为所需矩阵后传入子程序。该方法冗长且略显繁琐，px并不具备数组的“真实”形状，但至少可以正常运行且不会引发gfortran的报错。

是否存在更简单、更优的解决方案？

内容的提问来源于stack exchange，提问作者summentier