Anatomy of High-Performance Column-Pivoted QR Decomposition

In [DG:2017:QR], the authors show several benchmarks of prospective QRCP methods (written in Fortran $90$), compared against both the standard pivoted and unpivoted QR algorithms, provided by the contemporary (though unspecified) version of Intel Math Kernel Library (MKL). The RQRCP method, described in the pseudocode [DG:2017:QR, Algorithm 4] exhibits an order-of-magnitude speedup over the standard pivoted QR and achieves asymptotically up to $60\%$ of the performance of the standard unpivoted QR, when applied to full-rank square matrices. Their work also presents a TRQRCP scheme ([DG:2017:QR, Algorithm 6]), which proves faster than RQRCP in low-rank cases.

In the other work [MOHvdG:2017:QR], the authors present a randomized algorithm called HQRRP that offers an order-of-magnitude speedup over GEQP3 and that achieves up to $80$% of the performance of the standard unpivoted QR, GEQRF. In the benchmarks shown therein, the standard pivoted and unpivoted QR functions were taken from LAPACK version $3.4.0$, and the implementations were linked to BLAS from MKL version $11.2.3$. The performance results were reported for HQRRP that was implemented with libflame version $11104$ [FVE09]. In addition to the libflame-based implementation, the authors presented an LAPACK-compatible version of HQRRP.¹⁰¹⁰10Available at https://github.com/flame/hqrrp The emphasis on the practicality and portability of HQRRP can be seen throughout the [MOHvdG:2017:QR] paper, suggesting that, at its time, this algorithm qualified for an ideal replacement of GEQP3 in LAPACK.

The development of these algorithms, however, was not the final chapter in the quest for a high-performance QRCP method. As with its predecessors, neither HQRRP nor any of the variants of RQRCP made it into LAPACK. Many years of hardware development have passed, and benchmarking on modern machines shows that the gap in performance between the standard pivoted and unpivoted QR factorization algorithms has grown from roughly $10\times$ to near $100\times$. Simultaneously, while the LAPACK-compatible version of HQRRP proved to be easily portable to modern libraries, its performance did not increase commensurately. See Figure 1.

With that, following the steps of our predecessors, we accepted the challenge of formulating a modern approach for QR with column pivoting that would match the contemporary performance of the algorithms using Level 3 BLAS functionality.

Our newest installment in the QRCP family of algorithms is the BQRRP algorithmic framework. In this framework, one implements a BQRRP instantiation by selecting a set of subroutines, most suitable for a system on which BQRRP is to be used. Our main contributions are: (a) a detailed description of the most suitable practical choices for such subroutines on two modern systems; and (b) CPU and GPU implementations of BQRRP, as part of an open-source RandLAPACK library, that allows users to configure the algorithm in a way that is most suitable for their needs. This manuscript describes both CPU and GPU versions of BQRRP. As we show in Section 7, an implementation of BQRRP_CPU not only outperforms HQRRP, but also it is up to two orders of magnitude faster than the standard pivoted QR (GEQP3) implementation available in MKL 2025. BQRRP_GPU exhibits excellent throughput and shows reasonable performance relative to an unpivoted QR (GEQRF) implementation available in NVIDIA’s cuSOLVER 11.6.1.9. This makes BQRRP_GPU a strong candidate for adoption as a standard general QRCP method by GPU linear algebra library vendors.

The dominant cost of BQRRP_CPU comes from a routine that applies an implicitly-stored orthonormal matrix to a portion of the input data in a loop. On a GPU, in addition to this bottleneck, there is the added complication in the form of the cost of permuting the columns of a portion of the input matrix at every iteration of the main loop. This cost can be mitigated by employing advanced pivoting strategies that allow for additional levels of parallelism, such as the “parallel pivots” approach. We discuss these in Section 4. The rest of the major operations in BQRRP are all highly efficient, reaching the Level 3 BLAS performance.

BQRRP_CPU is designed to be storage-efficient, applying all of its subsequent operations in place (modulo comparatively tiny workspace buffers). Furthermore, the output format of both BQRRP_CPU and BQRRP_GPU is identical to that of GEQP3.

Whether our novel approach remains relevant in the long run is an open question. However, even in the worst-case scenario, we emphasize that the value of our work extends beyond delivering a modern high-performance algorithm for QRCP. It also lies in establishing a solid foundation for analyzing the performance and implementation details of future QRCP methods, as well as exposing the underlying structure of such algorithms, which can be invaluable for future developers if the scourge of slow QRCP returns.

