您不能使用全数组操作来实现这一点，因此您需要一个循环。

integer :: m,n
integer :: x((m-1)*(n-1))
integer :: y((m-1)*n+m*(n-1))

integer :: offset, block
integer :: xstart, ystart

offset = n*(m-1)
block = m-1

y = 0
y(:n*block) = y(:n*block) + x
y(m:(n+1)*block) = y(m:(n+1)*block) - x
do i=1,n-1
  xstart = (m-1)*(i-1)+1
  ystart = offset+m*(i-1)+1
  y(ystart  :ystart+block  ) = y(ystart  :ystart+block  ) - x(xstart:xstart+block)
  y(ystart+1:ystart+block+1) = y(ystart+1:ystart+block+1) + x(xstart:xstart+block)
enddo

在矩阵中存储零从来都不是一个好主意。如果是这样，就意味着您正在使用的数据结构对于您试图解决的问题来说不是最优的。
由于您的非方矩阵结构是由对角线或带状块组成的，似乎其中的大多数元素都将始终为零。我建议您使用像压缩稀疏行(CSR)这样的稀疏格式来存储矩阵(例如，参见Fortran中的稀疏矩阵向量积)。在这个例子中，转置后的矩阵向量积看起来非常简单

function transposed_matvec(A,x) result(b)
   type(CSRMatrix), intent(in) :: A
   real(real64), intent(in) :: x(:) 
   real(real64) :: b(A%n),aij
   integer :: j1,j2,row,col

   if (size(x)/=A%m) stop ' transposed_matvec: invalid array size '

   ! Initialize b
   b = 0.0_real64  
 
   ! Compute product by rows because data is row-ordered
   do row=1,A%m
     j1 = A%rowPtr(row)
     j2 = A%rowPtr(row+1)-1; 

     do j=j1,j2
        col = A%colPtr(j)
        aij = A%aij(j)
        b(col) = MxV(col) + aij*x(col) 
     end do
   end do   

end function transposed_matvec

如果改为使用压缩稀疏列存储(CSC)，则矩阵将已经按列存储，因此转置后的matvec例程看起来与链接中示例中CSR的直接matvec例程相同。