Y = M*X 
where M is a sparse matrix stored by Compressed Row Storage format(CRS)

//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>//
void sparse_matrix_multiplication(argument n...)
{
        for(int i=0; i<size; i++)
        {
           Y[i] += M[ia[i]] * X[i];
	}//diagonal element

	for(int i=0; i<size; i++)
	{
		for(int j=ia[i]+1; j<ia[i+1];j++)
		{
			Y[i] += M[j] * X[ja[j]];  //Upper element
			Y[ja[j]] += M[j] * X[i];  //lower element
		}
	}
}
//<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<//

this code run on the CPU gets 5.213ms(0.38 Gflops).

(M : half bandwidth is 404 and size is 20808 with 518436 nonzero elements)