## 2 Methods
### CUDA Implementation of Greedy Agglomerative Procedures  
Compute Unified Device Architecture (CUDA) is a hardware-software parallel computing architecture that can significantly increase computational performance [98] designed for the effective use of non-graphical computing on the graphics processing units (GPUs). A GPU is a set of multiprocessors that simultaneously execute parallel threads grouped into data blocks. The threads execute the same instructions on different data in parallel. Communication and synchronization in blocks is impossible. A thread is assigned an identifier (_threadIdx.x_) depending on its position in the block, and a block identifier (_blockIdx.x_) is also assigned to a thread depending on its position in the grid. The thread and block identifiers are available at run time, which allows us to set specific memory access patterns. Each thread in the GPU performs the same procedure, known as the kernel [98]. CUDA versions of the algorithms for solving discrete p-median problems [99; 100; 101] (which actually solve rather small discrete problems from the OR (Operation Research) Library) as well as parallel implementations of the ALA algorithm for the k-means problems are described in the modern literature [74; 98; 102]. In our research, we take into account the specificity of the continuous problems. For large amounts of data, calculating the distances \(L\big{(}X_{i},A_{j}\big{)}\) from demand points to cluster centers is the most computationally expensive part of the ALA algorithm. The distances should be calculated to all centers, since the algorithm does not know in advance which center is the nearest one. Traditionally used in solving the k-means problem, this approach can be accelerated by using the triangle inequality [103; 104] in the case of p-median problems: if the displacement of the ith center relative to the previous iteration was \(\Delta X_{i}=L\big{(}X_{i},X_{i}^{\prime}\big{)}\), and its distance to the jth demand point was \(L_{ij}=L\big{(}X_{i},A_{j}\big{)}\), then the new distance lies within \(L_{ij}\pm\Delta X_{i}\). If \(\exists i^{*}\in\big{\{}\overline{1,p}\big{\}}\):  
\[L_{i^{*}j}+\Delta X_{i^{*}}<L_{ij}-\Delta X_{i}, \tag{7}\]  
such that for new distances, \(L\big{(}X_{i^{*}}^{\prime},A_{j}\big{)}<L\big{(}X_{i}^{\prime},A_{j}\big{)}\). Knowing all distances \(L_{ij}\) from the previous iteration and the shift of the centers \(X_{i}\), we can avoid calculating a significant part of the distances at each iteration, which is especially important when processing multidimensional data. Nevertheless, for large-scale problems, matrix \(\big{\{}L_{ij}\big{\}}\) requires a significant additional amount of memory (\(pN\) values). Moreover, checking the condition (7) for each pair \(\big{(}A_{j},X_{i}\big{)}\) which requires significantly less computational resources than distance calculation  
\[L_{ij}=L\big{(}A_{j},X_{i}\big{)}=\sqrt{\sum_{l=1}^{d}\Big{(}a_{j,l}-x_{i,l} \Big{)}^{2}} \tag{8}\]  
will not accelerate the execution of the algorithm in the case of its CUDA implementation, which assumes strictly simultaneous execution of all instructions by all threads. Therefore, in this work, we do not use of the triangle inequality to accelerate the CUDA version of the ALA algorithm, although more complex implementations for CUDA could probably take advantage of the triangle inequality. Note that the distances (8) are calculated both in Step 1 of Algorithm 1, and in its Step 2 which includes the iteration of the Weiszfeld procedure (4). We used the following CUDA implementation of Step 1 of Algorithm 1. The step is divided into two parts. First, variables common for all threads are initialized Algorithm 4:```
\(X^{\prime}_{j}\gets 0\forall j=\overline{1,p}\)// Here, \(X^{\prime}_{j}\) are vectors used for recalculation of centers. \(counter_{j}\gets 0\forall j=\overline{1,p}\)// Counters of objects assigned to each of centers.