QR Factorization#

GEQRF (GEneral QR Factorization) function computes batched QR factorization of a general matrix:

\[A = Q R\]

using Householder reflection transformations, where:

  • \(A\) is a batched \(M \times N\) matrix, with leading dimension \(\mathrm{lda} \geq M\) if matrix A is in column-major layout, or \(\mathrm{lda} \geq N\) if matrix A is row-major,

  • \(Q\) is a unitary \(M \times M\) matrix, and

  • \(R\) is an upper triangular matrix if \(M \geq N\), or upper trapezoidal matrix if \(M < N\).

cuSolverDx geqrf device functions are (see Execution Methods):

__device__ void execute(data_type* A, data_type* tau);
// with the runtime lda
__device__ void execute(data_type* A, const unsigned int lda, data_type* tau);

After the function, the upper triangular or upper trapezoidal part of input A, including diagonal elements, is replaced by the matrix R.

Matrix \(Q\) is not explicitly formed. The elements below the diagonal of A, with the array tau, represent the matrix \(Q\) as a product of \(\min(M, N)\) Householder vectors:

\[Q = H(0) * H(1) * . . . * H(\min(M, N) - 1).\]

Each Householder vector has the form \(H(i) = I - \mathrm{tau}[i] v v^H\), where:

  • tau is an array of size \(\min(M, N)\) for each batch, and

  • v is a vector of size M for each batch, with v[0 : i-1] = 0, v[i] = 1, and v[i+1 : M] is stored on exit in A[i+1 : M, i].

The functions support A being either column- or row-major memory layout, see Arrangement operator.

Warning

Unlike LAPACK’s geqrf, cuSolverDx geqrf functions do not guarantee the diagonal of R is always real valued for complex inputs and applies reflectors in some places they could be skipped.

Cluster execution (experimental support starting from version 0.6.0)#

Starting from cuSolverDx 0.6.0, Cluster operator provides an experimental GEQRF path for matrices that do not fit in a single thread block’s shared memory. A thread block cluster stores a tile of A in distributed shared memory (DSMEM) and runs a binary-tree reduction across BlocksPerCluster (BPC) cooperating CTAs, then reconstructs compact Householder vectors matching the standard GEQRF output layout (Algorithm 9 of Ballard et al. [1]).

Row tiling the global matrix

The tall M x N matrix is partitioned into BPC tiles. Each CTA owns one tile of size TileSize x N, where:

\[\mathrm{TileSize} = \mathrm{ceil}(M / BPC)\]

To be able to use the cluster execution:

  • the matrix size must satisfy M >= N * BPC;

  • BlocksPerCluster must be a power of two; and

  • required shared memory size must not exceed the total shared memory size available on the selected GPU.

When BPC = 1, cluster GEQRF falls back to the ordinary block GEQRF path.

Each CTA allocates Solver::shared_memory_size bytes for its local TileSize x N tile, the length-N tau buffer, and the internal workspace (Solver::workspace_size elements). Remote tiles are accessed through DSMEM via cluster.map_shared_rank.

Example ( M = 1023, N = 32, BPC = 4 ):

NB = ceil(1023 / 4) = 256

block 0: rows [0, 256)
block 1: rows [256, 512)
block 2: rows [512, 768)
block 3: rows [768, 1023)

Load/store helpers

Before calling the cluster execution device function, load each CTA’s row tile from global memory into DSMEM, then store the factored tiles back afterward. Use the description member functions:

  • Solver::load_global_to_dsmem(A, lda, a_local) — each cluster block loads its TileSize x N row tile from global memory into local shared memory

  • Solver::store_dsmem_to_global(A, lda, a_local) — the reverse store after execute

Note

Cluster execute expects the tile workspace in shared memory to use exactly the same distribution and layout as produced by load_global_to_dsmem (and read back by store_dsmem_to_global): same TileSize = ceil(M / BPC) row ownership and arrangement. If you load or arrange tiles differently, execute will not factor the matrix correctly.

Important

Output tau array of the entire matrix is written by block 0 into its local shared memory tau buffer during the execute function. Only copy the local tau shared memory buffer in the first thread block to global memory after the cluster execution function returns.

Device function API

Cluster execute function signature is:

__device__ void execute(data_type* A, data_type* tau, data_type* workspace);

  • A must point to the local TileSize x N tile in shared memory (not the full global matrix). There is no runtime lda argument.

  • tau is a length-N buffer in shared memory.

  • workspace is a temporary buffer of Solver::workspace_size elements in shared memory (included in Solver::shared_memory_size).

TileSize is fixed for GEQRF cluster: it must equal TileSize = ceil(M / BPC) (also exposed as Solver::suggested_tile_size). See geqrf_batched_cluster for a complete kernel.

References#