Usually you would do one of two things:

1. Compute $z = Ax$, then $y = y+Az$;
2. If you anticipate needing to do this operation many times, you might compute the matrix $B = AA'$ and then evaluate $y = y+Bx$, but it will be much less sparse than the original matrix $A$.

I don't know of any libraries that will compute $AA'x$ directly, and if there are any that do so, they probably use (1) behind the scenes.

SPARSKIT has routines for both $A\cdot x$ and $A'\cdot x$, they're called amux and atmux respectively; these are for matrices in the compressed sparse row format, which is described in the SPARSKIT documentation.

Do you have a language preference? If you're just experimenting and don't need your program to be super fast, you could try the Scipy sparse matrix modules in python, which would be less of a hassle than using Fortran. These will also have the functionality you're looking for.

Usually you would do one of two things:

1. Compute $z = Ax$, then $y = y+Az$;
2. If you anticipate needing to do this operation many times, you might compute the matrix $B = AA'$ and then evaluate $y = y+Bx$, but it will be much less sparse than the original matrix $A$.

I don't know of any libraries that will compute $AA'x$ directly, and if there are any that do so, they probably use (1) behind the scenes.

SPARSKIT has routines for both $A\cdot x$ and $A'\cdot x$, they're called amux and atmux respectively; these are for matrices in the compressed sparse row format, which is described in the SPARSKIT documentation.

Do you have a language preference? If you're just experimenting and don't need your program to be super fast, you could try the Scipy sparse matrix modules in python, which would be less of a hassle than using Fortran. These will also have the functionality you're looking for.

EDIT: Since you want fast and you mentioned GPU computing, you might want to consider the CUSP library. The ellpack matrix format is particularly amenable to usage on GPU machines, although there are exceptional cases.

Usually you would do one of two things:

1. Compute $z = Ax$, then $y = y+Az$;
2. If you anticipate needing to do this operation many times, you might computed $B = AA'$ and then evaluate $y = y+Bx$.

I don't know of any libraries that will computed $AA'x$ directly, and if there are any that do so, they probably use (1) behind the scenes.

SPARSKIT has routines for both $A\cdot x$ and $A'\cdot x$, they're called amux and atmux respectively; these are for matrices in the compressed sparse row format, which is described in the SPARSKIT documentation.

Do you have a language preference? If you're just experimenting and don't need your program to be super fast, you could try the sparse matrix modules from the numpy package in python, which would be less of a hassle than using Fortran. These will also have the functionality you're looking for.

Usually you would do one of two things:

1. Compute $z = Ax$, then $y = y+Az$;
2. If you anticipate needing to do this operation many times, you might compute the matrix $B = AA'$ and then evaluate $y = y+Bx$, but it will be much less sparse than the original matrix $A$.

I don't know of any libraries that will compute $AA'x$ directly, and if there are any that do so, they probably use (1) behind the scenes.

SPARSKIT has routines for both $A\cdot x$ and $A'\cdot x$, they're called amux and atmux respectively; these are for matrices in the compressed sparse row format, which is described in the SPARSKIT documentation.

Do you have a language preference? If you're just experimenting and don't need your program to be super fast, you could try the Scipy sparse matrix modules in python, which would be less of a hassle than using Fortran. These will also have the functionality you're looking for.