将块对角矩阵乘以向量的有效方法答案

【问题标题】：efficient way of multiplying block diagonal matrix by vector将块对角矩阵乘以向量的有效方法
【发布时间】：2022-08-21 14:00:01
【问题描述】：

我有一个矩阵 C 结构如下：

需要将其转置乘以向量x。

上半部分清晰 - 取矢量前半部分的切片说：

假设索引从 1 开始。

x1 = x(1:(n-1)*(m-1))

x2 = -x(m:n*(m-1))

然后部分增加：

x(1:(n-1)*(m-1)) += x1

x(m:n*(m-1))+=x2

但是如何处理下部（转置后的左侧）部分？有什么建议么？

这些矩阵有多大？这里没有很多非零元素。你在使用sparse 矩阵吗？这是您代码中的瓶颈吗？ x 是水平的还是垂直的？我，e，是x*C\'还是C\'*x？
@StewieGriffin 很大。 ‘x’是向量，所以是垂直的。没有“稀疏”的意义，因为矩阵的结构是已知的。甚至不需要存储矩阵。

标签： algorithm performance matlab fortran coding-efficiency

【解决方案1】：

你不能用整个数组操作来做到这一点，所以你需要一个循环。

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

【讨论】：

稍微改变了你的，某事是错误的，让我更新
看起来像列方向的块大小是m，因为转置
nvm，它的 m-1 实际上你是对的
哦，是的，哎呀，我忘了转调:)

【解决方案2】：

在矩阵中存储零绝不是一个好主意。如果这样做，则意味着您使用的数据结构对于您要解决的问题不是最佳的。

由于您的非方阵结构由对角线或带状块组成，因此大多数元素似乎始终为零。我建议您使用 Compressed Sparse Row (CSR) 之类的稀疏格式存储矩阵（例如，参见 Sparse Matrix vector product in 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 例程相同。

【讨论】：