# GPU Algorithms for Structured Sparse Matrix Multiplication with Diagonal Storage Schemes +
## 6 Algorithms  
### Banded Matrix-Matrix Multiplication  
We implemented banded matrix-matrix multiplication using the CDM data structure. Banded matrices that arise in real-life applications usually have narrow bands compared to their dimensions [37], which means that the CDM representation of the nonzero diagonals of narrow banded matrices are packed with only a small number of zeros. Thus, the compact diagonal scheme is more suitable for this case. On the other hand, we implemented structured sparse matrix-matrix multiplication using the HM data structure, which is described in the next section. As a banded matrix is a special kind of structured sparse matrix, the implementation that we are going to describe in the next subsection can also be used for banded matrix-matrix multiplication. In the case that the width of a banded matrix is large with respect to \(n\), the number of zeros for padding increases, so the HM data structure is more suitable than CDM. In our GPU algorithm, a thread will compute one entry of matrix \(C\). As matrix \(C\) is also a banded matrix, it is natural to assign a 2-D thread block of size \(T\times T\) to compute \(T^{2}\) entries of \(C\). So, the total number of threads in a thread block is \(T^{2}\), which should be a multiple of warp size. We choose \(T=2^{d}\), where \(d=4\). Before describing GPU multiplication in greater detail, we first provide a high-level overview of the overall algorithm. Figure 4 displays the logical mapping of GPU threads to the computed entries of the product matrix \(C\). A row corresponds to a compact diagonal of matrix \(C\), while the entries of matrix \(C\) are computed in blocks (tiles) as shown. Suppose the elements in the highlighted tile of the \(C\) matrix of Figure 5 are to be computed using the corresponding thread block. The tiles marked \(0,1,2,3\) from matrix \(A\) are needed to compute the elements in the highlighted tile of \(C\). The GPU algorithm works in a block-iterative manner. In the \(i\)-th block iteration, where \(0\leq i<4\), of Figure 5, a thread block brings the \(i\)-th tile of matrix \(A\) to the shared memory, performs all required multiplications between the tile elements from \(A\) and the nonzeros from \(B\), such that the multiplications contribute to the elements in the highlighted tile of matrix \(C\). The algorithm terminates once all the tiles of matrix \(\mathrm{Cb}\) are computed. A banded matrix \(X\) is represented by the following two lists and one integer. Here, a diagonal of \(X\) will be stored as a row in the CDM data structure. As all the diagonals in \(X\) are centered around the principal diagonal, in the CDM storage scheme, the (short) diagonal of zeros can be appended either at the beginning or at the end of the nonzero diagonals. 1. Integer \(dgX\) representing the number of compact diagonals in \(X\). 2. List of integers \(DsX\) of length \(dgX\), which is the list of the indices of the diagonals of \(X\). 3. _Xlist_: One-dimensional list of length \(dgX\times n\) containing the values of \(X\). We use \(nzX\) to denote the number of nonzeros in matrix \(X\). Given two banded matrices \(A\) and \(B\), matrix \(C\) will also be a banded matrix whose compact diagonals are computed via the following formula. Compact diagonal \(i\) of \(C\) is computed by multiplying the \(k\)-th compact diagonal of \(A\) and the \((i-k)\)-th compact  
process, we store the starting index of each diagonal in \(Blist\) in a dense format. We pass this information as a separate array called _BdiagDense_ to our algorithm. The algorithm is given in Algorithm 1.