# GPU Algorithms for Structured Sparse Matrix Multiplication with Diagonal Storage Schemes +
## 6 Algorithms
### Structured Sparse Matrix-Matrix Multiplication
We use \(nzX\) to present the number of nonzeros in sparse matrix \(X\). A thread that handles all multiplications associated with an entry from \(A\) will run in a loop for each diagonal of \(B\). An entry from \(A\) can be from a super-, sub-, or principal diagonal of \(A\). In the same way, a diagonal from \(B\) can also be a super-, sub-, or principal diagonal of \(B\). We enumerate thirteen different cases based on the different choices of the diagonals of \(A\) and \(B\). These thirteen different cases are grouped into three categories. Ineach of the thirteen cases, given an entry from \(A\) and a diagonal from \(B\), our objective is to find the entry in the diagonal of \(B\) that can be multiplied with the entry of \(A\). After a multiplication, the thread will find the entry in \(C\) to which this multiplication will contribute. In our descriptions of the thirteen cases below, we will use three index variables \(i,j,k\), where \(n>i,j,k\geq 0\). For an entry of \(A\), it is possible that there might be no entry from a particular diagonal of \(B\) with which it can be multiplied. In short, we describe the situation by saying that such an entry in \(B\) does not _exist_. In the first category of cases, we describe the cases that can arise while computing the principal diagonal of \(C\). 1. The multiplications between the principal diagonals of \(A\) and \(B\) contribute to the principal diagonal of \(C\). Each multiplication can be described in the following way, with the left index \(0\) indicating the (principal) diagonal: \[C[0][i]=A[0][i]*B[0][i].\]
2. The entries of \(A\) and \(B\) are from the \(i\)-th superdiagonal and \((-i)\)-th subdiagonal, respectively. Each of the multiplications can be described in the following way: \(C[0][j]=A[i][j]*B[-i][j]\). 3. The entries of \(A\) and \(B\) are from the \((-i)\)-th subdiagonal and \(i\)-th superdiagonal, respectively. Each of the multiplications can be described in the following way: \(C[0][i+j]=A[-i][j]*B[i][j]\). In the second category, we have five different cases that arise while computing the entries in the superdiagonals of \(C\). 1. Let \(A[i][k]\) and \(j\) be an entry from \(A\) and an index of a diagonal from \(B\), respectively. If \(B[j][i+k]\) exists, then \(B[j][i+k]\) will be multiplied with \(A[i][k]\), and this multiplication will contribute to \(C[i+j][k]\). 2. Let \(A[0][i]\) and \(j\) be an  
\(A\) and \(B\)) contribute to a superdiagonal of \(C\) (second category case 1). In the same way, by changing the inequality among \(i,j\), and \(k\), we can verify the rest of the 12 cases. Algorithm 2 has a number of nested _if...else_ blocks, which in practice can increase the code divergence among threads in the same thread block. In the worst case, this code divergence might significantly affect performance. One effective solution to overcome this code divergence issue is to have a few threads (for example, \(k\), where \(1<k\leq 13\)) per nonzero entry of \(A\) that will take care of all multiplications associated with that nonzero entry of \(A\). The _if...else_ part of the code can be restructured such that each of the _ifcondition_)...is will be mapped to one of the \(k\) threads exclusively. As all threads in a thread block will take care of the same _ifcondition_)(...), we expect to have lower code divergence. ```
__global__voidsparseStructured( MatrixA, MatrixB, MatrixC){ int i = blockIdx.x * T + threadIdx.x; //nzAisthenumberofnonzerosinA if(i < nzA){ //FindDiagonalByBSusesbinarysearch Ad = FindDiagonalByBS(StartDgA,i) //indisttheindexofAlist[i] //inAd-thdiagonal for(eachdiagonalBdofB){ Cd = Ad + Bd //computethenonzeroindices //iandj j = FindBInd(Ad,Aind,Bd) k = FindCInd(Ad,Aind,Bd,j) //thefollowingupdateisperformed //viaanatomicoperation. CList[k] +=AList[i]*BList[j] } } }
```  
**Algorithm 2**CUDA Algorithm for Structured Sparse Matrix-Matrix Multiplication