4.2.2 Operations with sparse matrices on a vector unit

Different storage schemes
$\Longrightarrow$ different implementations for the same operation.

Ex.: $\underline{y} = A_{n\times n}\cdot\underline{x}$ in CRS

     DO 
$i  =  1,  n$ 


$y(i) \;:=\; 0$ 


DO 
$j  =  row\_ptr(i),  row\_ptr(i+1)-1$ 


$y(i)\;=\;y(i) + val(j) \ast x(col\_ind(j))$


END DO 


END DO

How to vectorize the double loop ?

1. Blow up row $i$ of $A$ with Zeros
   $\Longrightarrow$ inner loop $\widehat{=}$ $y_i \;=\; y_i + (A_{i,\ast},\underline{x})$
   kontrary nonsense in 3D-FEM calculations.
2. GATHER/SCATTER (collect/distribute) a matrix row via column index vector.

   
        y = 0 
   
   
   DO 
   $i  =  1,  n$ 
   
   		      
   $length \;:=\; row\_ptr(i+1)+1 - row\_ptr(i)\qquad$              [Register length]
   
   
   $j1 \;:=\; row\_ptr(i)$ 
   
   
   $y(i) \;:=\;y(i) + val(j1) \ast x(col\_ind(j1))$
   
   
   END DO
   

   The vector unit loads index vector $col\_ind(j1)$ and calculates the proper entry in $\underline{x}$ (GATHER). The storage of a vector works similar (SCATTER).
   • Matrix-Matrix operations are performed in a similar way.
   • By using the indirect addressing it does not seem to give a possibility of avoiding memory access conflicts or page faults.
   • The parallelization of Matrix-by-Vector product should take into account the (local) structure of the matrix pattern (Sect. 4.3).