When should I write a matrix-vector function to handle the sparse matrix vector multiplication?

Question

This semster, I have been studying the iterative methods for large sparse matrix system. But I have some questions. For large sparse matrix, we must use an economic storage to store them. The most popular way is the compressed row storage or compressed column storage, i.e., CRS,CCS format which have been proved as the most economic way. When we use the iterative method to solve Ax=b, we must repeatedly compute the matrix-vector mulitiplication Ax. So, first we must store the sparse matrix $A$ in CRS,or CCS format, then we have to write a subroutine to compute the matrix-vector, which is so troublsome during our coding. But when I read the matlab routines about iterative methods, e.g., GMRES,MINRES,PCG, downloaded from website, I found that the author just use the matlab operation ' * ', to compute the matrix-vector multiplication, Ax, where matrix $A$ is a sparse format in matlab. For example, to compute $y=A*x$:

n=4;
A = gallery('poisson',n);
x=rand(n^2,1);
y=A*x

As for the similar code in matlab, I found that there is nobody write a single subroutine to compute the matrix-vector $A*x$. I do not understand that because in many monographys about iterative methods for large sparse matrix, they always discuss the sparse matrix storage and the economic way of computing matrix-vector. But I found no one do that when using matlab. Why this? because of using matlab? or other reasons?

If we do not use matlab, e.g., use C,Python, or Fortran language, can we still use the simple way y=A*x like matlab? any suggestions are welcome.

Answer 1

Multiplying matrix and a vector using sparse matrix format, such as CSR is not hard, and is a basic and common operation if you write a linear solvers library for sparse linear systems. Take SPARSEKIT of Yousef Saad as an example. In the code one of the first subroutines you will find is what he calls 'amux'. You can look at it below. As one can see it is fairly simple. Hope you will able to do the same in e.g. MATLAB.

So you should do a special sparse matrix - vector product function for this operation when you deal with sparse matrices, and probably they have it in MATLAB. Are they able to overload the '*' operator - probably they are, but I am not sure is this the case.

subroutine amux ( n, x, y, a, ja, ia )

!*****************************************************************************80
!
!! AMUX multiplies a CSR matrix A times a vector.
!
!  Discussion:
!
!    This routine multiplies a matrix by a vector using the dot product form.
!    Matrix A is stored in compressed sparse row storage.
!
!  Modified:
!
!    07 January 2004
!
!  Author:
!
!    Youcef Saad
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the row dimension of the matrix.
!
!    Input, real X(*), and array of length equal to the column dimension 
!    of A.
!
!    Input, real A(*), integer ( kind = 4 ) JA(*), IA(NROW+1), the matrix in CSR
!    Compressed Sparse Row format.
!
!    Output, real Y(N), the product A * X.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) a(*)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ia(*)
  integer ( kind = 4 ) ja(*)
  integer ( kind = 4 ) k
  real ( kind = 8 ) t
  real ( kind = 8 ) x(*)
  real ( kind = 8 ) y(n)

  do i = 1, n
!
!  Compute the inner product of row I with vector X.
!
    t = 0.0D+00
    do k = ia(i), ia(i+1)-1
      t = t + a(k) * x(ja(k))
    end do

    y(i) = t

  end do

  return
end

Answer 2

Thanks for all your reply, I have done some experiments in matlab and confirmed that matlab indeed uses the sparse matrix-vector multiplication automatically optimally. For example,

clc;clear;
rng(0);
n=160;
A = gallery('poisson',n);%  A is n^2 * n^2 and sparse matrix
x=rand(n^2,1);

tic;A*x;toc;

tic;full(A)*x;toc;

And the results are

Elapsed time is 0.001107 seconds.
Elapsed time is 7.549824 seconds.

my CPU is 8GB and matlab 2018b. From the above results, we can see that when coding in matab, the matrix-vector multiplication automatically chooses the best way to perform. So, it is so convient that we do not need to write a subroutine to compute Ax seperately.