Matrices appear in boldface sans-serif capital letters. The transpose of a matrix $\bm{\mathsf{X}}$ is given by $\bm{\mathsf{X}}^{\text{\tiny{T}}}$. Numerical vectors appear as boldface lowercase letters, while index vectors appear as uppercase letters. The identity matrix is denoted by $\bm{\mathsf{I}}$. We enumerate components of matrices and vectors with indices starting from zero, rather than starting from one. We extract the leading $k$ columns of $\bm{\mathsf{X}}$ by writing $\bm{\mathsf{X}}(\hskip 1.25pt{}{:}\hskip 0.24994pt{}{},{0}{:}{k})$. This range excludes column $k$, while its trailing $n-k$ columns are extracted by writing $\bm{\mathsf{X}}(\hskip 1.25pt{}{:}\hskip 0.24994pt{}{},{k}{:}{n})$. The $(i,j)^{\text{th}}$ entry of $\bm{\mathsf{X}}$ is $\bm{\mathsf{X}}(i,j)$. Similar conventions apply to extracting the rows of a matrix or components of a vector.

The matrix we ultimately aim to decompose is denoted by $\bm{\mathsf{M}}$ and has dimensions $m$-by-$n$, with $m$ and $n$ not related in any particular way. We use the letters $i$ and $j$ as zero-based integer indices between 0 and $\min\{m,n\}-1$. In some contexts, $i$ denotes the iteration of an algorithm loop. The symbol $b$ is used to refer to the block size parameter used in a given algorithm; and $b$ is expected to be small relative to the matrix size, i.e., $b\ll\min\{m,n\}$. The symbol $\ell$ denotes the estimated rank of a given matrix. The symbol $k$ denotes the block rank, or the rank of the given submatrix, $k\leq b$.

Given an iteration $i\in{0}{:}{(\lceil n/b\rceil-1)}$ of an arbitrary algorithm that parses the given matrix in increments of block size $b$, we use $\mathit{s}=i\times b$ to denote the start of the current row/column block range; $\mathit{r}=\min\{m,(i+1)b\}$ denotes the (exclusive) end of the current row block range; and $\mathit{c}=\min\{n,(i+1)b\}$ denotes the (exclusive) end of the current column block range.

Throughout this paper, we aim to clearly identify the purpose of each function name, which are presented in typewriter-style font.

We often refer to linear algebraic functions using their conventional BLAS and LAPACK names in uppercase. The explanation of the standard naming conventions can be found at https://www.netlib.org/lapack/lug/node24.html. Whenever we mention a given LAPACK function name, we avoid using the letter that denotes the precision, assuming that all computations are performed in double precision (for example, we use “GEQP3” instead of “DGEQP3”). We use BQRRP (fixed-width font) when referring to any practical algorithm developed from the BQRRP (standard font) algorithmic framework. We use BQRRP_CPU and BQRRP_GPU to specifically denote the respective CPU and GPU versions of BQRRP.

A BQRRP algorithm intends to provide output in a format that is identical to the standard LAPACK QRCP routine, GEQP3. Assuming that a matrix $\bm{\mathsf{M}}\in\mathbb{R}^{m\times n}$ is passed as input into GEQP3, the output format is described as follows:

• •

  Vector $\bm{\tau}$ of length $\min\{m,n\}$ that holds the scalar factors of the elementary reflectors.

• •

  Modified in-place $\bm{\mathsf{M}}$ with its above-diagonal portion occupied by an explicit upper-trapezoidal $\bm{\mathsf{R}}$ of size $\min\{m,n\}\times n$. The below-diagonal portion of output $\bm{\mathsf{M}}$ stores Householder vectors $\bm{v}_{i}$ for $i={0}{:}{\min\{m,n\}}$, which, together with the vector $\bm{\tau}$, compactly represents the orthonormal matrix $\bm{\mathsf{Q}}$ as a product of $\min\{m,n\}$ elementary reflectors:

  	$\displaystyle\bm{\mathsf{Q}}=\bm{\mathsf{H}}_{1}\bm{\mathsf{H}}_{2}\dots\bm{\mathsf{H}}_{\min\{m,n\}}$
  	$\displaystyle\bm{\mathsf{H}}_{i}=\bm{\mathsf{I}}-\tau\bm{v}_{i}\bm{v}_{i}^{\text{\tiny{T}}}.$

• •

  Permutation vector $\bm{J}$ of size $n$ such that if $\bm{J}(j)=i+1$, then the $j^{\mathrm{th}}$ column of $\bm{\mathsf{M}}(:,\bm{J})$ was the $i^{\mathrm{th}}$ column of $\bm{\mathsf{M}}$. Note that in LAPACK convention, the permutation vector $\bm{J}$ stores entries in a one-based index format used by Fortran and MATLAB (hence the presence of “$+1$” term in the $\bm{J}(j)=i+1$ expression).

We refer back to the above output format throughout the paper when talking about the “intended” output format for BQRRP. The described format of the $\bm{\mathsf{Q}}$ factor is referred to as the implicit economical storage format. By contrast, some of the algorithms described in this manuscript use an explicit economical storage format, where $\bm{\mathsf{Q}}$ is comprised of $\min\{m,n\}$ explicitly-defined orthogonal column vectors.

