## 2 Methods
### CUDA Implementation of Greedy Agglomerative Procedures
\(D_{j}\gets 0\forall j=\overline{1,p}\)// \(D_{j}\) are used to calculate \(\sum\limits_{j\in C_{j},L(A_{j},X_{i})>\epsilon_{2}}\frac{1}{L(X_{i},A_{j})}\). \(D^{\prime}_{j}\gets 0\forall j=\overline{1,p}\)// \(d\)-dimensional vectors to calculate \(\sum\limits_{j\in C_{j},L(A_{j},X_{i})>\epsilon_{2}}\frac{A_{j}}{L(X_{i},A_{j})}\). ```

**Algorithm 4** CUDA implementation of Step 1 in the ALA algorithm, part 1 (initialization)

For the rest of the CUDA implementation of Step 1 of the ALA algorithm, we used the number of threads \(N_{threads}=512\) for each CUDA block. The number of blocks is calculated as

\[N_{blocks}=(N+N_{threads}-1)/N_{threads}. \tag{9}\]

Thus, each thread processes a single demand point in Algorithm 5:

```
\(i\gets blockIdx.x\)times\(timesblockDim.x+threadIdx.x\) if\(i>N\)then return endif \(j^{\prime}\leftarrow\operatorname*{argmin}_{j}\parallel A_{j}-X_{i}\parallel^{2}\)// number of a center for the \(i\)th data point \(d^{\prime}\leftarrow\min_{j}\parallel A_{j}-X_{i}\parallel^{2}\); \(C_{i}\gets j^{\prime}\)// Assign \(A_{i}\) to center \(j^{\prime}\). \(D\leftarrow(\sqrt{d^{\prime}})^{-1}\); \(D_{j}\gets D_{j}+D\); \(D^{\prime}_{j}\gets D^{\prime}_{j}+A_{j}D\); \(counter_{j^{\prime}}\gets counter_{j^{\prime}}+1\);  Synchronize threads. ```

**Algorithm 5** CUDA implementation of Step 1 in the ALA algorithm

The parallel implementation of Step 2 of the ALA algorithm (see Algorithm 6) uses \(N_{threads}=512\) threads for each CUDA block. The number of blocks is \(N_{blocks2}=(p+N_{threads}-1)\) /\(N_{threads}\). The changed variable, common to all threads, is initialized to 0. ```
\(j\gets blockIdx.x\times blockDim.x+threadIdx.x\) if\(j>k\)then return; endif; \(X^{\prime}_{j}\leftarrow\frac{D^{\prime}_{j}}{A_{j}}\); if\(L\Big{(}X^{\prime}_{j},X_{j}\Big{)}\geqslant\epsilon_{2}\)then \(changed\gets 1\); endif; \(X_{j}\gets X^{\prime}_{j}\);  Synchronize threads. ```

**Algorithm 6** CUDA implementation of Step 2 in the ALA algorithm

Variable changed is used at the last step of the ALA algorithm and signals the need to continue the iterations. In addition, we implemented Step 3 of BasicAgggl algorithm on the GPU. At this step, Algorithm 2 calculates the sum of the distances after removing a single center with index \(i^{\prime}\): \(F^{\prime}_{i^{\prime}}=F(S\backslash\{X_{i^{\prime}}\})\). Knowing the objective function value \(F(S)\) for set \(S\) of centers, we can calculate its new value

\[F^{\prime}_{i^{\prime}}\gets F(S\backslash\{X_{i^{\prime}}\})=F(S)+\sum \limits_{l=1}^{N}\Delta D_{l} \tag{10}\]

where\[\Delta D_{l}=\left\{ \begin{array}{c}0,\quad\quad C_{i^{\prime}}\
eq l,\\ \left(\min_{j\in\{\frac{1}{l,p\},\,j\
eq l^{\prime}}L\left(A_{j},X_{C_{j}} \right)\right)}-L\left(A_{j},X_{i}\right),C_{i^{\prime}}=l.\end{array} \right. \tag{11}\]

Here, \(C_{i^{\prime}}\) is the center number which \(A_{i^{\prime}}\) is assigned to (nearest to \(X_{i^{\prime}}\)). We also used 512 threads per block, the number of blocks is calculated according to (9). After initialization of common variable \(Dsum\)\(\leftarrow\)0, the Algorithm 7 is started in a separate thread for each data point:

```
\(l\gets blockIdx.x\)times\(blockDim.x+threadIdx.x\); if\(l>k\)then return; endif;  Calculate \(\Delta D_{l}\) in accordance with (2.8); if\(\Delta D_{l}>0\)then \(D_{sum}\gets D_{sum}+\Delta D_{l}\); endif Synchronize threads. ```

**Algorithm 7** CUDA implementation of Step 3 in Algorithm 2

After parallel running this algorithm, \(F_{i}\) is calculated in accordance with (10) on Step 3 of Algorithm 2: \(F^{\prime}_{i^{\prime}}\gets F(S)+D_{sum}\). The remaining parts of our GA are implemented by the central processor unit. The remaining parts performing the selection do not have any significant effect on the calculation speed. This CUDA implementation of the ALA algorithm enables us to use both the ALA algorithm and more complex algorithms that include it, on large amounts of data, up to millions of demand points.