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:

  1. A being either column- or row-major memory layout, see Arrangement operator, and

  2. A being in lower or upper fill mode, see FillMode operator.

Cluster execution (experimental support in 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, one has to load the matrix A from the global memory into the DSMEM tile workspace, then store the factored tiles back afterward. The load/store help functions can be found in potrf_cluster_io.cuh:

  • load_global_to_cluster_tiles — each owner block loads its tiles from global memory (via matrix_tile) into local shared memory

  • store_cluster_tiles_to_global — 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_cluster_tiles (and read back by store_cluster_tiles_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 signatures match Thread Operator and Block Operator, but A must point to the cluster tile workspace in shared memory, not the original global matrix, and the runtime leading dimension lda is ignored. See potrf_batched_cluster and potrf_cluster_io.cuh for a complete kernel.