efficient way of multiplying block diagonal matrix by vector

user_3121723 0
Viewed 105

I have a matrix C structured as following: C

Need to multiply its transpose by vector x.

with the upper part its clear - take slices of the first half of the vector say:

suppose indexation starts at 1.

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

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

then increment partially:

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

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

but how to deal with the lower (left after transpose) part? any suggestions?

coding-efficiency fortran matlab algorithm performance
2 Answers

You can't do it with whole-array operations, so you will need a loop.

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
user_85593000

It's never a good idea to store zeros in a matrix. If you do, it means that the data structure you're using is sub-optimal for the problem you're trying to solve.

Since your non-square matrix structure is made of diagonal or banded blocks, it seems like most of it elements will always be zero. I suggest you store your matrix using a sparse format like Compressed Sparse Row (CSR) (see for example Sparse Matrix vector product in Fortran). With that example, a transposed matrix-vector product would look pretty simply


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

If you used compressed sparse column storage (CSC) instead, the matrix would be stored by columns already, so the transposed matvec routine would look identical to the direct matvec routine for CSR at the example in the link.

user_47028700
Related