 |
Eigen-Contrib
5.0.1
|
|
|
Eigen-Contrib
5.0.1
|
|
GPU-accelerated linear algebra for Eigen users, dispatching to NVIDIA CUDA Math Libraries (cuBLAS, cuSOLVER, cuFFT, cuSPARSE, cuDSS). Requires CUDA 11.4+; cuDSS features require CUDA 12.0+ and a separate cuDSS install. Header-only.
This module dispatches rather than reimplements, so numerical behavior, supported shapes and scalar types, and performance characteristics are the vendor libraries'. Their documentation is the reference for anything this file does not state:
| Library | Used for | Documentation |
|---|---|---|
| CUDA Toolkit | streams, memory, error codes | https://docs.nvidia.com/cuda/ |
| cuBLAS (incl. cuBLASLt) | DeviceMatrix products, BLAS-1 | https://docs.nvidia.com/cuda/cublas/ |
| cuSOLVER | dense LLT / LU / QR / SVD / EVD | https://docs.nvidia.com/cuda/cusolver/ |
| cuSPARSE | SpMV / SpMM (CSC and BSR) | https://docs.nvidia.com/cuda/cusparse/ |
| cuFFT | gpu::FFT | https://docs.nvidia.com/cuda/cufft/ |
| NPP | device-side scalar and coefficient-wise arithmetic | https://docs.nvidia.com/cuda/npp/ |
| cuDSS | sparse direct solvers (separate install) | https://docs.nvidia.com/cuda/cudss/ |
Eigen is the linear algebra foundation for a large ecosystem of C++ projects in robotics (ROS, Drake, MoveIt, Pinocchio), computer vision (OpenCV, COLMAP, Open3D), scientific computing (Ceres, Stan), and beyond. Many of these projects run on GPU-equipped hardware but cannot use GPUs for Eigen operations without dropping down to raw CUDA library APIs.
GPU sparse solvers are a particularly acute gap. Sparse factorization is the bottleneck in SLAM, bundle adjustment, FEM, and nonlinear optimization – exactly the workloads where GPU acceleration matters most. Downstream projects like Ceres and COLMAP have open requests for GPU-accelerated sparse solvers, and third-party projects like cholespy exist specifically because Eigen lacks them. The contrib/Eigen/GPU module provides GPU sparse Cholesky, LDL^T, and LU factorization via cuDSS, alongside dense solvers (cuSOLVER), matrix products (cuBLAS), FFT (cuFFT), and sparse matrix-vector products (cuSPARSE).
Existing Eigen users should be able to move performance-critical dense or sparse linear algebra to the GPU with minimal code changes and without learning CUDA library APIs directly.
CPU and GPU coexist. There is no global compile-time switch that replaces CPU implementations (unlike EIGEN_USE_LAPACKE). Users choose GPU solvers explicitly – gpu::LLT<double> vs Eigen::LLT<MatrixXd>, gpu::SparseLLT<double> vs SimplicialLLT<SparseMatrix<double>> – and both coexist in the same binary. This also lets users keep the factored matrix on device across multiple solves, something impossible with compile-time replacement.
Familiar syntax. GPU operations use the same expression patterns as CPU Eigen. Here is a side-by-side comparison:
The GPU version reads like CPU Eigen with explicit upload/download for dense operations, and an almost identical API for sparse solvers. Expressions can copy-initialize a DeviceMatrix directly (as above), scalar factors accept plain literals (2 * d_A, d_A / 2, -d_A), and unsupported expressions are compile errors.
Standalone module. contrib/Eigen/GPU does not modify or depend on Eigen's Core expression template system (MatrixBase, CwiseBinaryOp, etc.). DeviceMatrix is not an Eigen expression type and does not inherit from MatrixBase. The expression layer is a thin compile-time dispatch where every supported expression maps to a single NVIDIA library call. There is no coefficient-level evaluation, lazy fusion, or packet operations.
Interoperability where useful. DeviceMatrix provides the same operator signatures as Matrix for common vector operations: +=, -=, *=, /=, dot(), squaredNorm(), norm(), stableNorm(), setZero(), noalias(), and copy construction and assignment (a device-to-device copy). This makes DeviceMatrix usable as a drop-in VectorType in Eigen algorithm templates that rely on these operations, and DeviceSparseView is a matrix-free matrix type for Eigen's iterative solvers: ConjugateGradient<gpu::DeviceSparseView<double>, Lower | Upper> runs Eigen's own algorithm, unmodified, on device vectors (see Eigen algorithm interop). Conjugate gradient is just the motivating example; we are open to expanding operator coverage as needed to support other high-level Eigen algorithms on the GPU.
Explicit over implicit. Host-device transfers, stream management, and library handle lifetimes are visible in the API. There are no hidden allocations or synchronizations except where documented (e.g., toHost() must synchronize to deliver data to the host).
gpu::DeviceMatrix<Scalar>A typed RAII wrapper for a dense column-major matrix in GPU device memory. This is the GPU counterpart of Eigen's MatrixX<Scalar>. A vector is simply a DeviceMatrix with one column. All public GPU classes live in namespace Eigen::gpu.
DeviceMatrix supports expression methods that mirror Eigen's API: adjoint(), transpose(), triangularView<UpLo>(), selfadjointView<UpLo>(), llt(), lu(). These return lightweight expression objects that are evaluated when assigned.
For BLAS Level-1 operations, DeviceMatrix also provides dot(), norm(), squaredNorm(), setZero(), noalias(), and arithmetic operators (+=, -=, *=) that dispatch to cuBLAS axpy, nrm2, dot, scal, and geam. These are the operations needed by iterative solvers.
gpu::DeviceScalar<Scalar>A device-resident scalar value. Reductions like dot(), norm(), and squaredNorm() return DeviceScalar instead of a host scalar, deferring the host synchronization until the value is actually needed:
Every reduction also has an overload that writes into an existing DeviceScalar on a given Context, d_x.dot(ctx, d_y, s), where s lives on ctx.stream(): it reuses the scalar's storage, so a loop, or a captured CUDA graph, repeats the reduction without allocating.
Division between DeviceScalar values (real types only) is performed on device via NPP, avoiding extra synchronizations. Small device allocations (including DeviceScalar) go through the stream-ordered allocator like every other block when the device has memory pools; on the cudaMalloc fallback path they are recycled through a thread-local DeviceBufferPool instead, to avoid cudaMalloc/cudaFree overhead in tight loops. Pool contract: a released block is recycled only after the device has retired every operation enqueued before the release on any blocking stream, so pooled buffers may move between the streams of one thread. The release is tracked by an event on the legacy default stream, the same ordering the stream-ordered allocator relies on. The pool is thread-local, so sharing a pooled buffer across threads needs external synchronization, and cudaStreamNonBlocking streams are outside the guarantee.
gpu::ContextEvery GPU operation needs a CUDA stream and library handles (cuBLAS eagerly, cuSOLVER / cuBLASLt / cuSPARSE lazily on first use). gpu::Context bundles these together. A single Context is not thread-safe – use one per thread (or external synchronization), since the underlying NVIDIA library handles are not thread-safe per handle.
For simple usage, you don't need to create one – a per-thread default context is created lazily on first use:
For concurrent multi-stream execution, create explicit contexts:
To integrate with existing CUDA code, borrow an existing stream:
To override the thread-local default (e.g., in CG where all ops share one context):
The module is header-only, but each feature pulls in the corresponding NVIDIA library at link time. cuSOLVER, cuBLASLt, and cuSPARSE are created lazily on first use, so a translation unit that only uses cuBLAS or cuFFT does not need to link the others:
| Feature | Link flags |
|---|---|
DeviceMatrix, GEMM, TRSM, SYMM, SYRK | -lcublas -lcublasLt |
| Dense solvers (LLT, LU, QR, SVD, EVD) | -lcusolver -lcublas |
FFT (gpu::FFT) | -lcufft -lcublas |
SpMV / SpMM (gpu::SparseContext) | -lcusparse -lcublas |
norm(), DeviceScalar arithmetic, /=, cwiseProduct | -lnpps -lnppc |
| Sparse direct solvers (cuDSS) | -lcudss -lcublas |
cuBLAS is required by DeviceMatrix itself (every Context creates a cuBLAS handle eagerly) and is also a runtime dependency of cuDSS, so it is the one constant. cuDSS additionally requires EIGEN_CUDSS to be defined before including contrib/Eigen/GPU.
Products dispatch to cuBLAS, GEMM through its cuBLASLt API; see Precision control for which compute type that selects.
Dot products, norms and vector arithmetic map to the corresponding cuBLAS Level-1 routines, except for device-side scalar arithmetic, which uses the signal-processing functions of NPP.
Backed by the dense part of cuSOLVER (cuSolverDN), whose documentation defines what each factorization returns and when it reports a numerical failure.
One-shot expression syntax – "one-shot" means factorization and solve run as a single fused call with no persistent factorization object; each evaluation re-factorizes:
Scratch for the one-shot form (factor copy, cuSOLVER workspace, info words) lives in the gpu::Context and grows monotonically — repeated one-shot solves perform no per-call allocations. In debug builds each call verifies the factorization status (one stream synchronization); release builds (EIGEN_NO_DEBUG/NDEBUG) skip the check and the sync, making the expression fully asynchronous — use the cached gpu::LLT / gpu::LU classes and info() when numerical failure must be detected.
cuSOLVER 11.4.1 and 11.4.2 (CUDA 11.8 and 12.0) have a defect here: cusolverDnXpotrf reports success on a matrix that is not positive definite and returns a factor full of NaN, so gpu::LLT::info() cannot detect that failure on those versions. cuSOLVER 11.4.4 (CUDA 12.1) and later report it.
Cached factorization – Factor once, solve many times:
The cached API keeps the factored matrix on device, avoiding redundant host-device transfers and re-factorizations. All five solvers accept compute(DeviceMatrix&&) to adopt the input and factor it in place with no copy (for QR/SVD with m < n the internal transpose still copies, and a view is always copied because its storage belongs to another object), and all five can bind to a gpu::Context to share its stream and handles. All solvers also accept host dense expressions directly as a convenience (e.g., gpu::LLT<double> llt(A) or qr.solve(B)), which handles upload/download internally. Host compute() finishes its upload before returning, while factorization remains asynchronous. The d_* accessors on gpu::SVD and gpu::SelfAdjointEigenSolver return non-owning DeviceMatrix views so downstream cuBLAS/cuSOLVER work can chain without round-tripping through host memory.
Requires cuDSS (separate install, CUDA 12.0+), which is distributed outside the CUDA Toolkit and versioned separately from it. Define EIGEN_CUDSS before including contrib/Eigen/GPU; see Linking for link flags.
gpu::SparseSolverConfig passes cuDSS tuning knobs through to the solver: fill-reducing reordering, matching, pivoting strategy / threshold / epsilon, iterative refinement, and the hybrid host/device memory and execute modes. Fields left at their defaults keep the cuDSS defaults, which favor speed over maximum robustness — for badly scaled or nearly singular systems, consider enabling matching and iterative refinement:
Each knob is consumed by the phase it affects (reordering/matching by analyzePattern(), pivoting by factorize(), refinement by solve()), so setConfig() must run before the first phase whose behavior it changes.
Every field is a pass-through, so which values are admissible for a given matrix type — and what each one does — is cuDSS's contract, not ours: cuDSS Data Types documents cudssConfigParam_t and the cudssReorderingAlg_t / cudssMatchingAlg_t / cudssPivotType_t values these enums mirror, and cuDSS Advanced Features describes the hybrid host/device memory and execute modes.
cuDSS < 0.8 names none of these algorithms. There the SparseReordering, SparseMatching and SparsePivoting enumerators other than Default are not declared, so selecting one is a compile error rather than a request the linked cuDSS cannot honor. The remaining fields — thresholds, refinement, the hybrid modes — still exist, and setConfig() refuses any non-default value of them: it asserts, and info() reports InvalidInput until the config is reset to default, so the request cannot be silently downgraded to the cuDSS defaults. EIGEN_HAS_CUDSS_SOLVER_CONFIG is 1 or 0 accordingly, for callers that need to branch at compile time.
Plans and data layouts are cuFFT's. The scaling convention is not: cuFFT leaves its transforms unnormalized, and gpu::FFT applies the 1/n on the inverse so that inv(fwd(x)) == x, matching contrib/Eigen/FFT.
Uses the cuSPARSE generic API (cusparseSpMV / cusparseSpMM), which fixes the supported index and value type combinations.
Host-input calls re-upload the sparse values and index arrays on every call (host pointer identity cannot detect a pattern rewritten in place or assigned into the same allocations, so the structure is never assumed unchanged). The cuSPARSE descriptors and workspace-size queries are cached across calls with matching shapes; deviceView() is the upload-once path. A DeviceSparseView carries a generation counter — using a view after any later upload through its context asserts instead of silently multiplying by the wrong matrix.
Every SparseContext entry point above also accepts a BlockSparseMatrix. Square blocks of size at least 2 upload in cuSPARSE's BSR (block sparse row) format on cuSPARSE 12.6.3 (CUDA 13.0 Update 1) or newer, where the generic SpMV and SpMM run on BSR descriptors; EIGEN_HAS_CUSPARSE_BSR is 1 there and 0 on older toolkits. Any other block shape, and every BlockSparseMatrix on an older toolkit, takes the CSC path instead: the matrix is expanded with toSparse() on the host, once per host-input call or once per deviceView(), so source compatibility does not depend on the toolkit.
cuSPARSE multiplies a BSR descriptor only as op == NoTrans with row-major blocks (CUSPARSE_ORDER_ROW), so the op is applied on the host rather than passed to the library. A RowMajor matrix is BSR of itself and uploads without a host copy for NoTrans; a ColMajor one — column-major blocks in block-column order — is BSR of its transpose and uploads without a copy for Trans (and ConjTrans on real scalars). Every other op / storage-order combination transposes (and conjugates) the matrix on the host first, once per host-input call or once per deviceView(). spmv_device_exec() and spmm_device_exec() consequently accept only NoTrans against a BSR upload (debug builds assert); a BlockSparseMatrix on the CSC path passes the op to cuSPARSE like a SparseMatrix.
The CSC fallback exists because cuSPARSE rejects rectangular blocks at descriptor creation and runs no BSR SpMV on 1 x 1 blocks; internal::use_cusparse_bsr<BlockRows, BlockCols> is the exact selector. The int index type and the int limits on dimensions and nonzeros are as for SparseMatrix.
Eigen's ConjugateGradient runs on the GPU types. DeviceSparseView is a matrix-free matrix type (it inherits EigenBase and carries SparseMatrix traits, so IterativeSolverBase holds it by pointer), the vectors are DeviceMatrix, and the algorithm's VectorType (Dest::PlainObject) is DeviceMatrix itself. solve() and solveWithGuess() return Eigen expressions and need Eigen operands; solveWithGuessInPlace(b, x) is the entry point for device vectors, with x holding the initial guess. The identity preconditioner works as is; DiagonalPreconditioner and the other Eigen preconditioners evaluate host expressions and do not. compute() stores a pointer to the view, so the view and the SparseContext behind it must outlive the solver. The class form is real-Scalar only: DeviceScalar arithmetic covers real types, and numext::real() of a complex DeviceScalar has no host conversion to RealScalar.
The algorithm reads three values on the host per iteration – alpha, the stableNorm() convergence check and absNew – so this form synchronizes three times per iteration. The division absNew / p.dot(tmp) resolves to the device-side operator/(Scalar, DeviceScalar): absNew is uploaded into a fresh DeviceScalar, divided through NPP and read back, two small allocations and a kernel launch per iteration on top of the sync. The hand-written loop below is the same algorithm with one host sync per iteration (the convergence check); all scalar intermediates (alpha, beta, absNew) stay on device as DeviceScalar values. Its convergence test squares the residual norm (squaredNorm(), a cuBLAS dot), which overflows for ||r|| > sqrt(max) and underflows to 0 for ||r|| < sqrt(min), so it lacks the extreme-scale robustness the template gets from stableNorm() and its residual scaling:
GEMM dispatch routes through cublasLtMatmul. The compute type is selected per scalar via the cuda_compute_type trait in CuBlasSupport.h, gated by two compile-time macros:
| Macro | Effect |
|---|---|
| (default) | CUBLAS_COMPUTE_32F / CUBLAS_COMPUTE_64F. cublasLt heuristics may pick tensor-core algorithms; on sm_80+ doubles can land on Ozaki-emulated tensor cores. |
EIGEN_CUDA_TF32 | CUBLAS_COMPUTE_32F_FAST_TF32 for float and complex<float> (~2x faster, 10-bit mantissa). No effect on double / complex<double>. |
EIGEN_NO_CUDA_TENSOR_OPS | Pedantic compute types (CUBLAS_COMPUTE_*_PEDANTIC) for every scalar — disables tensor-core algorithms. Use for bit-exact reproducibility. Takes precedence over EIGEN_CUDA_TF32. |
These are independent of cuBLAS's runtime cublasSetMathMode() / CUBLAS_TF32_OVERRIDE controls; the cublasLt path keys off the compile-time compute type instead. The cublasGemmEx fallback (used when cublasLt's heuristic returns no candidate) honors EIGEN_NO_CUDA_TENSOR_OPS via its algorithm hint (CUBLAS_GEMM_DEFAULT vs CUBLAS_GEMM_DEFAULT_TENSOR_OP).
Operations are asynchronous by default. The compute-solve chain runs without host synchronization until you need a result on the host:
Mandatory sync points:
fromHost() – Synchronizes to complete the upload before returningtoHost() / HostTransfer::get() – Must deliver data to hostinfo() – Must read the factorization statusDeviceScalar implicit conversion – Downloads scalar from deviceDebug-only sync points (compiled out under EIGEN_NO_DEBUG/NDEBUG): every solver solve() and accessor verifies info() == Success via eigen_assert, which forces one stream synchronization the first time after each compute()/factorize(). Release builds perform no such check — call info() explicitly where failure detection matters. gpu::SVD's device solve additionally downloads the singular values once per (truncation, lambda) setting to build its cached inverse diagonal.
Cross-stream safety is automatic. DeviceMatrix tracks write completion via CUDA events. When a matrix written on stream A is read on stream B, the module automatically inserts cudaStreamWaitEvent. Same-stream operations skip the wait (CUDA guarantees in-order execution within a stream).
Device memory allocation is stream-ordered. All module allocations go through cudaMallocAsync / cudaFreeAsync on devices that support memory pools (detected at runtime; cudaMalloc/cudaFree fallback otherwise, or force the fallback with EIGEN_GPU_NO_STREAM_ORDERED_ALLOC — required when borrowing cudaStreamNonBlocking streams, which do not synchronize with the legacy stream the allocator uses for ordering). Consequences:
DeviceMatrix temporaries no longer performs a device-wide synchronization; freed blocks recycle through the driver pool (the pool's release threshold is raised so steady-state loops reallocate at user-space speed).DeviceMatrix) with work still in flight is safe and async: the stream-ordered free waits for previously enqueued work without stalling the host.DeviceMatrix::resize() is capacity-aware: shrinking or same-size reshapes reuse the existing allocation (contents are still discarded).Every CUDA runtime and library call the module makes is checked, in release builds as in debug builds. A failed call prints file:line: call: error to stderr and stops the program:
std::abort() when EIGEN_NO_DEBUG (or NDEBUG) is defined;eigen_assert otherwise.error is the status name where the library provides one (cudaErrorInvalidValue, CUBLAS_STATUS_INVALID_VALUE, ...) and <library> status <code> for cuFFT, cuDSS, NPP and cuBLAS before 11.6.1. There is no mode that ignores a failure: the failed call has not done its work, and a sticky error (say, an illegal address in a kernel) makes every later call on the device fail as well.
To handle failures yourself, define EIGEN_GPU_CHECK_FAILED(error, expression, file, line) before including the module. error, expression and file are C strings; line is an int. For example, to turn failures into exceptions:
Destructors release their resources without the checks, so a throwing handler never runs inside one. A throw also restores the library-handle state an operation changes temporarily (the cuBLAS pointer mode), so the context stays usable; only the interrupted operation's output is unspecified. A handler that returns lets execution continue past the failed call, which is only useful in tests.
Numerical failures (a matrix that is not positive definite, a singular factorization) are not call failures: they are reported by info() as described above.
float, double, std::complex<float>, std::complex<double> (unless noted otherwise).
| DeviceMatrix expression | Library call | Parameters |
|---|---|---|
C = A * B | cublasLtMatmul (with cublasGemmEx fallback) | transA=N, transB=N, alpha=1, beta=0 |
C = A.adjoint() * B | cublasLtMatmul | transA=C, transB=N |
C = A.transpose() * B | cublasLtMatmul | transA=T, transB=N |
C = A * B.adjoint() | cublasLtMatmul | transA=N, transB=C |
C = A * B.transpose() | cublasLtMatmul | transA=N, transB=T |
C = alpha * A * B | cublasLtMatmul | alpha from LHS |
C = A * (alpha * B) | cublasLtMatmul | alpha from RHS |
C += A * B | cublasLtMatmul | alpha=1, beta=1 |
C -= A * B | cublasLtMatmul | alpha=-1, beta=1 |
X = A.llt().solve(B) | cusolverDnXpotrf + Xpotrs | uplo, n, nrhs |
X = A.llt<Upper>().solve(B) | same | uplo=Upper |
X = A.lu().solve(B) | cusolverDnXgetrf + Xgetrs | n, nrhs |
X = A.triangularView<L>().solve(B) | cublasXtrsm | side=L, uplo, diag=NonUnit |
C = A.selfadjointView<L>() * B | cublasXsymm / cublasXhemm | side=L, uplo |
C.selfadjointView<L>().rankUpdate(A) | cublasXsyrk / cublasXherk | uplo, trans=N |
C = A + B | cublasXgeam | alpha=1, beta=1 |
C = A + alpha * B | cublasXgeam | alpha=1, beta from scaled |
C = A - B | cublasXgeam | alpha=1, beta=-1 |
C = A - alpha * B | cublasXgeam | alpha=1, beta=-scaled |
C = alpha * A + beta * B | cublasXgeam | both sides scaled |
C = alpha * A, C = -A, C = A / alpha | cublasXgeam | beta=0, aliasing-safe |
x += alpha * y | cublasXaxpy | alpha (host scalar) |
x += dAlpha * y | cublasXaxpy | alpha (DeviceScalar, device pointer mode) |
x -= alpha * y | cublasXaxpy | alpha negated |
x *= alpha | cublasXscal | alpha (host or DeviceScalar) |
x.dot(y) | cublasXdot / cublasXdotc | returns DeviceScalar |
x.norm() | cublasXdot(x, x), then nppsSqrt | as squaredNorm(), then its square root on device |
x.stableNorm() | cublasXnrm2 | returns DeviceScalar<RealScalar> |
x.squaredNorm() | cublasXdot(x, x) | real dot over the 2n real and imaginary parts for complex x; returns DeviceScalar<RealScalar> |
d_y = view * d_x | cusparseSpMV | device-resident SpMV |
d_Y = view * d_X | cusparseSpMM | device-resident SpMM (RHS with >1 column) |
same, view of a BlockSparseMatrix | cusparseSpMV / cusparseSpMM on a BSR descriptor | opA=N, row-major blocks; op(A) formed on the host |
DeviceMatrix<Scalar>Typed RAII wrapper for a dense column-major matrix in GPU device memory. Always dense (leading dimension = rows). A vector is a DeviceMatrix with one column.
DeviceScalar<Scalar>Device-resident scalar. Returned by dot(), norm(), and squaredNorm(). Implicit conversion to Scalar triggers cudaStreamSynchronize + download.
gpu::ContextUnified GPU execution context owning a CUDA stream and library handles. Not thread-safe – use one Context per thread, or external synchronization across threads.
Non-copyable, non-movable (owns library handles). Translation units that never call cusolverHandle() do not pull cuSOLVER symbols at link time – see Linking.
gpu::LLT<Scalar, UpLo> – Dense Cholesky (cuSOLVER)Caches the Cholesky factor on device for repeated solves.
gpu::LU<Scalar> – Dense LU (cuSOLVER)Same pattern as gpu::LLT. Adds a gpu::GpuOp parameter on solve().
gpu::GpuOp: NoTrans, Trans, ConjTrans.
gpu::QR<Scalar> – Dense QR (cuSOLVER)QR factorization via cusolverDnXgeqrf. Solve uses ORMQR (apply Q^H) + TRSM (back-substitute on R) – Q is never formed explicitly.
gpu::SVD<Scalar> – Dense SVD (cuSOLVER)SVD via cusolverDnXgesvd. Supports ComputeThinU | ComputeThinV, ComputeFullU | ComputeFullV, or 0 (values only). Wide matrices (m < n) handled by internal transpose.
Note: singularValues(), matrixU(), matrixV(), and matrixVT() download to host on each call. The d_* accessors return non-owning DeviceMatrix views into the solver's internal buffers; the gpu::SVD object must outlive any view derived from it. For wide matrices (m < n) the U/V^T views are owning (one cublasXgeam adjoint pass).
gpu::SelfAdjointEigenSolver<Scalar> – Eigendecomposition (cuSOLVER)Symmetric/Hermitian eigenvalue decomposition via cusolverDnXsyevd. options: ComputeEigenvectors (the default) or EigenvaluesOnly.
Note: eigenvalues() and eigenvectors() download to host on each call. The d_* accessors return non-owning DeviceMatrix views into the solver's internal buffers; the gpu::SelfAdjointEigenSolver object must outlive any view derived from it.
HostTransfer<Scalar>Future for async device-to-host transfer. Returned by DeviceMatrix::toHostAsync().
gpu::SparseLLT<Scalar, UpLo> – Sparse Cholesky (cuDSS)Requires cuDSS (CUDA 12.0+, #define EIGEN_CUDSS). Three-phase workflow with symbolic reuse. Accepts SparseMatrix<Scalar, ColMajor, int> (CSC). Matrix dimensions and nonzero count must fit in int (cuDSS limitation; debug builds assert).
All three cuDSS solvers also accept a gpu::Context& as first constructor argument to borrow its stream (gpu::SparseLLT<double> llt(ctx) or llt(ctx, A)).
gpu::SparseLDLT<Scalar, UpLo> – Sparse LDL^T (cuDSS)Symmetric indefinite. Same API as gpu::SparseLLT.
gpu::SparseLU<Scalar> – Sparse LU (cuDSS)General non-symmetric. Same API as gpu::SparseLLT (without UpLo).
gpu::FFT<Scalar> – FFT (cuFFT)Plans cached by (size, type) in a bounded LRU and reused; the least-recently-used plan is destroyed via cufftDestroy on overflow. Cache capacity is set at construction (default kDefaultCufftPlanCacheCapacity = 16). Inverse transforms scaled so inv(fwd(x)) == x. Supported scalars: float, double. Stream and cuBLAS handle borrowed from a gpu::Context (default: Context::threadLocal()), so by default the FFT shares a stream with other GPU operations on the same thread.
All FFT methods accept host data and return host data. Upload/download is handled internally. The C2C and R2C overloads of fwd() are distinguished by the input scalar type (complex vs real).
gpu::SparseContext<Scalar> – SpMV/SpMM (cuSPARSE)Accepts SparseMatrix<Scalar, ColMajor> and BlockSparseMatrix<Scalar, Options, BlockRows, BlockCols, int> (the BlockSpMat<Options, BlockRows, BlockCols> alias; see Block sparse matrices for which block shapes upload as BSR and which op / storage-order combinations do so without a host copy). Host-input methods accept host data and return host data; device-input methods (deviceView(), multiply(A, d_x, d_y)) operate on DeviceMatrix. Matrix dimensions and nonzero count must fit in int (cuSPARSE limitation; debug builds assert).
DeviceSparseView<Scalar> – Device-resident sparse matrixReturned by gpu::SparseContext::deviceView(). Holds a sparse matrix on device for repeated SpMV without re-uploading.
Unlike Eigen's Matrix, where omitting .noalias() triggers a copy to a temporary, DeviceMatrix dispatches directly to NVIDIA library calls which have no built-in aliasing protection. All operations are implicitly noalias. The caller must ensure operands don't alias the destination for GEMM, TRSM, SYMM/HEMM, and SYRK/HERK. Debug builds assert on these violations before dispatching to cuBLAS. geam expressions (d_C = d_A + alpha * d_B) are safe with aliasing. The .noalias() method exists as a no-op for Eigen template compatibility.
compute(MatrixXd), solve(MatrixXd)) and device-input (compute(DeviceMatrix), solve(DeviceMatrix)) overloads, plus host- and device-side accessors (matrixU() vs d_matrixU()). This eases migration from CPU Eigen but may invite accidental host ↔ device round-trips when users mix the two without realising the cost. Revisit once the module is in users' hands; if the convenience overloads cause more confusion than they save, narrow toward a single explicit fromHost / toHost boundary.gpu::SparseSolverConfig exposes matching, pivoting, refinement, and the hybrid modes, but a default-constructed solver still runs with cuDSS's performance-tuned defaults (matching off). Consider flipping the shipped defaults toward robustness now that users can override them.gpu::SparseLDLT treats complex inputs as Hermitian (matching Eigen::SimplicialLDLT). cuDSS also supports CUDSS_MTYPE_SYMMETRIC for complex matrices (A = A^T, no conjugation); exposing this would need a separate solver mode.| File | Depends on | Contents |
|---|---|---|
GpuSupport.h | <cuda_runtime.h> | EIGEN_GPU_CHECK_FAILED, runtime error macro, DeviceBuffer, DeviceBufferPool, cuda_data_type<> |
DeviceMatrix.h | GpuSupport.h | gpu::DeviceMatrix<>, gpu::HostTransfer<> |
DeviceExpr.h | DeviceMatrix.h | GEMM, geam, and device-scalar expression wrappers |
DeviceBlasExpr.h | DeviceMatrix.h | TRSM, SYMM, SYRK expression wrappers |
DeviceSolverExpr.h | DeviceMatrix.h | Solver expression wrappers (LLT, LU) |
DeviceScalar.h | GpuSupport.h, DeviceScalarOps.h | gpu::DeviceScalar<> (device-resident scalar) |
DeviceScalarOps.h | <npps_*.h> | Scalar div/neg/sqrt/cwiseProduct via NPP, NPP error macro |
DeviceDispatch.h | all above | All dispatch functions, BLAS-1 out-of-line defs, gpu::Assignment |
GpuContext.h | CuBlasSupport.h, CuSolverSupport.h, CuSparseSupport.h | gpu::Context |
CuBlasSupport.h | GpuSupport.h, <cublas_v2.h>, <cublasLt.h> | cuBLAS error macro, type-specific wrappers |
CuSolverSupport.h | GpuSupport.h, <cusolverDn.h> | cuSOLVER params, fill-mode mapping |
GpuSolverContext.h | CuSolverSupport.h, CuBlasSupport.h | Shared solver context (stream, handles, scratch) |
GpuLLT.h | GpuSolverContext.h | gpu::LLT<> – Cached dense Cholesky factorization |
GpuLU.h | GpuSolverContext.h | gpu::LU<> – Cached dense LU factorization |
GpuQR.h | GpuSolverContext.h | gpu::QR<> – Dense QR decomposition |
GpuSVD.h | GpuSolverContext.h | gpu::SVD<> – Dense SVD decomposition |
GpuEigenSolver.h | GpuSolverContext.h | gpu::SelfAdjointEigenSolver<> |
CuFftSupport.h | GpuSupport.h, <cufft.h> | cuFFT error macro, type-dispatch wrappers |
GpuFFT.h | CuFftSupport.h, CuBlasSupport.h, GpuContext.h | gpu::FFT<> – 1D/2D FFT with plan caching |
CuSparseSupport.h | GpuSupport.h, <cusparse.h> | cuSPARSE error macro, EIGEN_HAS_CUSPARSE_BSR |
GpuSparseContext.h | CuSparseSupport.h | gpu::SparseContext<>, gpu::DeviceSparseView<>, BSR binding of BlockSparseMatrix |
CuDssSupport.h | GpuSupport.h, <cudss.h> | cuDSS error macro, type traits (optional) |
GpuSparseSolverBase.h | CuDssSupport.h | CRTP base for sparse solvers (optional) |
GpuSparseLLT.h | GpuSparseSolverBase.h | gpu::SparseLLT<> – Sparse Cholesky via cuDSS (optional) |
GpuSparseLDLT.h | GpuSparseSolverBase.h | gpu::SparseLDLT<> – Sparse LDL^T via cuDSS (optional) |
GpuSparseLU.h | GpuSparseSolverBase.h | gpu::SparseLU<> – Sparse LU via cuDSS (optional) |
DeviceBatchMatrix). A strided batch of N identical-size matrices dispatching to cuBLAS/cuSOLVER batched APIs (cublasDgemmBatched, cusolverDnXpotrfBatched, etc.). This enables robotics and model-predictive control workloads where many small independent systems are solved in parallel.TensorContractionGpu.h / TensorReductionGpu.h) with cuTENSOR dispatch, following the same library-dispatch pattern used by contrib/Eigen/GPU.cudaMallocManaged or cudaHostAllocMapped to eliminate fromHost() / toHost() copies on integrated GPUs (Jetson) where CPU and GPU share DRAM.DeviceMatrix dispatch and device-side Eigen expression templates (Core + Tensor) running inside CUDA kernels. Raw-pointer + Map / TensorMap as the zero-copy interop surface.cudaMemPool_t per stream (cudaDeviceSetMempool / cudaMallocFromPoolAsync) could further reduce cross-stream allocator contention for workloads that fan out many concurrent solves. 
1.13.2