## 4 SURAA: The Proposed Method and Tool
### Setup Phase  
Algorithm 2 depicts the setup phase where the matrix is sorted, segmented, and grouped into various groups. The matrix is sorted row by row in the descending order of the number of non-zero elements per row. The rows with the same number of number of non-zeros are further sorted based on the column index of the first non-zero element, which further balances the load of vector \(x\). The sorting can be done on the GPUs or the CPUs. In our work, we use a GPU based parallel radix sort to perform the sort. The function \(radixSort\_Key(...)\) is used to perform the sorting of the \(nprv\) vector. It sorts the \(nprv\) vector and stores the permutations of the original \(nprv\) vector in the \(perm\) vector and the associated \(npr\) values in the \(nprSort\) vector. The input matrix is then rearranged using the \(perm\) vector. ```
0:\(A_{csr}(val,row_ptr,col_idx)\), \(vector:x\), \(vector:perm\), \(vector:nprv\)
0:\(set:G_{1}\), \(set:G_{2}\), \(set:G_{3}\)
1:functionSetup_Phase
2:Initialisation :\(perm\in seq[0,m]\)
3:\(nnzpr\leftarrow\) non-zero elements per row
4:\(nprSort,perm\leftarrow\)\(radixSort\_Key(nprv,perm)\)\(\triangleright\) Sort Descending
5: Rearrange \(A_{csr}\) using \(perm\)
6:\(q12\leftarrow\)\(nprSort[m*1/4]\)
7:\(q34\leftarrow\)\(nprSort[m*3/4]\)
8:\(Seg\_size\leftarrow(2\times(q34-q12))/\sqrt{m})\)
9:for Each segment \(s\)do
10:if\(32<nprv[s_{i}]<max_child_threads\)then
11:\(G_{1}\gets G_{1}\cup s_{i}\)
12:elseif\(nprv[s_{i}]\geq 32\)then
13:\(G_{2}\gets G_{2}\cup s_{i}\)
14:else
15:\(G_{3}\gets G_{3}\cup s_{i}\)
16:endif
17:endfor
18:endfunction
```  
**Algorithm 2** SURAA_Setup_Phase  
Secondly, the sorted matrix is divided into various segments. The division into segments is based on the Freedman-Diaconis rule [PERSON] (1992); [PERSON] (1992), which is given below:  
\[Seg\_size\gets 2\times(IQR(nprv)/\sqrt[3]{m}), \tag{1}\]  
where \(Seg\_size\) is the size of each segment, \(IQR\) is the Interquartile range, \(nprv\) has been defined already, and \(m\) is the number of rows in the matrix. Moreover, \(IQR(nprv)=q_{34}-q_{12}\), where \(q_{34}\) and \(q_{12}\) are the third and first quartiles of \(nprv\), respectively. In our algorithm, we calculate the quartiles of the sorted \(npr\) array (\(nprSort\)). Since the matrix is pre-sorted, we can get the optimal size of the segments as well as the number of segments using the Freedman-Diaconis rule for each matrix with \(\mathcal{O}(1)\) complexity. In the third step, the segments are allocated to various groups depending upon the average number of non-zero elements per row in each segment. This is also depicted in the left part of Figure 5. The number 0 to N indicates the \(npr\) for each row and the numbers 0 to M indicate the number of rows in the matrix. The red dashed line in the figure indicates the size of the warp (0 to 31) for GPUs, which is used as a threshold in our method. Since the matrix is sorted in the descending order of \(npr\), we add each segment \(s_{i}\) from the beginning of the matrix to Group 1 (\(G_{1}\)) until the mean number of non-zero per row for a segment becomes less than 32 or the total number of rows in \(G_{1}\) exceeds \(max\_child\_threads\) (The number of child grids in the queue ready for execution). Since the segment size might not be a multiple of group size, if \(n\times segment\_size>max\_child\_threads\), we will only include \((n-1)\times segment\_size>max\_child\_threads\). The rest of the segments which has mean \(npr\) greater than or equal to 32 is allocated to \(G_{2}\). The remaining segments that have a mean number of non-zero per row less than 32 are added to \(G_{3}\). The last segment might not have the same number of rows as the other segments and hence is executed separately. All the segments in \(G_{3}\) are executed asynchronously on different streams using a scalar CSR kernel, whereas, all the segments in \(G_{2}\) are executed by a modified vector CSR kernel asynchronously.