|
Uni20
Asynchronous tensor and runtime library.
|
Status: implemented BLAS and initial LAPACK checkpoints plus forward plan.
The direct descriptors below describe the native provider ABI. Optional MPFR/MPC kernels and binary80 kernels use the same BlasBackend/LapackBackend families, but pack logical values into owning provider objects instead of using native pointers directly. Their current coverage includes GEMM, square solves, reusable LU, QR/LQ, self-adjoint and general eigensystems, and SVD; the integration roadmap tracks further routine and workspace support.
This layer is below tensor front-end dispatch and above the existing raw BLAS/LAPACK provider facades. Direct mdspan GEMM/GEMV, accessor-respecting CPU fallbacks, explicit Tensor copy/materialization, and strict LAPACK projected eigensystem/Schur plus dense self-adjoint adapters now use the operation-tag dispatcher.
The implemented mdspan-facing wrappers accept resolved views, build provider-compatible descriptors, and call the existing Uni20 BLAS/LAPACK backend wrappers. Future work extends this adapter layer operation by operation; it does not replace the backend facades.
Related notes:
Keep the existing wrapper layers and add one new layer above them:
The existing backend wrappers remain provider and ABI facades:
src/uni20/backend/blas: BLAS discovery, vendor selection, LP64/ILP64 configuration, raw BLAS declarations, and vendor extension targets.src/uni20/backend/lapack: checked and unchecked LAPACK overloads over the configured provider.The new mdspan layer owns view interpretation:
Do not move shape/materialization policy into src/uni20/backend. That layer should stay close to provider ABI calls. The mdspan layer may call uni20::blas::* and uni20::lapack::*, but it should not redeclare Fortran symbols or include provider-specific headers directly.
The direct BLAS descriptor implementation lives under src/uni20/linalg/blas/:
blas_matrix.hppBlasReadableMatrix and BlasWritableMatrix operands, backend-independent transform algebra, and BLAS transpose lowering.blas_vector.hppBlasReadableVector and BlasWritableVector operands.mdspan_access.hppmdspan_matrix.hpptry_mdspan_matrix_stage(...), role-specific readable/writable lowering helpers, and strict column-major-compatible helpers for LAPACK update operands.mdspan_vector.hpptry_mdspan_vector_stage(...), BLAS increment lowering, and readable or writable rank-one descriptors, including strict contiguous writable LAPACK vectors.gemm.hppgemv.hppblas.hppThe first operation-tag dispatch slice adds:
operation_tags.hppdispatch.hppKernelTypeAcceptance, backend_list, selector normalization, type acceptance detection, and the runtime backend walk.ops/gemm.hppgemm(...).ops/contract.hppops/gemv.hppgemv(...) with rank-1 output/input and a rank-2 matrix.backends/blas/gemm.hppBlasBackend and try_kernel(BlasBackend, gemm_op, ...).contract_op grouping, host acquisition, and GEMM lowering.backends/blas/gemv.hppBlasBackend and try_kernel(BlasBackend, gemv_op, ...).backends/cpu/gemm.hppCpuReferenceBackend tensor-view adapter for reference GEMM.backends/cpu/gemv.hppCpuReferenceBackend tensor-view adapter for reference GEMV.cpu/gemm.hpp and cpu/gemv.hppLAPACK operation adapters live under src/uni20/linalg/backends/lapack/; currently this includes tridiagonal and nonsymmetric eigensystems, Schur and Hessenberg Schur decomposition, Schur reordering, and real symmetric/complex Hermitian syev/heev. Shared provider-descriptor construction remains under src/uni20/linalg/blas/. Keep that division so descriptor helpers remain testable without running a LAPACK operation.
Keep two descriptor levels separate:
Do not merge these roles. The first one is an ABI object; the second one is a planning object.
BlasWritableMatrix is the small object that can be passed to wrappers around uni20::blas::* and uni20::lapack::* for an output/update matrix. Its rows and cols are the matrix dimensions before applying any BLAS transpose character. They are not necessarily the logical mdspan dimensions.
BlasReadableMatrix has the same fields plus a MatrixTransform. This is the provider-ready input operand for BLAS calls.
Writable BLAS outputs do not carry a transform. If an output mdspan is row-major-like, the operation must be rewritten so the provider writes the untransformed BlasWritableMatrix, or the wrapper must use an explicit output temporary with copy-back policy.
Do not use inheritance or a common contained base for these two objects. The duplicated fields are intentional: the objects are small, and implicit conversion from a readable object to a common base would lose the transform.
BlasReadableVector and BlasWritableVector similarly duplicate data, size, and increment. The readable vector has no transform field because portable GEMV cannot conjugate x; MdspanVectorStage::needs_conjugation therefore causes a direct GEMV decline instead of being discarded.
The staging descriptor keeps the logical extents in mdspan axis order and records which mdspan axis has stride 1. The non-unit stride is canonicalized for extent-1 axes so it also satisfies the BLAS leading-dimension rule. It does not store BLAS transposition state.
This chooses the mdspan-storage-first policy rather than a BLAS-first policy. The alternative would store BLAS rows, cols, leading_dimension, and a transpose flag. That is too easy to make ambiguous: once a view is marked as transposed, it is unclear at the staging boundary whether rows and cols describe the logical matrix or the matrix seen by the BLAS ABI. In this plan, the staging descriptor always describes logical mdspan axes; helper functions derive the BLAS ABI reference.
There is only one staging descriptor. needs_conjugation defaults to false and records that the logical matrix is conjugated relative to the raw storage interpretation. For readable lowering, it composes into BlasReadableMatrix::transform. For writable lowering, it is not represented in BlasWritableMatrix itself. The current direct GEMM wrapper rejects such a stage as a user-provided output view. This is separate from a prepared operation internally conjugating the output storage before and after a provider call as an implementation trick.
The staging helper derives needs_conjugation from the mdspan accessor policy. It is not an independent option supplied by the BLAS adapter. The generic mdspan conj(...) helper in src/uni20/mdspan/conjugate_accessor.hpp returns a read-only conjugating accessor view for complex mdspans, cancels to the const original view when applied twice, and returns a const identity view for non-complex mdspans. The same header defines accessor_applies_conjugation_v<Accessor> and mdspan_needs_conjugation_v<View>. The implemented conj(Tensor) operation returns a tensor view whose mdspan uses that conjugation accessor rather than eagerly materializing conjugated storage. When that view reaches the BLAS staging helper, try_mdspan_matrix_stage(...) inspects the accessor and sets needs_conjugation = true.
The conjugating accessor follows the C++26 std::linalg::conjugated_accessor direction in WG21 P3050R3. Uni20 deliberately keeps uni20::conj(span) as the user-facing operation and does not adopt conj-if-needed terminology: for non-complex values, uni20::conj is the identity on values and the mdspan helper returns a read-only identity view.
The direct BLAS path should accept only accessors that can still expose the raw storage handle plus declarative metadata such as needs_conjugation. An accessor that computes arbitrary transformed values by value is a generic transform view, not directly BLAS-addressable; it needs materialization or a non-BLAS fallback. This follows the mdspan accessor model: mapping and accessor describe how logical elements are read from the same handle, and the BLAS adapter may lower only the subset it understands. The broader view policy lives in docs/tensor/dispatch_and_view_semantics_draft.md: structural views such as real(x) and imag(x) are slice-like address transformations, while semantic transform views such as conj(x) are read-only by default.
A pointer data_handle_type is not enough to prove direct BLAS readability or writeability. Direct BLAS/LAPACK wrappers may bypass access(...); that is only correct for stdex::default_accessor<T> or for accessors whose semantics the wrapper explicitly lowers. For the current BLAS adapter, readable inputs may use stdex::default_accessor<T> or Uni20's conjugated_accessor over a default accessor. Writable outputs and LAPACK update operands require default-accessor views.
Fixed-output GEMM requires the output element addresses to be disjoint from both readable inputs. Fixed-output GEMV similarly requires the output vector to be disjoint from the matrix and input vector. Readable operands may overlap each other. These are leaf-kernel preconditions: supporting an aliased output requires an explicit temporary and copy-back policy above the direct wrapper.
Do not reject operands merely because their enclosing allocation or bounding address ranges overlap. Two slices of one matrix may still address disjoint elements. A directly BLAS-compatible matrix is an ordered union of contiguous unit-stride intervals, one per provider column, so a diagnostic or binding layer can test exact overlap in linear time in the number of provider columns. Generic accessor-respecting CPU views do not necessarily admit an equally cheap exact test.
The C++ leaf API may rely on this documented precondition. Python and other runtime-erased front ends should validate overlap whenever it is exactly representable and should reject or materialize uncertain cases according to their user-facing policy.
unit_stride_axis is deliberately an integer mdspan axis, not a row/column enum. Names such as row and column are easy to misread because stride(0) == 1 is the column-major-compatible case, while stride(1) == 1 is the row-major-like case. The descriptor stores the facts about the view and leaves interpretation to helper functions.
For a logical matrix A(i, j):
extent0 is the logical row extent.extent1 is the logical column extent.unit_stride_axis == 0 means stride(0) == 1 and nonunit_stride == stride(1). The BLAS ABI sees the logical matrix directly.unit_stride_axis == 1 means stride(1) == 1 and nonunit_stride == stride(0). The BLAS ABI sees the logical transpose.When the corresponding extent is 1, the non-unit stride is not observable in the logical mdspan mapping. The staging helper may therefore choose a canonical representative that satisfies BLAS, rather than preserving the arbitrary value reported by the mapping.
Runtime descriptor construction should reject views where neither matrix stride is 1. The first pass should also reject negative strides for direct BLAS/LAPACK calls. Any extent or stride that cannot be represented by the configured signed blas_int ABI type must be rejected by the direct wrapper; larger logical problems require an ILP64 provider or a higher-level blocked operation.
For direct BLAS lowering, MdspanMatrixStage::nonunit_stride must satisfy the provider ABI lower bound:
unit_stride_axis == 0: nonunit_stride >= max(1, extent0);unit_stride_axis == 1: nonunit_stride >= max(1, extent1).Degenerate dimensions need explicit tests because a non-unit stride may be arbitrary when it corresponds to an extent-1 provider column axis. If the provider column count is at most one, normalize nonunit_stride upward to max(1, provider_rows) instead of declining the view. If the provider would actually step through the non-unit stride, do not repair it by changing the address arithmetic; reject the view unless the stride already satisfies the provider contract.
Provider-ready BLAS operands are still derived from the mdspan staging descriptor:
These helpers return provider-ready objects with rows = extent0, cols = extent1, and leading_dimension = nonunit_stride when unit_stride_axis == 0; they swap the extents when unit_stride_axis == 1. For readable operands, the returned transform is the effective transform after composing the storage transform and stage.needs_conjugation. Callers should express logical conjugation by passing uni20::conj(span) and should express logical transpose by passing a transposed mdspan view. For writable operands, blas_writable_matrix(...) returns the provider-writeable storage shape only. A stage with needs_conjugation == true is not accepted by the current direct GEMM wrappers: try_gemm(...) declines it and the checked gemm(...) wrapper aborts. In normal use this is unreachable because uni20::conj(mdspan) returns a const readable view, but the runtime guard keeps custom accessors honest. This keeps the staging descriptor unambiguously mdspan-like while making the BLAS ABI conversion local and testable.
Represent transforms independently of the backend-specific character, using the MatrixTransform enum above.
The transform helper layer should provide:
storage_transform(MdspanMatrixStage const&);transpose_result_transform(MatrixTransform);compose(MatrixTransform, MatrixTransform);logical_rows(...), logical_cols(...), abi_rows(...), abi_cols(...);blas_trans_char_is_supported(MatrixTransform) -> bool;blas_trans_char(MatrixTransform) -> char.For readable operands, compose the storage transform and stage.needs_conjugation. Operation-specific rewrites, such as row-major GEMM output normalization, may compose additional internal transforms before the provider call. User-facing mdspan wrappers should not take MatrixTransform parameters; transform intent comes from the mdspan view/accessor itself. The transform spelling helper maps the internal transform alphabet to ‘'N’,'T', 'C', and'R', where'R'means conjugate without transpose. The current generic direct GEMM wrapper still gates provider calls to the portable'N', 'T', and'C'subset, so complex conjugate-only states are declined before side effects. OpenBLAS documents the extension in its CBLAS enum as CblasConjNoTrans, and itsinterface/gemm.cdispatcher also recognizes the Fortran-style'R'spelling on the develop branch: <https://github.com/OpenMathLib/OpenBLAS/blob/develop/interface/gemm.c>. Intel MKL's Fortran GEMM rejected'R'in local testing, so any direct'R'fast path should be explicitly provider-gated rather than assumed as baseline BLAS. Otherwise a prepared wrapper must materialize the conjugated readable operand and dispatch through the ordinaryN/T/C` path.
storage_transform(view) returns normal when unit_stride_axis == 0 and transpose when unit_stride_axis == 1. The helper name should describe the derived BLAS interpretation; the stored descriptor remains mdspan-axis based. There should be no stored transposed field.
If L is the logical mdspan matrix and S is the untransformed BlasWritableMatrix shape derived from the same stage, then L = storage_transform(stage)(S). The readable provider operand is lowered as:
where o means function composition. For example, a row-major-like input stage has storage_transform == transpose, so it lowers to BLAS ‘'T’. If the same view is wrapped byuni20::conj(...)`, the accessor contributes conjugation and the effective transform becomes conjugate-transpose. A caller that wants a logical transpose should pass a structural transposed mdspan view rather than a separate transform flag.
Writable outputs cannot carry a BLAS transpose flag. A row-major-like writable matrix therefore needs an operation rewrite, not a flag on the output.
For GEMM:
if C has unit_stride_axis == 1, rewrite to:
The helper should swap readable operands and apply an internal transpose transform to both staged readable operands. This keeps the writable output untransformed while still describing the logical result.
This rewrite is appropriate for BLAS update-style operations. It is not a general rule for LAPACK drivers that overwrite inputs with structured outputs.
BlasWritableMatrix intentionally has no transform field. It is only the provider-writeable storage target. Operation-specific wrappers normalize any output-side transform before dispatch.
There are two different notions that should stay separate:
uni20::conj(mdspan), and that view is read-only.N/T/C provider character set. It is not user-visible and does not require a conjugate_output flag on the direct wrapper.For example, after composing row-major storage, input view conjugation, and operation-specific normalization, the best lowering may require a conjugate-only readable transform that the baseline provider cannot encode. A future prepared GEMM fallback may:
beta != 0, so the accumulated old-output term represents conj(C_old);alpha, beta, and all readable operands logically, which flips each readable transform's conjugation bit and may turn the unsupported conjugate-only case into an ordinary N/T/C provider call;When beta == 0, the pre-call output conjugation can be skipped because the old output is not read. Uni20 provides the accessor-respecting conjugate_inplace_op CPU and CUDA reference kernels for this eager storage rewrite. The CUDA path supports positive-strided cfloat and cdouble storage; real and integer conjugation is a no-op. The fallback must still be used only when the transformed readable operands are all representable by the selected provider. If conjugating everything merely moves an unsupported conjugate-only transform from one input to another, the prepared wrapper still needs input materialization or must decline.
Direct no-copy wrappers should not silently materialize output temporaries. If a future backend cannot implement the required provider transform directly, that backend-specific wrapper should either expose a clear fallback contract or decline before side effects.
Consider the logical operation:
Assume A and B are column-major-compatible and C is row-major-like:
The mdspan staging descriptors describe the logical operands. Because C is row-major-like, the wrapper first chooses the GEMM rewrite:
It then derives provider-ready operands for the rewritten call:
C_write is the untransformed writable matrix seen by BLAS, so it represents logical C^T.
The provider call is therefore:
The intermediate staging layer is what made this rewrite possible: it knew that the output was row-major-like before constructing the final BLAS call. The final BlasReadableMatrix and BlasWritableMatrix operands alone are not enough to decide whether to swap operands, apply output-normalization transforms, or allocate temporaries.
If A were row-major-like instead, its stage would derive a BlasReadableMatrix whose rows and cols describe the provider's untransformed A^T storage, and whose transform is ‘'T’for a logical input. If composing storage layout, accessor conjugation, and any internal operation rewrite producesconjugate`, the generic direct wrapper declines. Provider-specific wrappers may add a fast path for OpenBLAS-style conjugate-no-transpose support later.
There are two wrapper contracts:
Output/update temporaries are stricter. They require explicit copy-back policy, aliasing checks before side effects, and diagnostics that report the materialization.
LAPACK workspace is not operand materialization. Many LAPACK operations need work arrays whose size is chosen by a workspace query or by the wrapper's algorithm policy. A direct LAPACK operation wrapper may allocate and free those work arrays internally while still being a direct no-copy wrapper for its matrix/vector operands. The distinction is:
For example, a direct self_adjoint_eigh wrapper may query and allocate work/rwork arrays before calling syev or heev. It may not silently copy a row-major matrix into a column-major temporary and copy eigenvectors back unless the wrapper contract says that operand materialization is allowed.
Example prepared input shape:
The prepared wrapper points readable.data at either the original storage or scratch storage.
The direct GEMM and GEMV wrappers are wired below descriptor-level backend implementations. Fixed-output tensor overloads derive the default selector from tensor storage, normalize the existing operands to MdspecLike refinements, and enter operation-tag dispatch. The selected backend acquires host mdspan leases and calls the direct wrappers. Bare-mdspan convenience overloads pass const-normalized inputs through the same descriptor boundary. The remaining axes of progress are:
The implemented GEMM/GEMV slices prove the dispatcher, backend selector, type acceptance CPO, runtime decline path, matrix/vector descriptors, and CPU fallback shape without mixing in LAPACK workspace or overwrite semantics. Storage-derived default backend lists are implemented in the tensor front end. The mdspan implementation functions themselves are lower-level Uni20 module APIs rather than operation-tag dispatch customization points.
The direct wrapper layer still exposes functions that are easy to test:
try_gemm(...) and try_gemv(...) return unsupported_layout when an operand cannot be staged as a direct no-copy BLAS view and unsupported_transform when the provider cannot express the required transform. GEMV therefore declines a conjugating input vector, while a conjugating row-major matrix can lower to conjugate-transpose. Once staging succeeds, dimension consistency is a checked precondition: invalid operation parameters are logic errors and should abort through CHECK/CHECK_EQUAL, not report ordinary kernel non-acceptance.
GEMM shape validation precedes provider staging. An empty output succeeds without inspecting operand layouts. A zero inner extent is mathematically the update C = beta*C. The host BLAS backend declines that runtime instance before side effects and CpuReferenceBackend performs the update. cuBLAS accepts null zero-sized input buffers and applies the update directly.
Operation-tag dispatch wraps the direct GEMM leaf like this:
kernel_accepts_types(...) remains type-level: scalar family, rank, mutability, and host raw-addressability. Descriptor construction remains runtime: dynamic extents, strides, pointer values, and aliasing.
For backends that provide this static CPO, it is the authoritative type filter. The corresponding try_kernel(...) should not repeat the same long requires clause; it may use a local static_assert against the shared support predicate to catch accidental direct calls with rejected types.
For GEMM and GEMV, try_kernel(BlasBackend, operation, ...) is a thin adapter over the corresponding direct wrapper. The operation-tag layer does not duplicate mdspan stride analysis. If the direct wrapper declines because an instance cannot be represented as a no-copy BLAS call, the dispatch walk continues to the next backend, normally the CPU reference oracle.
Completed checkpoint:
try_mdspan_matrix_stage(...) and role-specific readable/writable lowering helpers.uni20::blas::gemm.gemm_op, BlasBackend, and CpuReferenceBackend dispatch tests.try_kernel(BlasBackend, gemm_op, ...) by calling uni20::linalg::blas::try_gemm(...).beta == 0 no-read semantics.gemv_op, BLAS and CPU backend adapters, Tensor forwarding, and clean fallback for conjugating vector accessors.Completed LAPACK checkpoint:
blas_int.sterf/steqr, geev, gees, hseqr, and trexc adapters.self_adjoint_eigh adapter over syev/heev; the preserving eigh value API explicitly copies to a column-major work matrix first.The next wrapper should be driven by an active algorithm or example. It must define overwrite, aliasing, workspace, and any operand-materialization policy before its provider call is formed.
This order keeps the generic dispatch path tested by the GEMM and GEMV leaf kernels while separating tensor default-backend policy from LAPACK workspace and update-operand semantics.
Descriptor tests:
layout_left, layout_right, and layout_stride rank-2 views.1.1 and no singleton axis can supply the provider unit-stride dimension.blas_int.needs_conjugation is derived from the accessor policy.needs_conjugation == false.uni20::conj(mdspan) derives needs_conjugation == true and composes into readable transforms.needs_conjugation is rejected by the direct GEMM wrapper.GEMM tests:
GEMV tests:
alpha == 0, beta == 0, and empty-inner-dimension no-read semantics;LAPACK tests:
float, double, uni20::complex<float>, and uni20::complex<double>.prepared_*.hpp namespace/header? 1.9.8