Versions of BQRRP algorithm discussed in this work, as well as their subcomponent functions, are implemented as part of an open-source C++ library called RandLAPACK. Together with its counterpart, RandBLAS, RandLAPACK provides a platform for developing, testing, and benchmarking high-performance RandNLA algorithms. This software was produced as part of the BALLISTIC project [BALLISTIC], and it is actively developed by the authors of this paper. It is important to note that RandBLAS and RandLAPACK rely on BLAS++ and LAPACK++ [gates2022portable] for the purpose of obtaining basic linear algebra functions; BLAS++ and LAPACK++ themselves serve as wrappers around the low-level vendor-optimized BLAS and LAPACK packages. As such, many vendor-optimized libraries (Intel’s MKL, AMD’s AOCL, Apple Accelerate, etc.) can be used to supply our software with basic linear algebra capabilities. Furthermore, RandBLAS relies on the Random123¹¹¹¹11Available at https://github.com/DEShawResearch/random123. library for the purpose of exporting random number generators. For a detailed description of RandBLAS, visit:

https://randblas.readthedocs.io/en/latest/.

In contrast to RandBLAS, RandLAPACK’s scope still continues to change and expand. Because of this, we have yet to define explicit project documentation.

All experiments in this manuscript were run using the following version of our software:

https://github.com/BallisticLA/RandLAPACK/releases/tag/BQRRP-benchmark.

