## 4 Experiments
Thus, in the Algorithm 3 (in p. 11), the value of \(\rho\) is equal to 1 because, in our implementation, we scale the original matrix so that we can always have that value. Vector \(\lambda\) elements are calculated as \((d_{k}-d_{l})-(\lambda_{i}-d_{l})\), with \(I=i\) or \(i+1\), to have good accuracy in floating point operations. This procedure is shown in Algorithm 4. ```
0:A vector \(d\in\mathbb{R}^{n}\) whose entries are in ascendant order without repeated values, a vector \(v\in\mathbb{R}^{n}\) with \(\left\lVert v\right\rVert_{2}=1\), the machine precision \(e\in\mathbb{R}\)
0:The \(n\) roots of the secular equation in vector \(\lambda\in\mathbb{R}^{n}\), a matrix \(drmlc\in\mathbb{R}^{n\times n}\) whose value in row \(r\) and column \(c\) is \(d_{r}-\lambda_{c}\)
1:for\(i=1\)a\(n\)do
2: Use Algorithm 3 with \(d\), \(v\), \(i\), \(e\) as input and \(\lambda_{i}\), the column \(i\) of \(drmlc\) as output
3:endfor
```

**Algorithm 4** A Comprehensive algorithm for solving the secular equation to determine its roots

Also, with this slight modification, we allow the CPU or GPU threads to have a more homogeneous execution since they have the same execution path. We implement both a sequential and a parallel version. The sequential version was programmed in C. The parallel version was programmed in C, CUDA, OpenMP, and MPI. Data structures in both of them use a small amount of memory. The tridiagonal matrix \(T\) uses two arrays: length \(n\) for the main diagonal and length \(n-1\) for the sub-diagonal. The rank-one modification matrix \(D+\mathbf{v}\mathbf{v}^{T}\) uses two arrays of length \(n\); one for the main diagonal of \(D\) and another for the entries of vector \(v\), given that the rotation matrix \(G\) uses four arrays of a length equal to the maximum number of rotations that can be applied to \(n\); two arrays with the positions \(i,j\) of the rotations; and two with the values \(s,c\) of the rotations. Dense square matrices use an array of length \(n\times n\) and are stored in _column-major order_. Based on [17; 18], we apply deflation in our implementation if

\[A = \frac{1}{\alpha}B, \tag{11}\]3

We use \(t\) as tolerance because it is an approximation to \(\|D+\mathbf{v}\mathbf{v}^{T}\|\), and \(\epsilon\) is the machine precision. We take it from the macro DBL_EPSILON in file float.h of C language. In the sequential version, we calculate matrix multiplication by using dgemm LAPACK (version 3.4.2). We use dstec LAPACK function to solve sub-problems of the divide-and-conquer algorithm. The main idea of using the divide-and-conquer algorithm is to split matrices into two independent parts. Thus, in our parallel implementation, we send half of the principal diagonal vector and half of the subdiagonal vector to another node with MPI. We solve each of the subproblems generated due to this partitioning by using the dstedx MAGMA (version 1.5.0) routine (version 1.5.0); before doing so, we apply some previously mentioned transformations (scaling and \(\rho\) subtraction) using CUDA and OpenMP. Then, we returned the solution to the principal node to apply deflation and calculate the rank one modification. So now we can combine the two solutions to finally build the complete solution to the original problem. We developed several functions that employ CUDA, OpenMP, and MPI to make these steps. More relevant ones are described next in this piece of writing. There are functions we do not parallelize because we do not improve times; we can take deflation, which is inherently sequential, as an example. Figure 2 shows the different stages of the parallel implementation. Blue sections represent the CPU, and green ones represent the GPU. We execute each phase of the algorithm on CPU, GPU, or both, taking care of the size of the problem. To know when we should change between each, we previously performed several experiments measuring how much time each approach took with different sizes. Another important factor was the suitability of the work to be performed.