|
Uni20
Asynchronous tensor and runtime library.
|
Status: implemented first vertical slices plus forward plan. The mdspan linalg dispatch layer now covers GEMM/GEMV, accessor-respecting CPU fallbacks, full and partial sum reductions, full inner-product and stable-norm reductions, Tensor output-shape preparation, matrix products, explicit copy/materialization, matrix initialization and exponential front ends, and the LAPACK projected eigensystem/Schur operations used by native Krylov. Broader BLAS/LAPACK coverage and domain-aware prepared operands remain design work.
This note records the implemented mdspan wrappers over the existing Uni20 BLAS/LAPACK backend layer and the next steps for extending them. These wrappers are the first real consumer of the kernel-dispatch design in `../architecture/kernel_dispatch.md`.
The static CPO name in this note is kernel_accepts_types(...). It returns a KernelTypeAcceptance value (no, maybe, or yes) for the C++ argument types. It does not inspect extents, strides, pointer values, library availability, or device state. Those checks belong to the runtime attempt, try_kernel(...).
Related notes:
The tensor/linalg frontend should accept tensor views rather than a Krylov-only matrix class or vendor-specific pointer API. Once it has selected a backend list, fixed-operation dispatch should use normalized mdspecs. The selected backend lowers those descriptors to execution-domain mdspans for its lower-level implementation.
The intended contract is:
For LAPACK, this means the type-level check verifies host-accessible scalar storage with a supported scalar type and the required rank/mutability. The runtime check verifies details such as extents fitting blas_int and matrix strides being representable by LAPACK.
Matrix operands for this layer are rank-2 views. Vector operands such as eigenvalue arrays, Householder coefficients, pivots, and work/result vectors are rank-1 views. Rank is therefore an immediate compile-time filter in the linalg concepts, not a runtime LAPACK-layout check.
Uni20 already has the pieces to build on:
uni20::StridedMdspanLike and uni20::MutableStridedMdspanLike in src/uni20/mdspan/concepts.hpp.src/uni20/backend/lapack/lapack.hpp.src/uni20/backend/lapack/reference and src/uni20/backend/lapack/mplapack.uni20::DenseMatrix<T, LayoutPolicy> Tensor alias in src/uni20/tensor/tensor.hpp.src/uni20/krylov/dense_subspace.hpp.The abandoned ../cytnx-mdspan prototype is useful for operation-tag examples and diagnostic formatting, but Uni20 should not copy its std::variant dispatch funnel for the first pass. Uni20's current draft design uses an ordered backend list with a type-level capability CPO and a runtime attempt CPO.
uni20::lapack::* wrappers.uni20::lapack for supported binary128 types.try_kernel(...) wrappers as a surprise fallback. Prepared-view helpers may materialize input-only operands into explicit scratch matrices, and higher-level tensor operations may own output copy-back policy deliberately.std::variant backend selection yet. It can be layered on later as another backend-list entry if Python or plugin boundaries need it.The final layering should look like this:
kernel_accepts_types(...) / detected try_kernel(...).try_kernel(...) succeeds.MdspecLike refinements.prepare_output(...); a later backend may reuse or replace it after a decline.Descriptor normalization preserves rather than projects the operand. Future symmetry-aware, distributed, file-backed, and remote-storage descriptors must retain their sector metadata, distribution, communicator, file mapping, or resource identity through dispatch. Extracting a descriptor never implies a dense conversion, host transfer, or data acquisition.
Public bare-mdspan convenience calls may require an explicit backend selector because they cannot derive a default Uni20 backend stack. They pass a const-normalized input and mutable output directly through the same mdspec dispatch boundary. Resolved-mdspan functions below the backend are lower-level Uni20 module or provider interfaces and do not participate in the operation-tag dispatch walk.
The first implementation refines the existing Uni20 mdspan concepts instead of introducing unrelated concept names. Generic rank filters are not specific to linalg and live with the mdspan concepts. The linalg layer should then add scalar-family and host-addressability refinements.
Implemented generic mdspan refinement:
Candidate linalg refinements:
The exact names can change, but the split matters:
StridedMdspanLike: structural mdspan-like API with runtime strides.RankedStridedMdspanLike<View, Rank>: structural mdspan-like API plus a static rank requirement; this belongs in the generic mdspan layer and is implemented in src/uni20/mdspan/concepts.hpp.LinalgMatrixView and LinalgVectorView: rank-specific real-or-complex linalg operands.RealLinalgMatrixView, ComplexLinalgMatrixView, RealLinalgVectorView, and ComplexLinalgVectorView: longer names for operations that need a specific scalar family.HostRawAddressableRankedView: memory can be passed to a host pointer ABI. A pointer-like data_handle_type alone is not enough to prove this. Accessor semantics matter: a conjugating accessor can preserve the original pointer handle while changing the values returned by access(...). Provider-facing direct concepts should therefore distinguish "storage handle can be passed to
the ABI" from "the accessor is identity-like or otherwise explicitly lowered
into provider transform metadata". The default direct-memory case is an mdspan with stdex::default_accessor<T>. Other accessors are direct only when the wrapper recognizes their semantics and lowers them deliberately. The BLAS adapter currently treats Uni20's C++26-style conjugated_accessor as lowerable for readable inputs because it maps to MatrixTransform; arbitrary transform accessors still require materialization or a generic path.BlasScalarRankedView and LapackScalarRankedView: element type is exactly supported by the configured BLAS or LAPACK layer. These are backend refinements, not base linalg concepts. Keep both names because BLAS and LAPACK coverage may diverge for extension scalar types.accessor_applies_conjugation_v and mdspan_needs_conjugation_v live in src/uni20/mdspan/conjugate_accessor.hpp. The generic mdspan conj(...) helper returns a read-only conjugating accessor view for complex mdspans, cancels to the const original view when applied twice, and is a no-op for non-complex values while still returning a const identity mdspan view. A conj(Tensor) returns a tensor view whose mdspan accessor advertises the same trait, so the mdspan-to-BLAS adapter discovers conjugation from the accessor instead of from tensor-specific side metadata. The accessor follows the C++26 std::linalg::conjugated_accessor direction in WG21 P3050R3, but Uni20 keeps uni20::conj as the value-level customization point instead of adopting conj-if-needed terminology.real(x) and imag(x) are slice-like structural views with adjusted handles and strides, not proxy accessor adaptors. Writable proxy or semantic-transform accessors are special operation-specific cases that need an explicit assignment law and backend lowering. The general policy is in docs/tensor/dispatch_and_view_semantics_draft.md.The implementation can add named mutable real/complex aliases such as MutableRealLinalgVectorView when call sites become clearer, but the important policy is that the short names are real-or-complex and the longer names state a specific scalar family.
The compile-time check should reject device accessors, arbitrary transform accessors, unsupported scalar types, rank mismatches, and const outputs before the function body is instantiated. A transform accessor can still be eligible when it exposes raw host storage and declares a BLAS-lowerable operation such as accessor-level conjugation through a trait.
New Uni20 linalg and kernel APIs should put mutable outputs first. This applies to mdspan leaf wrappers, tensor-facing linalg wrappers, and backend runtime-attempt overloads.
Examples:
For update operations, the updated object is both input and output, but it still appears first. API tags and backend selectors appear before outputs; options that are ordinary configuration values follow the required operands unless they must prefix a parameter pack. Older draft documents may still contain BLAS/LAPACK-style output-last examples; new code should use the prefix-tag, output-first convention unless an external ABI boundary forces another order.
The first concrete implementation target is a new mdspan wrapper layer above the existing BLAS/LAPACK backend facades. The detailed plan lives in `blas_lapack_wrappers.md`.
The important split is:
src/uni20/backend/blas and src/uni20/backend/lapack remain provider and ABI facades. They own vendor detection, LP64/ILP64 details, raw declarations, checked LAPACK status handling, and MPLAPACK provider overloads.src/uni20/linalg mdspan wrapper layer owns view interpretation: ranked mdspan concepts, mdspan-storage-first MdspanMatrixStage descriptors, derived provider-ready BlasReadableMatrix and BlasWritableMatrix objects, transform helpers, direct no-copy runtime acceptance or decline, prepared input temporaries, and operation-specific rewrites such as row-major GEMM output handling.This keeps dense linalg policy out of the raw provider wrappers while giving the kernel-dispatch design a real, testable leaf-kernel boundary.
The staging object should store logical mdspan extents, the non-unit stride, and the integer mdspan axis whose stride is 1. It should not store a BLAS transpose flag. Helpers derive the provider-ready BLAS rows, columns, leading dimension, and readable-input transform from those storage facts plus accessor metadata such as conjugation.
The mdspan linalg layer supplies lower-level implementations used by the operation-tag model. The fixed-operand operation families select from tensor policy, normalize each fixed operand once, and dispatch those descriptors to the candidate backends. This includes copy, GEMM, GEMV, elementwise transforms, reductions, matrix initialization and exponentiation, and the fixed-output LAPACK eigensystem, Schur, and SVD operations. Each selected host backend acquires access and calls its ordinary mdspan implementation. GEMV additionally lowers rank-one BLAS increments before LAPACK workspace policy enters the picture.
Backend CPOs put the backend value first, then the operation tag, then ordinary reference parameters matching the public call order. A separate state_type<State> / backend_type<Backend> tag layer is not needed for the first pass.
try_kernel(...) assumes dispatch has already obtained a non-no result from kernel_accepts_types(...); it does not repeat that type predicate.
The tensor frontend normalizes fixed operands before calling the generic dispatcher:
Here c, a, and b model the appropriate tensor-view concepts, while the operation-tag customization receives their normalized mdspecs. Direct calls such as blas::try_gemm(...) operate below dispatch on resolved mdspans.
The first LAPACK wrapper can then use the same pattern with richer operand rules:
The workspace query and allocation in this example are LAPACK work storage, not operand materialization. The wrapper still counts as direct if it does not copy, pack, transpose, or conjugate the user-visible matrix/vector operands.
kernel_accepts_types(...) should stay type-level. It should not inspect runtime strides, sizes, pointer values, or backend state. The dispatcher inspects its result type through unevaluated decltype and std::declval expressions. The runtime attempt does the value checks and returns a specific non-success KernelAttempt if the selected backend cannot represent this particular view.
The acceptance constants encode the decision in their types, so the dispatcher does not need synthetic objects, a tuple, or an explicit type-pack token:
kernel_accepts_types(...) owns type eligibility and may be narrowly constrained. If it is not callable for the exact argument types, that backend is a hard no; a callable try_kernel(...) does not change the result. try_kernel(...) can therefore remain unconstrained, and the dispatcher does not instantiate it for rejected types. This single-backend query remains a detail helper. Public code uses probe_dispatch_kernel(backends, op, args...), which returns yes if any candidate is yes, otherwise maybe if any is maybe, otherwise no.
The standard dense symmetric/Hermitian eigensolver is now implemented. It exercises real and complex scalar paths and shows the distinction between type eligibility and runtime view representability.
Use an output-first, in-place work-matrix signature:
matrix_work is mutable because LAPACK overwrites its input. If compute_vectors is true, matrix_work contains the eigenvectors on return. If compute_vectors is false, matrix_work should be treated as destroyed workspace. The higher-level eigh(matrix) value API has a preserving lvalue form that copies into an explicit column-major temporary and a consuming owning rvalue form that may reuse the input allocation:
Reuse is valid because the matrix is a single update operand: LAPACK's A argument is both the input matrix and, for ‘jobz = 'V’, the eigenvector output. For directly addressable column-major storage, the consuming wrapper transfers the allocation unchanged. For square row-major storage, it adopts the same allocation under a column-major{1, LDA}mapping, exchangesUpperand Lower`, and conjugates complex eigenvectors after LAPACK. This preserves a padded leading dimension and any unused storage tail. A mapping without a unit-stride axis is copied into column-major work storage instead.
This should not be represented as arbitrary aliasing between independent parameters. If a higher-level API has separate input and eigenvectors arguments, then it should either:
input into a work matrix before calling LAPACK.Partial overlap must be rejected. The eigenvalue vector must not alias the matrix work storage or LAPACK workspace.
The scalar relation is:
So ssyev/dsyev/Rsyev cover real matrices, and cheev/zheev/Cheev cover complex Hermitian matrices.
The operation-specific concept can be written in terms of the generic rank and scalar-family concepts:
For the LAPACK backend this is not quite enough, because LAPACK also needs raw host-addressable storage:
Then the backend type check is concise:
This rejects, at compile time:
The strict first LAPACK implementation should require direct writable matrix representation with unit stride along mdspan axis 0:
In descriptor terms, this is unit_stride_axis == 0 and nonunit_stride == stride(1). stride(1) == 1 row-major-like matrices should initially decline. Real symmetric value-only solves can probably support row-major-as-transposed by flipping uplo, but complex Hermitian eigenvectors are easy to mishandle because row-major storage interpreted as column-major storage represents a transpose. Keeping the first wrapper strict avoids silently returning conjugated/transposed eigenvector data. A higher-level tensor wrapper can make an explicit column-major temporary when it wants to support arbitrary strided inputs.
The runtime attempt is schematic but should have this shape:
The actual helper names can differ, but the responsibilities should stay split:
kernel_accepts_types(...) says this backend can compile for these view types, returning maybe for direct LAPACK views because extents, strides, aliasing, and provider availability are runtime facts.try_blas_column_major_writable_stage(...) checks the particular extents, strides, and writable storage orientation.Matrix scratch before calling a direct LAPACK wrapper.uni20::lapack::syev/heev perform the provider-dispatched LAPACK call.Minimal tests:
float, double, uni20::complex<float>, and uni20::complex<double>.static_assert.1.unit_stride_axis == 1 row-major-like matrices until that case is explicitly implemented.||A v_i - lambda_i v_i|| / scale within the scalar-specific tolerance.Low-level try_kernel(...) functions return a specific non-success KernelAttempt for clean runtime decline before execution effects. Provisional preparation of an operation-declared replaceable output is permitted; input, fixed-output, and update-output mutation is not. Once a LAPACK routine has been called, a terminal positive INFO raises a structured LapackError; dispatch must not continue because update operands may already have been overwritten. A negative INFO is an unconditional invariant failure: it means the provider rejected arguments constructed by the checked Uni20 wrapper. Positive values explicitly documented as usable operation statuses remain operation-specific return values. Decline categories are stable enum values rather than backend-constructed strings.
Public convenience wrappers turn "all backends declined" into a structured KernelDispatchError. Its presentation report lists the operation and each backend's no/maybe/yes type acceptance, together with its structured runtime result or not eligible. This records the dispatch decision without requiring each backend to construct text.
Operand metadata such as extents, strides, scalar type, layout type, and accessor kind may be added to optional structured diagnostics later.
Direct mdspan GEMM proves the BLAS matrix descriptor, transform helpers, and row-major writable output rewrite. GEMV additionally proves rank-one descriptors and mixed-rank dispatch. The first LAPACK mdspan wrappers cover the operations needed by the active native Krylov projected problems. The remaining rows are candidates to add only when an algorithm or example needs them:
| Operation | LAPACK routines | Status/reason |
|---|---|---|
| LU conditioning | gecon | lu_rcond_op estimates one-/infinity-norm reciprocal conditioning from reusable factors and captured original norms. Native and configured MPLAPACK providers share the API; no CPU estimator fallback. |
| symmetric tridiagonal eigensystem | sterf, steqr | Implemented for Lanczos projected diagonalization. |
| real and complex Schur decomposition | gees, hseqr, trexc | Implemented for Arnoldi restarts and reordering. |
| nonsymmetric eigensystem | geev | Implemented for projected Ritz extraction and validation. |
| symmetric/Hermitian dense eigensystem | syev, heev, then syevd/heevd | syev/heev implemented with in-place and preserving value APIs; divide-and-conquer variants remain future work. |
| QR/LQ factorization | geqrf, orgqr, gelqf, orglq | Reduced real/complex QR and LQ are dispatched through destructive-workspace, preserving and consuming Tensor APIs. Complex Q uses ungqr/unglq; MPFR/MPC adds explicit-precision provider packing. |
| SVD and least squares | gesvd, gesdd, gesvdx, gelss, gelsd | Exact real and complex SVD is dispatched through gesvd. Tensor APIs cover singular values only, either factor separately, or both factors; preserving and consuming Async forms mirror them. Reduced consuming factors may reuse the input allocation through ‘JOBU='O’orJOBVT='O'. Densetruncated_svdis a Tensor policy layer over reduced exact SVD rather than a new mdspan kernel. [Full-column-rank least squares](least_squares.md) dispatches through reduced QR with recoverable diagnostics.svd_least_squaresadds rank-aware minimum-norm solves throughgelss`, with epsilon-based truncation, residuals and recoverable nonconvergence. Divide-and-conquer and pivoted-QR least-squares methods remain future coverage. |
Do not expose every backend wrapper through mdspan immediately. Add wrappers when a Uni20 algorithm or example needs them, and keep unused wrappers in a separate tested header if we want to preserve prototype coverage.
Implement the wrapper layer described in `blas_lapack_wrappers.md`.
The completed BLAS slices are:
MdspanMatrixStage helpers, derived BlasReadableMatrix/BlasWritableMatrix helpers, and transform helpers;uni20::blas::gemm, including row-major writable output rewrite tests;MdspanVectorStage and provider-ready vector operands;The completed LAPACK slice adds checked workspace conversion and strict direct column-major adapters for tridiagonal eigensystems, nonsymmetric eigensystems, Schur decomposition, Hessenberg Schur reduction, and Schur reordering. Runtime layout decline occurs before mutation; provider INFO failures are terminal.
This layer should call uni20::blas::* and uni20::lapack::*, never raw Fortran symbols, and should keep strict direct wrappers non-owning and no-copy. Prepared wrappers may later own input-only materialization by contract.
The generic pieces from ../architecture/kernel_dispatch.md are implemented for the mdspan dense-linalg layer:
backend_list<...>.kernel_accepts_types(backend const&, op const&, args&...) detection using private type-probe lvalues.try_dispatch_kernel(...) returning whether a runtime candidate ran.dispatch_kernel(...), which reports an exhausted list through ERROR.dynamic_dispatch_kernel(...) for Python and runtime-erased boundaries that must convert a type-level no into a runtime ERROR.Tests should use fake backends to verify:
Bare-mdspan convenience overloads adapt their operands to tensor views before generic dispatch:
Add operation-specific convenience functions only where the Tensor layer has real policy work to perform, such as allocating results, selecting storage backends, or resolving Tensor mdspans. Mutable outputs should appear first. For in-place LAPACK-style operations, the overwritten work matrix is a mutable output/update operand and should be grouped with the other output operands before read-only inputs.
This checkpoint is implemented. uni20::DenseMatrix<T, LayoutPolicy> is an owning rank-2 alias of the concrete Tensor class, column-major by default, and exposes its mdspan through the ordinary Tensor interface. The active Krylov projected eigensystem, Schur, reorder, matrix-set, and matrix-exponential paths now use matrix-level linalg front ends rather than raw pointer-level LAPACK calls. Explicit row/column layout conversion remains where a Krylov algorithm genuinely needs a differently ordered representation.
The Tensor front-end checkpoints now:
mdspec() before dispatch; the generic CPU path does not require stridedness, while BLAS lowering does.gemm_op/gemv and add_product updates from replaceable-output assign_product_op overwrites. gemm_op remains fixed even when beta is numerically zero; overwrite permission is carried by the operation type rather than inferred from a scalar value.conj(tensor) views and explicit eager copy and make_tensor operations. make_tensor(conj(input)) allocates first and then dispatches copy_op.Replaceable-output backend adapters call prepare_output(...) before resolving the writable mdspan. They may still return a clean decline after preparation; the next backend receives the prepared output and may reuse or replace it. Fixed-output GEMM and GEMV only validate their existing outputs and remain unchanged on decline.
The current BLAS assign_product_op adapter provides a stronger guarantee by resolving its input mdspans and probing BLAS with prospective output mapping/accessor metadata. The probe uses no usable data handle and performs no provider call. This avoids provisional host-output preparation, but that extra preflight is an implementation choice rather than a dispatch requirement.
This is where CUDA, MPI/block tensor placement, async scheduling, and temporary allocation policy enter. They should not be forced into the first LAPACK mdspan wrapper step.
Use three classes of tests:
mdspan_needs_conjugation_v == false.uni20::conj(mdspan) to verify accessor-derived conjugation metadata and double-conjugation cancellation.layout_left, layout_right, and layout_stride descriptors.1.1.float, double, uni20::complex<float>, uni20::complex<double>, and binary128 types when MPLAPACK is enabled.For kernel dispatch itself, use fake backends first. That keeps dispatch semantics independent of LAPACK availability and makes compile-time and runtime decline behavior easy to test.
1 before a direct LAPACK call is even a candidate.uni20::lapack.std::variant dispatch is deferred.try_kernel(...) wrappers do not hide layout copies. Higher linalg wrappers may use prepared-view helpers that own input scratch matrices or explicit output copy-back semantics.conj(...) remains a view. Explicit copy_op observes accessor and layout semantics; a future BLAS implementation may lower rank-two copy to provider omatcopy-style extensions without making conjugation eager.copy_op for transpose/conjugate combinations (omatcopy, omatcopy2, or equivalent)?CblasConjNoTrans or its Fortran-style R spelling, and where should fallback materialization live?input_temporaries, and which should require direct_only unless the caller opts into materialization?try_* and *_or_throw forms?uni20::linalg, uni20::kernel, and uni20::backend? 1.9.8