The code for BQRRP can be found in /RandLAPACK/drivers/rl_bqrrp*. The code for constructing and dispatching the experiments can be found in /benchmark/bench_BQRRP/*, as well as /test/drivers/bench_bqrrp*.

An automated script for re-running all of the experiments shown in this paper, as well as the MATLAB plotting scripts for replicating our plots, can be found in:

https://github.com/BallisticLA/BQRRP_benchmarking.

Throughout this paper, we measure the algorithm performance via the canonical FLOP rate. This metric is obtained by dividing the FLOP count of a standard algorithm for a given matrix size by the wall clock time required to run a comparable algorithm under consideration. For example, for comparing canonical FLOP rates of various versions of QR and QRCP factorizations, we use the flop count of the standard LAPACK unpivoted QR called GEQRF [LAWN41:1994, Page 121] and divide it by the wall clock time of any given QR or QRCP algorithm being compared. The canonical FLOP rate gives a clear way to compare algorithm performance by using a standard FLOP count divided by the algorithm’s actual runtime. Unlike the raw FLOP rate, which gives credit to unnecessary work, the canonical FLOP rate sets a fair comparison by holding all algorithms to the same amount of work. While runtime alone could be used for comparisons, the FLOP rates are useful because they connect to familiar performance standards, such as the peak FLOP rates advertised for the hardware or standard algorithms like GEMM for matrix-matrix multiply. FLOP rates also help reveal scalability issues across problem sizes: if an algorithm’s FLOP rate drops a lot as the problem size grows, it may reveal an inefficient method or implementation. Canonical FLOP rates, therefore, give a balanced measure that considers both runtime and standard performance benchmarks.

All test matrices were generated using RandBLAS. For the performance experiments, each matrix entry is independently sampled from the standard normal distribution. However, it should be noted that such matrices are appropriate for our performance tests as being representative of “typical” user inputs, and they have predictable numerical and performance properties. For the numerically challenging matrices, see Section 6. The average-case performance, on which we focus extensively here, is much better observed with our primary choice of matrix elements drawn from i.i.d. distribution, giving a representative behavior in pivot column interchanges.

We focus our experiments on large matrices where the performance gap between Level-2 and Level-3 BLAS based algorithms is readily apparent. Simultaneously, we want to ensure that the experiments we set up terminate in reasonable time. Hence, we do not conduct runs on input matrices of the absolute largest sizes that can fit on our systems. Once we finish analyzing the established versions of the BQRRP algorithm for the larger matrix sizes, we present the investigation of the performance of various methods when applied to input matrices of a wide range of sizes, using various numbers of OpenMP threads, in Section 7.

There are some known performance concerns when it comes to selecting the specific sizes of the input matrices. Specifically, matrices with column sizes being powers of 2 experience increased demand for the main memory bandwidth based on how the modern cache memory works: elements from nearby columns are mapped to the exact cache line because power-of-2 memory addresses have all lower bits set to 0. Note that all modern libraries use a packed storage format for an efficient in-cache matrix multiply kernel, but they still have to read and write the cached data from/to the main memory, causing false sharing of cache lines. Due to a much-simplified structure of caches, GPUs do not experience this problem, and thus, we needed to compromise on the matrix size selection for comparable results between the types of computing devices. Because of this, we chose to run our experiments on matrices with dimensions that are powers of two and multiples of ten. Additionally, in our main performance experiments with various QR and QRCP algorithms, we focus on square input matrices. This choice provides a meaningful baseline for performance comparisons, while ensuring a balanced computational workload and even distribution of work across threads. Nonetheless, we do present the performance results captured on both tall and wide matrices in Section C.2.

All of our tests used double-precision arithmetic, all code was compiled with the -O3 flag. The performance of each algorithm is determined by selecting its best execution time from five consecutive runs. Instead of running each algorithm individually five times, we execute the entire set of compared algorithms in sequence, repeating this set five times to capture comparatively consistent performance.

The detailed configurations of the machines that we used for conducting the CPU experiments are shown in Table 1.

When running benchmarks on CPUs described in Table 1, we set the number of OpenMP threads (with OMP_NUM_THREADS environment variable) to the maximum value of threads available, unless specified otherwise. While using the maximum number of available threads may not always be the optimal approach to achieve peak performance, we advocate for harnessing all available computational resources. Both tested Intel and AMD systems featured dual-socket configurations. We launched tests with the numactl --interleave=all command when running our experiments to balance the memory usage across the two memory controllers.

The system setup used for the GPU experiments can be found in Table 2. Since the GPU listed in Table 2 has “only” $80$ GB of memory, we are unable to run double-precision experiments with the input matrices of certain sizes on this GPU.

Step 6 in Algorithm 2 is crucial, since the format in which a pivoted LU represents the permutation vector $\bm{J}$ is different from the pivoted QR format, and hence the format conversion procedure is required. In the context of pivoted LU factorization, row $i$ of the input matrix was interchanged with row $\bm{J}^{\mathrm{lu}}(i)$. The format conversion consists of first creating a vector $\bm{J}^{\mathrm{qr}}$ of length $n$ with entries from $1$ to $n$ and then serially swapping elements in it according to the entries in $\bm{J}^{\mathrm{lu}}$. Simply put, for element at index $i\in{0}{:}{(n-1)}$, element at $\bm{J}^{\mathrm{qr}}(i)$ is to swap positions with the element at $\bm{J}^{\mathrm{qr}}(\bm{J}^{\mathrm{lu}}(i)-1)$.

When Algorithm 2 is in use in BQRRP, the first $b$ components of $\bm{J}^{\text{qr}}$ are the same for $\gamma=1$ and $\gamma>1$. Therefore the rank-revealing properties of BQRRP using Algorithm 2 for qrcp_wide cannot be improved by using $\gamma>1$.

Observe a practical comparison of the performance of the two approaches in Figure 2. The performance superiority of Algorithm 2 to the standard QRCP approach makes it the best option to be used in a BQRRP implementation.

LAPACK offers several algorithms suitable for tall QR. The natural choice in this context is LATSQR [lapack_latsqr]: a specialized Householder QR factorization for tall matrices. Alternatively, one could use a standard Householder QR factorization, GEQRF [lapack_geqrf]. As the third option, GEQRT [lapack_geqrt], is a blocked algorithm based on compact WY representation [Bischof:1987:WRP:37316.37324]. Finally, LAPACK’s GEQR calls either LATSQR or GEQRT, depending on the size of the input matrix and the output from LAPACK’s ILAENV inquiry function [lapack_ilaenv]. In addition to listing the standard LAPACK algorithms, we consider the use of Cholesky QR in this context, following our approach from prior work [MBM2024]. Background on Cholesky QR is given in Appendix B.

Although Cholesky QR can be safe to use in the context of BQRRP, its use produces the following complication: steps 24 and 27 require access to the fully-formed $m$-by-$m$ representation of the orthonormal factor computed at the current iteration, while Cholesky QR only outputs an explicit economical representation $\bm{\mathsf{Q}}^{\mathrm{chol}}\in\mathbb{R}^{m\times k}$. We can acquire the full representation by using a specialized LAPACK routine ORHR_COL. This routine transforms an explicit economical $\bm{\mathsf{Q}}^{\mathrm{chol}}$ output by Cholesky QR into an implicit representation via Householder reconstruction. The description of the basic approach for performing the Householder reconstruction can be found in [BD2015, Algorithms 5, 6]. An advanced recursive implementation of ORHR_COL is described in its respective Netlib LAPACK documentation [lapack_orhr_col]; see also [F1997].

Algorithm 3 shows a practical implementation of Cholesky QR and its dependencies, offering a method that can be used in step 20.

Figure 3 illustrates how the use of preconditioning and Householder restoration affect the performance of Cholesky QR, relative to alternative methods for tall QR factorization. Figure 3 shows that GEQRF is the best alternative method for tall QR factorization across both systems. Additionally, favoring the use of Cholesky QR on the Intel system for the input matrices of select sizes can be a reasonable option.

Only GEQRF outputs the representation of $\bm{\mathsf{Q}}$ that precisely matches the one we are after (described in Section 1.3). The format of $\bm{\mathsf{Q}}$ produced by ORHR_COL differs from the intended one in the following way: instead of generating a vector with scalar factors of the elementary reflectors $\bm{\tau}$, it produces an upper-trapezoidal matrix $\bm{\mathsf{T}}$ that represents a sequence of upper triangular block reflectors stored compactly. Each block in $\bm{\mathsf{T}}$ is of size $n_{b}\times b$ (except possibly in the last iteration if $n$ is not evenly divisible by $b$), with $n_{b}$ being the block size parameter. The scalar factors of the elementary reflectors, originally intended for $\bm{\tau}$, are stored along the diagonals of the block matrix $\bm{\mathsf{T}}$.

For details of the representations used by GEQR, LATSQR, and GEQRT, we refer the readers to the official LAPACK documentation.

There are two main methods for applying $(\bm{\mathsf{Q}}^{\mathrm{curr}})^{\text{\tiny{T}}}$ to a given matrix; the choice of the right method relates to the decision made in the tall QR step.

The first available method is LAPACK’s ORMQR function [lapack_ormqr], which takes reflectors produced by GEQRF or GEQP3. Recall that Section 2.3 states that the GEQP3-compatible representation of $\bm{\mathsf{Q}}$ must be obtained regardless of the algorithm used in the tall QR step. Therefore ORMQR is always an option after obtaining the GEQP3-compatible representation. The cost of ORMQR is $4nmk-2nk^{2}+3nk$ FLOPs [LAWN41:1994, Page 122].

The second option to use in this context is the GEMQRT function [lapack_gemqrt]. The input format of GEMQRT matches the output format of ORHR_COL. Since Cholesky QR must use ORHR_COL, this implies pairing Cholesky QR with GEMQRT can be an efficient choice. Alternatively, one could pair GEMQRT with any tall QR method, first ensuring that the representation of the output $\bm{\mathsf{Q}}$-factor is made compatible with the input representation of GEMQRT.

One more alternative to the methods described above is avoiding the trailing update altogether, as described in [DG:2017:QR, Algorithm 5.1]. This approach is, however, only applicable to low-rank data, and in an implementation, it would further increase the complexity of a rather nontrivial BQRRP scheme.

The relative performance of the alternative ways of applying an orthonormal factor $(\bm{\mathsf{Q}}^{\mathrm{curr}})^{\text{\tiny{T}}}$ to an arbitrary matrix of size $m$-by-$(n-b)$ is shown in Figure 4. The performance superiority of ORMQR function to the alternatives in the majority of the explored cases makes it the best option to be used in a BQRRP implementation.

At step 12, performing QRCP on the sketch as described in Algorithm 2 requires an $n\times d$ buffer for transposing ${\bm{\mathsf{M}}}^{\mathrm{sk}}$ (the full buffer size is only needed at $i=0$). The important data computed in this step at the iteration $i$ is stored as follows: the upper-trapezoidal ${\bm{\mathsf{R}}}^{\mathrm{sk}}\in\mathbb{R}^{d\times(n-\mathit{s})}$ is stored at ${\bm{\mathsf{M}}}^{\mathrm{sk}}(:,{\mathit{s}}{:})$; and a vector ${\bm{J}}^{\mathrm{sk}}\in\mathbb{R}^{(n-\mathit{s})}$ is stored in a buffer of length $n$. In our further elaboration, we partition ${\bm{\mathsf{R}}}^{\mathrm{sk}}$ into three components: upper-triangular ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{11}={\bm{\mathsf{M}}}^{\mathrm{sk}}({0}{:}{b},{\mathit{s}}{:}\mathit{c})$; rectangular ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{12}={\bm{\mathsf{M}}}^{\mathrm{sk}}({0}{:}{b},{\mathit{c}}{:})$; and an upper-trapezoidal ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{22}={\bm{\mathsf{M}}}^{\mathrm{sk}}({b}{:},{\mathit{c}}{:})$.

Our implementation of BQRRP allows for alternatively performing QRCP via GEQP3, in which case a buffer for transposing ${\bm{\mathsf{M}}}^{\mathrm{sk}}$ is not needed.

In step 14, if the block rank $k$ (computed as described in Section 2.2) is not equal to $\min\{b,(n-\mathit{s}),(m-\mathit{s})\}$ (when the current block is not full rank), then the given iteration would be the final one, as the termination criteria at step 32 states. In that case, the trailing portion of the matrix $\bm{\mathsf{M}}$ is not updated at step 29.

If BQRRP is set to use Cholesky QR in step 20, then a rank estimate is employed to ensure no infinite or not-a-number values appear when the preconditioning is performed.

Since the intended output format of BQRRP matches that of GEQP3, the output $\bm{\mathsf{R}}$ factor will be stored in the upper-triangular portion of $\bm{\mathsf{M}}$’s memory space. At iteration $i\geq 1$, step 16 permutes the trailing columns in the computed rectangular portion of the $\bm{\mathsf{R}}$-factor, implying that $\bm{\mathsf{M}}({0}{:}{\mathit{s}},{\mathit{s}}{:})$ is to be permuted in accordance with ${\bm{J}}^{\mathrm{sk}}$. Meanwhile, step 18 requires $\bm{\mathsf{M}}({\mathit{s}}{:},{\mathit{s}}{:})$ to be permuted by the same vector. Therefore, steps 16 and 18 can be combined into permuting $\bm{\mathsf{M}}(\hskip 1.25pt{}{:}\hskip 0.24994pt{}{},{\mathit{s}}{:})$ at every iteration of the BQRRP’s main loop.

After the columns of $\bm{\mathsf{M}}(\hskip 1.25pt{}{:}\hskip 0.24994pt{}{},{\mathit{s}}{:})$ have been permuted, we may perform an early termination check by verifying whether the first column in $\bm{\mathsf{M}}(\hskip 1.25pt{}{:}\hskip 0.24994pt{}{},{\mathit{s}}{:})$ consists of all zeros. In that case, BQRRP terminates immediately. If at iteration $i$ we have $k\neq\min\{b,n-\mathit{s},m-\mathit{s}\}$, then BQRRP should avoid permuting $\bm{\mathsf{M}}({\mathit{s}+k}{:},{\mathit{c}}{:})$, since this trailing portion of the input matrix would not be used (i.e., the trailing update to $\bm{\mathsf{M}}$ would be avoided in step 24). This was not shown in Algorithm 1 for simplicity of presentation.

By step 19, the pivot vector ${\bm{J}}^{\mathrm{sk}}$, computed at the current iteration, has been used for all the necessary internal operations. It can now be incorporated into the vector $\bm{J}$ that is to be output by BQRRP. When $i$ is $0$, this is done by copying ${\bm{J}}^{\mathrm{sk}}$ into $\bm{J}$; at any subsequent iteration, the trailing portion of $\bm{J}$, $\bm{J}({\mathit{s}}{:})$, will be permuted in accordance with ${\bm{J}}^{\mathrm{sk}}$ via (a vector version of) Algorithm 4. Since ${\bm{J}}^{\mathrm{sk}}$ is not needed after this step, we do not have to create a copy of it in Algorithm 4.

When using GEQRF in step 20, we perform the factorization on all columns of ${\bm{\mathsf{M}}}^{\mathrm{perm}}=\bm{\mathsf{M}}({\mathit{s}}{:},{\mathit{s}}{:}\mathit{c})$, regardless of whether $k$, the current column block rank, is equal to $k_{\mathrm{max}}=\min\{b,(n-\mathit{s}),(m-\mathit{s})\}$. This is because we want $\bm{\mathsf{R}}_{11}$ to have $(\mathit{c}-\mathit{s})=\min\{b,(n-\mathit{s})\}$ columns. Despite the fact that $(\mathit{c}-\mathit{s})$ reflectors are computed in that case, we use only $k$ of them outside of this step. Using GEQRF on ${\bm{\mathsf{M}}}^{\mathrm{perm}}$ in this step requires no additional storage. This function produces $\bm{\mathsf{Q}}^{\mathrm{curr}}$, represented by $(\mathit{c}-\mathit{s})$ reflectors of length $(m-\mathit{s})$ (stored in the portion of $\bm{\mathsf{M}}({\mathit{s}}{:},{\mathit{s}}{:}\mathit{c})$ below the diagonal) and the vector of scalar factors of the elementary reflectors (stored in $\bm{\tau}({\mathit{s}}{:}\mathit{c})$). The explicit $\bm{\mathsf{R}}_{11}$ of size $(\mathit{r}-\mathit{s})\times(\mathit{c}-\mathit{s})$ is stored right above the $\bm{\mathsf{Q}}^{\mathrm{curr}}$, in the upper-triangular portion of $\bm{\mathsf{M}}({\mathit{s}}{:}\mathit{r},{\mathit{s}}{:}\mathit{c})$. Note that in this case, step 22 is performed implicitly, since $\bm{\mathsf{R}}_{11}$ will be stored where it needs to be on output from GEQRF.

Using the Cholesky QR approach in step 20 is comprised of four parts: preconditioning; performing Cholesky QR on the preconditioned matrix; implicit $\bm{\mathsf{Q}}^{\mathrm{curr}}$ reconstruction; and computing $\bm{\mathsf{R}}_{11}$. All of these are done in accordance with Algorithm 3.

The preconditioning is to be done on a portion of the matrix ${\bm{\mathsf{M}}}^{\mathrm{perm}}$. We perform ${\bm{\mathsf{M}}}^{\mathrm{pre}}={\bm{\mathsf{M}}}^{\mathrm{perm}}(\hskip 1.25pt{}{:}\hskip 0.24994pt{}{},{0}{:}{k})({\bm{\mathsf{R}}}^{\mathrm{sk}}_{11}({0}{:}{k},{0}{:}{k}))^{-1}$, using ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{11}$ stored in ${\bm{\mathsf{M}}}^{\mathrm{sk}}({0}{:}{b},{\mathit{s}}{:}\mathit{c})$. This step requires no workspace buffers, and the matrix ${\bm{\mathsf{M}}}^{\mathrm{pre}}\in\mathbb{R}^{(m-\mathit{s})\times k}$ is stored at $\bm{\mathsf{M}}({\mathit{s}}{:},{\mathit{s}}{:}(\mathit{s}+k))$ after the preconditioning. The next step is to perform Cholesky QR on ${\bm{\mathsf{M}}}^{\mathrm{pre}}$.

As explained in Section 2.3, Cholesky QR is comprised of: computing the Gram matrix (SYRK), performing the Cholesky factorization (POTRF), and forming the explicit economical version of the orthonormal factor (TRSM). These functions require (at most) a single $b\times b$ workspace buffer to store the Gram matrix and the output upper-triangular $\bm{\mathsf{R}}^{\mathrm{chol}}$ factor, both of size $k\times k$, $k\leq\min\{b,(n-\mathit{s}),(m-\mathit{s})\}$. The explicit orthonormal factor $\bm{\mathsf{Q}}^{\mathrm{chol}}\in\mathbb{R}^{(m-\mathit{s})\times k}$ is stored at $\bm{\mathsf{M}}({\mathit{s}}{:},{\mathit{s}}{:}(\mathit{s}+k))$, in place of ${\bm{\mathsf{M}}}^{\mathrm{pre}}$.

ORHR_COL used for the Householder reconstruction returns an implicit representation of $\bm{\mathsf{Q}}^{\mathrm{curr}}$ in the form of $k$ reflectors and an upper-trapezoidal matrix $\bm{\mathsf{T}}$ of size $n_{b}\times k$. The $k$ reflectors fit exactly in place of the portion of $\bm{\mathsf{Q}}^{\mathrm{chol}}$ below the diagonal of $\bm{\mathsf{M}}({\mathit{s}}{:},{\mathit{s}}{:}(\mathit{s}+k))$, while $\bm{\mathsf{T}}$ requires an additional buffer of size $n_{b}\times b$. Furthermore, ORHR_COL produces a sign vector $\bm{D}$ of length $k$ that requires a storage buffer of length $b$.

Computing $\bm{\mathsf{R}}_{11}$ requires first applying the entries from the sign vector $\bm{D}$ to the columns of $\bm{\mathsf{R}}^{\mathrm{chol}}$ (this is done so that $\bm{\mathsf{R}}^{\mathrm{chol}}$ is in line with $\bm{\mathsf{Q}}^{\mathrm{curr}}$, computed by ORHR_COL). The next step amounts to undoing the preconditioning on $\bm{\mathsf{R}}^{\mathrm{chol}}$ via multiplying it by an upper-triangular ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{11}({0}{:}{k},{0}{:}{(\mathit{c}-\mathit{s})})$, stored at ${\bm{\mathsf{M}}}^{\mathrm{sk}}({0}{:}{k},{\mathit{s}}{:}\mathit{c})$. The product is temporarily stored in place of $\bm{\mathsf{R}}^{\mathrm{chol}}$ and then is copied into the upper-triangular portion of $\bm{\mathsf{M}}({\mathit{s}}{:}(\mathit{s}+k),{\mathit{s}}{:}\mathit{c})$, placing it right above the implicitly-stored Householder reflectors (step 22).

For the output format of BQRRP to match that described in Section 1.3, we would need to copy the entries from the block diagonals of $\bm{\mathsf{T}}$ into $\bm{\tau}({\mathit{s}}{:}(\mathit{s}+k))$.

The termination criteria check (step 32) simply amounts to verifying whether the maximum number of iterations has been exceeded ($i$ reached $\lceil n/b\rceil-1$) or whether the block rank $k$ does not match the column block size $\min\{b,n-\mathit{s},m-\mathit{s}\}$. Before the termination, the rank parameter $\ell$ is updated to $\mathit{s}+k$ (step 33).

If the maximum number of iterations has not been reached at this point, the sketch ${\bm{\mathsf{M}}}^{\mathrm{sk}}$ is updated at step 35. This step first computes the expression ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{11}(\bm{\mathsf{R}}_{11})^{-1}$, where ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{11}$ is stored in the upper-triangular portion of ${\bm{\mathsf{M}}}^{\mathrm{sk}}({0}{:}{b},{\mathit{s}}{:}\mathit{c})$ and $\bm{\mathsf{R}}_{11}$ is stored in the upper-triangular portion of $\bm{\mathsf{M}}({\mathit{s}}{:}\mathit{r},{\mathit{s}}{:}\mathit{c})$. The result is written into the space of ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{11}$, at ${\bm{\mathsf{M}}}^{\mathrm{sk}}({0}{:}{b},{\mathit{s}}{:}\mathit{c})$; we explicitly zero out the entries in this space below the diagonal. Next, the full expression ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{12}-{\bm{\mathsf{R}}}^{\mathrm{sk}}_{11}(\bm{\mathsf{R}}_{11})^{-1}\bm{\mathsf{R}}_{12}$ is computed. The matrix ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{12}$ is stored in ${\bm{\mathsf{M}}}^{\mathrm{sk}}({0}{:}{b},{\mathit{c}}{:})$ and $\bm{\mathsf{R}}_{12}$ is stored at $\bm{\mathsf{M}}({\mathit{s}}{:}\mathit{r},{\mathit{c}}{:})$. The result of this expression is placed into the space of ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{12}$, at ${\bm{\mathsf{M}}}^{\mathrm{sk}}({0}{:}{b},{\mathit{c}}{:})$. The upper-trapezoidal matrix ${\bm{\mathsf{R}}}^{\mathrm{sk}}_{22}$ is already stored at its intended location, in ${\bm{\mathsf{M}}}^{\mathrm{sk}}({b}{:}(d-b),{\mathit{c}}{:})$, but we need to make sure that the entries are zeroed out below the diagonal. After the sketch update is completed, the “working submatrix” of ${\bm{\mathsf{M}}}^{\mathrm{sk}}$ is implicitly located at ${\bm{\mathsf{M}}}^{\mathrm{sk}}(\hskip 1.25pt{}{:}\hskip 0.24994pt{}{},{\mathit{c}}{:})$.

Recall that Algorithm 4 illustrated a sequential approach to permuting columns. Despite this approach being sequential, we do not anticipate it having a major negative effect on the overall performance of BQRRP_CPU, as described in Section 2.5. By contrast, on a GPU, there are very few hardware resources devoted to dealing with inherently sequential instruction streams or even tolerating the high latency of the main memory transactions. The initially presented sequential loop for pivoting is anathema to the parallel GPU hardware.

The widely-used alternative solution is the “parallel pivoting” strategy, shown in Algorithm 5. By contrast to Algorithm 4, Algorithm 5 does not swap columns within the memory space of a single matrix. Instead, it copies columns from one space into another at the new location from the permutation vector, which allows for parallelizing the main loop, significantly increasing the performance of this column permutation approach. The performance gain offered by Algorithm 5 comes at the cost of increased space usage, as it requires a copy of the input $\bm{\mathsf{M}}$.

In order to minimize the number of explicit copies performed as part of using Algorithm 5 in practice, one could swap the pointers that point to the memory spaces of $\bm{\mathsf{M}}$ and $\bm{\mathsf{M}}^{\mathrm{cpy}}$ before performing column permutation. This approach is safe as long as either the entire matrix $\bm{\mathsf{M}}$ is permuted in Algorithm 5, or the non-permuted part of $\bm{\mathsf{M}}$ is not used outside of Algorithm 5. As stated in Section 2.5, the column permutation is used in steps 16, 18, and 19 in Algorithm 1 (where steps 16 and 18 can be merged) and in step 7 of Algorithm 2. As such, we will be applying permutations to portions of matrices $\bm{\mathsf{M}}^{\mathrm{sk}}$ and $\bm{\mathsf{M}}$, as well as a vector ${\bm{J}}^{\mathrm{sk}}$ at every iteration of BQRRP_GPU.

Since the full $\bm{\mathsf{M}}^{\mathrm{sk}}$ is not used outside of BQRRP_GPU, the pointer-swapping strategy can be used without any memory safety concerns. This is, however, not the case with $\bm{\mathsf{M}}$ and ${\bm{J}}^{\mathrm{sk}}$, as both of their entire memory spaces will have to store correct data on output from the BQRRP_GPU. To illustrate the complication with the pointer-swapping strategy, suppose the main loop in BQRRP_GPU is executing starting with $i=0$, where index $0$ is considered “even.” Then, at the end of the $0^{th}$ iteration, the space of $\bm{\mathsf{M}}^{\mathrm{cpy}}$, will contain the correct data to be a part of the output of BQRRP_GPU since the pointer swapping took place. At iteration $1$, the space $\bm{\mathsf{M}}^{\mathrm{cpy}}$ will contain the first $b$ properly-computed columns, the rest of the correct entries will be contained in $\bm{\mathsf{M}}(:,{b}{:})$. Continuing with that logic, we conclude that the space $\bm{\mathsf{M}}^{\mathrm{cpy}}$ will contain the “correct” entries in $\bm{\mathsf{M}}^{\mathrm{cpy}}(:,{ib}{:}(i+1)b)$ for all even $i\leq n/b$, and the space $\bm{\mathsf{M}}$ will contain the correct entries for all odd $i$ in $\bm{\mathsf{M}}(:,{ib}{:}(i+1)b)$ when BQRRP_GPU terminates. If BQRRP_GPU terminated at an even iteration, the pointers to $\bm{\mathsf{M}}$ and $\bm{\mathsf{M}}^{\mathrm{cpy}}$ will need to be swapped back around. Regardless of the final parity status of the loop’s index, the “correct” columns will need to be copied from $\bm{\mathsf{M}}^{\mathrm{cpy}}$ to $\bm{\mathsf{M}}$. Note that if BQRRP_GPU was to terminate early at an even iteration, the trailing entries in range $(i+1)b$ to $m$ in $\bm{\mathsf{M}}^{\mathrm{cpy}}$ will need to be copied into $\bm{\mathsf{M}}$.

As per the vector $\bm{J}$, since it is not permuted at iteration $0$, the space of $\bm{J}^{\mathrm{cpy}}$ would contain the “correct” entries in $\bm{J}^{\mathrm{cpy}}({ib}{:}(i+1)b)$ for all odd $i\leq n/b$, while the space $\bm{J}$ will contain the correct entries for all even $i$ in $\bm{J}({ib}{:}(i+1)b)$ when BQRRP_GPU terminates. Otherwise, the vector $\bm{J}$ is to be processed similarly to the matrix $\bm{\mathsf{M}}$.