# GPU Algorithms for Structured Sparse Matrix Multiplication with Diagonal Storage Schemes +
## 6 Algorithms
### Structured Sparse Matrix-Matrix Multiplication  
We consider two main techniques for designing a structured sparse matrix-matrix parallel multiplication algorithm targeting GPU architecture: (i) each entry of \(C\) is computed using one thread, and (ii) all multiplications that are associated with an entry of either \(A\) or \(B\) are performed using one thread. Both of the above approaches have a similar outlook from a data locality point of view. That is, we expect that multiple threads in a warp are accessing entries from matrices \(A\), \(B\), or \(C\) that are physically stored closely in the global memory. The first approach is a pure parallel approach. The second approach requires atomic operations since two different threads might contribute to the same entry of \(C\) at the same time. In this work, we implement this second approach and claim that the second approach minimizes computational overhead. Let \(i\) and \(j\) be two diagonal indices of matrices \(A\) and \(B\), respectively, where \(-n<i,j<n\). The multiplications between the elements of the \(i\)-th diagonal of \(A\) with the elements of the \(j\)-th diagonal of \(B\) will contribute to the \((i+j)\)-th diagonal of \(C\), where \(-n<(i+j)<n\). From this observation, we can expect that \(C\) is potentially more dense compared to \(A\) or \(B\). The number of multiplications required to compute an entry of \(C\) varies widely due to the irregular distribution of nonzero diagonals in matrices \(A\) and \(B\). In the first approach, such disproportionate thread workload creates code divergence in thread blocks, affecting the performance negatively. Moreover, a thread that is assigned to compute an entry from \(C\) needs to examine every possible pair of diagonals from \(A\) and \(B\). If a pair of diagonals from \(A\) and \(B\) exist (i.e., nonzero), then a multiplication is performed and the associated entry of matrix \(C\) is updated. As both \(A\) and \(B\) are sparse, there will be fewer nonzero diagonal pairs compared to the quadratic search space, thereby increasing the overhead. In the second approach, a thread is assigned to perform all multiplications associated with one entry from one of the operands \(A\) or \(B\). Without loss of generality, let \(A\) be the selected operand. The number of threads required to compute the matrix product \(C=AB\) is equal to the number of nonzeros in matrix \(A\) (i.e., \(nzA\)). A thread first computes the index of the entry from each nonzero diagonal in \(B\) with which the entry from \(A\) will be multiplied. After multiplication, the thread also updates the corresponding entry in \(C\) via an atomic operation. In this approach, we expect to experience less code divergence than that of the first approach. In summary, we conclude that although both approaches have similar data locality challenges and disadvantages, we prefer the second approach due to its favorable (expected) code divergence compared with the first approach. A structured sparse matrix (say, \(X\)) is represented by the following three lists and one integer value. Here, a diagonal of \(X\) will be treated as a diagonal in the HM data structure. 1. An integer value \(dgX\) that represents the number of HM diagonals in \(X\). In structured sparse matrices, we know that only some of the diagonals have nonzeros. We will store those diagonals only. 2. A list of integers \(DsX\) (of length \(dgX\)) that stores the indices of the diagonals of \(X\). The first element of this array is \(0\) if the principal diagonal is dense. It stores the indices of the super diagonals first, followed by those of the subdiagonals. For superdiagonal indices, those are stored in ascending order. For subdiagonal indices, those are stored in descending order. 3. _Xlist_: This is a one-dimensional array that stores the entries of \(X\). The entries from the same diagonal are kept together according to their column indices. 4. A list of integers \(StartDgX\) (of length \(dgX\)) that stores the indices of the entries in _Xlist_ for each diagonal of \(X\).