Cholesky Factorization#

POTRF (POsitive definite TRiangular Factorization) function computes batched Cholesky factorization of a Hermitian positive-definite matrix:

\[ \begin{align}\begin{aligned}A = L L^H,\\A = U^H U\end{aligned}\end{align} \]

\(A\) is a batched \(M \times M\) Hermitian matrix (with leading dimension \(\mathrm{lda} \geq M\)), only the lower or upper part is meaningful. The input FillMode operator indicates which part of the matrix \(A\) is used.

cuSolverDx POTRF device functions (see Execution Methods):

__device__ void execute(data_type* A, status_type* info);
// with runtime leading dimension of the matrix A
__device__ void execute(data_type* A, const unsigned int lda, status_type* info);

For lower fill mode, only the lower triangular part of A is processed, and replaced by the lower triangular Cholesky factor L. The upper part of the matrix is untouched.

For upper fill mode, only upper triangular part of A is processed, and replaced by upper triangular Cholesky factor U. The lower part of the matrix is untouched.

If Cholesky factorization failed for any batches, i.e. some leading minor of A is not positive definite, or equivalently some diagonal elements of L or U is not a real number, the corresponding output status parameter, info[batch_id] would indicate smallest leading minor of A which is not positive definite.

The functions support:

A being either column- or row-major memory layout, see Arrangement operator, and
A being in lower or upper fill mode, see FillMode operator.

Cluster execution (experimental support starting from version 0.5.0)#

Starting from cuSolverDx 0.5.0, Cluster operator provides an experimental blocked POTRF path for matrices that do not fit in a single thread block’s shared memory. A thread block cluster stores the triangular factor tiles in distributed shared memory (DSMEM) and runs a blocked right-looking Cholesky factorization across BlocksPerCluster (BPC) cooperating CTAs.

Tiling the global matrix

The stored triangle of A (lower or upper, per FillMode) is partitioned into a grid of TileSize x TileSize tiles. For matrix size N and tile size TileSize, there are ceil(N / TileSize) tile rows and columns; only tiles inside the chosen triangle are kept. Each such tile is numbered based on its row and column index in the triangle, and assigned to cluster blocks round-robin by the tile number.

Each CTA allocates Solver::shared_memory_size bytes for its owned tiles (ceil(num_triangle_tiles / BPC) tiles of size TileSize x TileSize). Remote tiles are accessed through DSMEM via cluster.map_shared_rank.

Example (upper fill mode, n_tiles = 4, BPC = 4 ):

Global tile layout (X marks unused lower triangle):

| A00  A01  A02  A03 |
|  X   A11  A12  A13 |
|  X    X   A22  A23 |
|  X    X    X   A33 |

Round-robin owner (tile index → block):

A00 (0) → block 0    A01 (1) → block 1    A11 (2) → block 2    A02 (3) → block 3
A12 (4) → block 0    A22 (5) → block 1    A03 (6) → block 2    A13 (7) → block 3
A23 (8) → block 0    A33 (9) → block 1

Tiles owned by each block (local slot order in that block’s shared memory):

block 0: A00, A12, A23
block 1: A01, A22, A33
block 2: A11, A03
block 3: A02, A13

Load/store helpers

Before calling the cluster execution device function, load the matrix A from global memory into the DSMEM tile workspace, then store the factored tiles back afterward. Use the description member functions:

Solver::load_global_to_dsmem(A, lda, cluster_tiles) — each owner block loads its tiles from global memory into the cluster tile workspace
Solver::store_dsmem_to_global(A, lda, cluster_tiles) — the reverse store after execute

Note

Cluster execute expects the tile workspace in shared memory to use exactly the same tile distribution and layout as produced by load_global_to_dsmem (and read back by store_dsmem_to_global): same round-robin owner mapping, local slot order, and TileSize x TileSize packing per tile. If you load or arrange tiles differently, execute will not factor the matrix correctly.

Device function API

Cluster execute function signatures are:

__device__ void execute(data_type* A, status_type* info);

A must point to the cluster tile workspace in shared memory, not the original global matrix. The tile workspace is arranged by TileSize, not by a leading dimension, therefore there is no runtime lda argument.

TileSize must be between 1 and N. If omitted, the default is Solver::suggested_tile_size. When BlocksPerCluster is 1, TileSize must equal N. See potrf_batched_cluster example for a complete kernel.