BLA&DSP Subprogram Specification

[Mitchelldscott]

This issue tracks the foundational task of defining the generic traits that will form the core API of our library. The goal is to create a set of abstractions that align with standard BLAS (Levels 1, 2, 3) and common DSP interfaces. These traits will allow downstream algorithms to be written generically, independent of the underlying data storage (e.g., static no_std arrays vs. dynamic ndarray-backed matrices).

This is the first critical step outlined in Towards a Zero-Cost, Trait-Based Math Abstraction Layer for Control and Signal Processing in Rust 🚀

---

Tasks & Acceptance Criteria

A pull request that resolves this issue should include the following defined traits in a new module (e.g., src/traits.rs):

• 1. Generic Storage Traits:

  • Define base traits to abstract over vector and matrix-like containers. This is essential for writing generic operations.

  /// A trait for types that can be treated as a vector.
  pub trait Vector<T> {
      fn len(&self) -> usize;
      // Maybe an iterator or `as_slice` method?
  }
  
  /// A trait for types that can be treated as a matrix.
  pub trait Matrix<T> {
      fn shape(&self) -> (usize, usize); // (rows, cols)
      // Methods for row/col access?
  }

• 2. BLAS Level 1 (Vector-Vector) Traits:

  • trait Axpy<T, V: Vector<T>>: For scaled vector addition, $y \leftarrow \alpha x + y$.
  • trait Dot<V: Vector<T>>: For the dot product, $s \leftarrow x^T y$.
  • trait Nrm2<T>: For the Euclidean ($L_2$) norm, $|x|_2$.
  • trait Asum<T>: For the sum of absolute values ($L_1$ norm), $|x|_1$.
  • trait Scal<T>: For in-place vector scaling, $x \leftarrow \alpha x$.

• 3. BLAS Level 2 (Matrix-Vector) Traits:

  • trait Gemv<T, V: Vector<T>>: For general matrix-vector multiplication, $y \leftarrow \alpha A x + \beta y$.

• 4. BLAS Level 3 (Matrix-Matrix) Traits:

  • trait Gemm<T, M: Matrix<T>>: For general matrix-matrix multiplication, $C \leftarrow \alpha A B + \beta C$.

• 5. Core DSP Traits:

  • trait Conv<K: Vector<T>>: For 1D convolution.
  • trait Fft: For the Fast Fourier Transform.

---

Design Considerations & Open Questions

• Naming Convention: Should we stick to the traditional, terse BLAS names (gemm, axpy) for recognizability, or adopt more verbose, Rust-style names (general_matrix_multiply)?
  • Proposal: Stick to BLAS names. They are a well-understood standard in the domain.
• Scalar Generics: The traits should be generic over the scalar type T (e.g., f32, f64). What bounds should we use? T: num_traits::Float seems like a reasonable starting point.
• Mutability and Ownership: How should inputs and outputs be handled? Many BLAS operations are mutable. The method signatures should reflect this clearly. For example, axpy modifies its y vector, so it should probably be a mutable receiver or argument.

  // Example signature for AXPY
  fn axpy(&self, alpha: T, x: &impl Vector<T>, y: &mut impl Vector<T>);

• Error Handling: How should dimension mismatches be handled?
  • Proposal: For performance-critical code, mismatched dimensions are a logic error. We should panic in debug builds and potentially have unchecked versions for release builds. This avoids the overhead of returning a Result. Let's discuss this.
• Documentation: Each trait and method must include clear rustdoc comments explaining the operation it represents, including the mathematical formula in LaTeX.

[Mitchelldscott]

"However, converting arbitrary linear algebra expressions into an efficient sequence of well-matched calls to BLAS and LAPACK routines is non-trivial [8], [9]; manual conversion can be laborious and error-prone, and requires good understanding of the intricacies of BLAS and LAPACK, including various trade-offs across available routines and storage formats. A further downside of directly using BLAS/LAPACK routines is that the resultant source code is quite verbose, has little similarity to the original mathematical expressions, involves keeping track of many supporting variables, and requires manual memory management."[4]

Arm math performance libraries
Neon Programmers guide
CMSIS-DSP

Level 1

Operation	BLAS Name	Description (Formula)	Function Signature	Supported Data Types
Dot Product	dot	w = ∑ (xᵢ × yᵢ)	dot(N, x, x_stride, y, y_stride) -> w	single, double, complex (conjugated), complex (unconjugated)
AXPY	axpy	yᵢ = a × xᵢ + yᵢ	axpy(N, a, x, x_stride, y, y_stride)	single, double, complex
Vector Copy	copy	yᵢ = xᵢ	copy(N, x, x_stride, y, y_stride)	single, double, complex
Vector Swap	swap	yᵢ ↔ xᵢ	swap(N, x, x_stride, y, y_stride)	single, double, complex
Euclidean Norm	nrm2	sqrt(∑ xᵢ²)	nrm2(N, x, x_stride) -> w	(single, single), (double, double), (complex, single)
Sum of Magnitudes	asum	*∑	xᵢ	asum(N, x, x_stride) -> w
Vector Scaling	scal	x = a × x	scal(N, a, x, x_stride)	(single, single), (double, double), (single, complex), (complex, complex)
Index of Max	iamax	i such that xᵢ = max	iamax(N, x, x_stride) -> int	single, double, complex

Notes

1. For asum with complex numbers, the calculation is sum(abs(x_i.re) + abs(x_i.im)).
2. For iamax with complex numbers, the comparison is based on abs(x_j.re) + abs(x_j.im). The function returns the index of the first occurrence of the maximum value.

Level 2:

"Unless otherwise stated it is assumed that matrices are stored conventionally in a two-dimensional array with matrix element cij stored in array element A(I, J)." [2]

Category	BLAS Name	Description / Operation	Function Signature (Fortran-style)
General Matrix-Vector Products	gemv	General matrix-vector product: 1. y ← α·A·x + β·y if TRANS='N' 2. y ← α·Aᵀ·x + β·y if TRANS='T' 3. y ← α·Aᴴ·x + β·y if TRANS='C'	GEMV(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
	gbmv	General band matrix-vector product (with KL subdiagonals, KU superdiagonals)	GBMV(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
Symmetric / Hermitian Products	symv	Symmetric matrix-vector product:y ← α·A·x + β·y	SYMV(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
	hemv	Hermitian matrix-vector product:y ← α·A·x + β·y	HEMV(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
	spmv	Symmetric matrix-vector product (packed storage)	SPMV(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
	hpmv	Hermitian matrix-vector product (packed storage)	HPMV(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
	sbmv	Symmetric band matrix-vector product	SBMV(UPLO, N, K, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
	hbmv	Hermitian band matrix-vector product	HBMV(UPLO, N, K, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
Triangular Matrix-Vector Products	trmv	Triangular matrix-vector product:• x ← A·x if TRANS='N'• x ← Aᵀ·x if TRANS='T'• x ← Aᴴ·x if TRANS='C'	TRMV(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
	tpmv	Triangular matrix-vector product (packed storage)	TPMV(UPLO, TRANS, DIAG, N, AP, X, INCX)
	tbmv	Triangular band matrix-vector product	TBMV(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)

Notes

• TRANS parameter:
  • 'N': No transpose.
  • 'T': Transpose.
  • 'C': Conjugate transpose.
• UPLO parameter: Specifies whether matrix is upper or lower triangle ('U' or 'L').
• DIAG parameter: Specifies whether matrix has unit diagonal ('U') or non-unit diagonal ('N').
• Packed storage (AP): Only the relevant triangle stored as a contiguous array.

Level 3:

Category	BLAS Name	Description / Operation	Function Signature (Fortran-style)
General Matrix-Matrix Product	gemm	General matrix–matrix product: C ← α·op(A)·op(B) + β·C op(X) depends on TRANS:• 'N': X • 'T': Xᵀ • 'C': Xᴴ	GEMM(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Symmetric / Hermitian Matrix-Matrix Product	symm	Matrix product with symmetric matrix A: C ← α·A·B + β·C if SIDE='L' C ← α·B·A + β·C if SIDE='R'	SYMM(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
	hemm	Matrix product with Hermitian matrix A (complex): C ← α·A·B + β·C or C ← α·B·A + β·C	HEMM(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Rank-k Update (Symmetric)	syrk	Symmetric rank-k update: C ← α·A·Aᵀ + β·C if TRANS='N'C ← α·Aᵀ·A + β·C if TRANS='T'	SYRK(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
Rank-k Update (Hermitian)	herk	Hermitian rank-k update: C ← α·A·Aᴴ + β·C if TRANS='N'C ← α·Aᴴ·A + β·C if TRANS='C'	HERK(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)

Notes

• GEMM TRANSA / TRANSB:
  • 'N': No transpose
  • 'T': Transpose
  • 'C': Conjugate transpose (Hermitian)
• SYMM / HEMM SIDE:
  • 'L': Multiply A on the left
  • 'R': Multiply A on the right
• SYRK / HERK: For real data, 'T' and 'C' mean the same (transpose). For complex data:
  • In SYRK, TRANS='C' is not allowed.
  • In HERK, TRANS='T' is not allowed.
• UPLO:
  • 'U': Upper triangle of C is stored/updated.
  • 'L': Lower triangle.

[1] C. L. Lawson, R. J. Hanson, D. R. Kincaid, and F. T. Krogh, “Basic Linear Algebra Subprograms for Fortran Usage,” ACM Trans. Math. Softw., vol. 5, no. 3, pp. 308–323, Sep. 1979, doi: 10.1145/355841.355847.
[2] J. J. Dongarra, J. Du Croz, S. Hammarling, and R. J. Hanson, “An extended set of FORTRAN basic linear algebra subprograms,” ACM Trans. Math. Softw., vol. 14, no. 1, pp. 1–17, Mar. 1988, doi: 10.1145/42288.42291.
[3] J. J. Dongarra, J. Du Croz, S. Hammarling, and I. S. Duff, “A set of level 3 basic linear algebra subprograms,” ACM Trans. Math. Softw., vol. 16, no. 1, pp. 1–17, Mar. 1990, doi: 10.1145/77626.79170.
[4] C. Sanderson and R. Curtin, “Armadillo: an Efficient Framework for Numerical Linear Algebra,” in 2025 17th International Conference on Computer and Automation Engineering (ICCAE), Perth, Australia: IEEE, Mar. 2025, pp. 303–307. doi: 10.1109/ICCAE64891.2025.10980539.
[5] F. G. Van Zee et al., “The BLIS Framework: Experiments in Portability,” ACM Trans. Math. Softw., vol. 42, no. 2, pp. 1–19, Jun. 2016, doi: 10.1145/2755561.
[6] G. Frison, D. Kouzoupis, T. Sartor, A. Zanelli, and M. Diehl, “BLASFEO: basic linear algebra subroutines for embedded optimization,” ACM Trans. Math. Softw., vol. 44, no. 4, pp. 1–30, Dec. 2018, doi: 10.1145/3210754.

[Mitchelldscott]

Category	Method	Key Math Operations	BLAS	DSP
Classical	PID Control	Filtering, convolution, gain scaling	axpy	conv
	Root Locus	Polynomial root-finding, eigenvalue computation	eigen	poly_eval
	Bode/Nyquist Plots	Frequency response via FFT, magnitude/phase calc	norm2	fft
	Lead/Lag Compensation	Polynomial ratio manipulation, convolution	–	conv, poly_eval
Modern	LQR	Solve Riccati equation (matrix ops)	gemm, solve_triangular, lyapunov	–
	Kalman Filter	Matrix updates, gain computation, covariance	gemm, gemv, solve	–
	Observer Design (Luenberger)	Eigenvalue placement, matrix multiplication	eigen, gemm	–
	MPC	Quadratic programming, constraints	gemm, factorization, optimization	–
Robust	H∞ Control	Norm minimization, Riccati-like equations	norm2, gemm, solve	–
	µ-Synthesis	Iterative LMI solving, eigenvalue computations	eigen, factorizations	–
	Sensitivity Analysis	Frequency response, convolution	–	fft, conv
Nonlinear	Feedback Linearization	Symbolic differentiation, Jacobian/Hessian ops	gemm, gemv	–
	Sliding Mode Control	Sign functions, Lyapunov candidate evaluation	norm2, axpy	–
	Backstepping	Recursive Lyapunov equations, derivatives	gemm, solve	–
	Adaptive Control	Parameter update laws (gradient descent)	axpy, dot	–
System Identification	ARX/ARMAX Models	Least squares, convolution, regression	gemm, gemv	conv
	Subspace ID	SVD, eigen decomposition, projections	svd, eigen, gemm	–
	Maximum Likelihood Est.	Optimization, gradient/Hessian	gemm, grad	–
	State-Space Estimation	Matrix factorization, solve	factorizations, solve	–

[Mitchelldscott]

Here is the comprehensive guide, written from the perspective of a Computer Architect, detailing the physical cost of mathematical abstractions.

---

The Silicon Cost of Mathematics: A Hardware-First Taxonomy

Abstract:

Mathematical theory treats operations as abstract relationships between entities (e.g., $a + b = c$). Computer architecture treats operations as physical electrical signals traversing silicon pathways. This document maps theoretical mathematical structures to their concrete hardware implementations, organizing operations by "Cycle Cost"—the currency of computing.

---

Tier 0: The "Wire" Tier (Bitwise Logic)

Theoretical Equivalent: Identity, Negation, Set Membership.

Hardware Reality: 0–1 Cycles.

In the physical world, some mathematical operations require almost no logic gates to execute. They are matters of routing wires or flipping a single switch.

• The Operations: abs, copysign, signum, is_nan.

• The Hardware:

  • Register Renaming: Modern CPUs can sometimes handle moves and sign flips without even using an execution port.

  • Bit Masking: abs(f32) is not arithmetic. It is simply AND-ing the register with 0x7FFFFFFF to force the sign bit to 0.

• The Theoretical Map:

  • Math views |x| as a function of distance.

  • Hardware views |x| as a wire cut. This is why abs is distinct from "calculation."

---

Tier 1: The "ALU" Tier (Integer Arithmetic)

Theoretical Equivalent: Abelian Groups (Integers $\mathbb{Z}$), Modular Arithmetic.

Hardware Reality: 1 Cycle.

This is the domain of the Arithmetic Logic Unit (ALU). These operations are discrete, exact, and blistering fast. They represent the "native language" of the processor.

• The Operations: Integer add, sub, bitwise shl/shr, xor, and, or.

• The Hardware:

  • Carry-Lookahead Adders: Parallel circuits that calculate sums almost instantly.

  • Barrel Shifters: A grid of transistors that shifts bits left or right in a single clock tick.

• The Theoretical Map:

  • Theory defines $\mathbb{Z}$ (Integers) as a Ring.

  • Hardware implements this Ring via modulo $2^{32}$ or $2^{64}$.

  • Note: Integer multiplication (mul) used to be Tier 2, but modern silicon multipliers are so wide they are now effectively Tier 1 (often 3 cycles latency, but 1 cycle throughput).

---

Tier 2: The "Pipeline" Tier (Floating Point Basics)

Theoretical Equivalent: Fields ($\mathbb{R}$), Linear Combinations.

Hardware Reality: 3–5 Cycles (Latency), 1 Cycle (Throughput).

Here we enter the Floating Point Unit (FPU). Unlike integers, floats require a multi-stage assembly line: Unpack $\to$ Align Exponents $\to$ Add Significands $\to$ Normalize $\to$ Round $\to$ Repack.

• The Operations: Float add, sub, mul, fma (Fused Multiply-Add).

• The Hardware:

  • Pipelining: The FPU is like a factory line. It takes ~4 cycles to finish one addition, but because it is pipelined, you can shove a new addition in every single cycle.

  • FMA Units: Specialized hardware that computes $a \times b + c$ in a single step to preserve precision and save time.

• The Theoretical Map:

  • Math sees addition and multiplication as separate axioms of a Field.

  • Hardware often merges them (FMA) for efficiency.

  • Architect's Note: This is where associativity breaks. $(a+b)+c \neq a+(b+c)$ because the "Rounding" stage of the pipeline drops information.

---

Tier 3: The "Stall" Tier (Iterative Convergence)

Theoretical Equivalent: Multiplicative Inverses, Roots.

Hardware Reality: 10–25 Cycles. Non-Pipelined (Blocks Execution).

This is where the illusion of "instant math" breaks. Division and Square Roots are not calculated directly; they are solved. The hardware must guess the answer and refine it over several loops.

• The Operations: div ($x/y$), sqrt ($\sqrt{x}$), rem (remainder).

• The Hardware:

  • Iterative Solvers: Uses algorithms like Goldschmidt or Newton-Raphson baked into the silicon.

  • Pipeline Hazard: Because the unit is busy looping on one number, it often cannot accept new work. It "stalls" the pipeline, creating a traffic jam of instructions behind it.

• The Theoretical Map:

  • In Math, $x \cdot (1/x) = 1$ is a basic property.

  • In Hardware, mul is a Sprinter; div is a Heavy Lifter. They are not peers.

---

Tier 4: The "Microcode" Tier (Transcendentals)

Theoretical Equivalent: Real Analysis, Infinite Series, Geometry.

Hardware Reality: 50–150+ Cycles.

These operations are too complex for a single logic gate circuit. When you call these, the CPU effectively pauses your program to run a tiny, internal subroutine stored in Read-Only Memory (Microcode) or traps to a software library (libm).

• The Operations: sin, cos, exp, log, pow.

• The Hardware:

  • CORDIC / Polynomial Approximation: The CPU breaks the request down into a long sequence of adds, multiplies, and table lookups (LUTs).

  • Range Reduction: Calculating $\sin(100000.5)$ requires reducing the argument to the range $[-\pi, \pi]$ first, which is incredibly expensive and error-prone.

• The Theoretical Map:

  • Math defines these via infinite series (Taylor Series).

  • Hardware approximates them. The "Transcendental" nature of these numbers means exact representation is physically impossible; we are paying for bits of precision.

---

Summary Matrix: The Architect's Cheat Sheet

Tier	Mathematical Structure	Hardware Implementation	Cycle Cost (Approx)	Pipeline Status
0	Identity / Logic	Wires / Muxes	0–0.5	Instant
1	Ring (Integer)	ALU (Gates)	1	Full Speed
2	Field (Basic Float)	FPU (Pipeline)	4	Full Speed
3	Inverse / Roots	Iterative Logic	20	STALL
4	Analysis (Trig/Exp)	Microcode / ROM	100+	BLOCKING

Next Step

For a high-performance system, a "Fast Math" trait library should expose Tier 0, 1, and 2 as default methods, but gate Tier 3 and 4 behind explicit ComplexMath or Approximation traits to force the developer to acknowledge the cost.

[Mitchelldscott]

Project: math-rs

Hardware-Aware Mathematical Abstractions for Embedded Control

Abstract

safety_math is a no_std Rust library that maps abstract mathematical theory to the physical reality of processor pipelines. By segregating operations into traits based on their "Cycle Cost" (Bitwise, Integer, Pipeline, Stall, and Microcode), the library forces architectural awareness at the type level. This allows control algorithms (like PID or Kalman Filters) to be written generically, while the underlying hardware implementation (hardware FPU vs. software interpolation vs. DSP accelerator) is swapped via zero-cost marker types. This ensures that a control loop intended for a 10kHz cycle time cannot accidentally invoke a 100-cycle software emulation of a transcendental function without explicit developer opt-in.

1. The Tiered Trait Hierarchy

We define traits not just by mathematical function, but by hardware behavior. This prevents "Stall" tier operations (Division/Sqrt) from sneaking into "Pipeline" tier loops.

// Tier 1: Single Cycle / Pipeline Friendly
// Maps to: ALU ops, basic FMA.
// Cost: 1-4 cycles, 1 cycle throughput.
pub trait LinearOps<T>: Copy {
    fn add(self, other: Self) -> Self;
    fn mul(self, other: Self) -> Self;
    fn scale(self, factor: T) -> Self;
}

// Tier 3: The "Stall" Tier
// Maps to: Iterative hardware units (Div, Sqrt).
// Cost: 10-30 cycles (Pipeline bubbles).
// Note: We return Result to force handling of Domain Errors (div by zero).
pub trait NonLinearOps<T>: LinearOps<T> {
    fn div(self, other: Self) -> Result<Self, MathError>;
    fn sqrt(self) -> Result<Self, MathError>;
}

// Tier 4: The "Microcode" Tier
// Maps to: libm calls, software emulation, polynomial approximation.
// Cost: 50-150+ cycles (Function call overhead).
pub trait TranscendentalOps<T>: NonLinearOps<T> {
    fn sin(self) -> Self;
    fn cos(self) -> Self;
    fn atan2(self, other: Self) -> Self;
}

2. Hardware Marker Types (The Backends)

These are zero-sized types (ZSTs) that serve as compilation configurations. They implement the traits using the specific assembly or intrinsic instructions available on the target core.

/// Marker for a standard ARM Cortex-M4F (Hardware FPU present).
/// Implements LinearOps via VADD/VMUL (Hardware).
/// Implements NonLinearOps via VDIV/VSQRT (Hardware, but stalling).
pub struct CortexM4F;

/// Marker for a low-power RISC-V core (RV32I) without F Extension.
/// Implements LinearOps via software library (__aeabi_fadd).
/// Does NOT implement TranscendentalOps to prevent accidental bloat.
pub struct SoftFloatRV32;

/// Marker for a system with a dedicated CORDIC coprocessor.
/// Implements TranscendentalOps via memory-mapped IO writes to the accelerator.
pub struct CordicAccelerator;

// Example Implementation
impl LinearOps<f32> for CortexM4F {
    #[inline(always)]
    fn add(a: f32, b: f32) -> f32 { a + b } // Compiles to vadd.f32
    #[inline(always)]
    fn mul(a: f32, b: f32) -> f32 { a * b } // Compiles to vmul.f32
}

3. The Control System Interface

Higher-level algorithms use the Marker Type as a default generic. This allows the logic to remain pure, while the implementation can be hot-swapped for testing or porting.

/// A PID Controller agnostic of the underlying hardware.
/// 
/// Generics:
/// T: The numeric type (f32, f64, fixed_point)
/// Backend: The hardware implementation (defaults to the target's best guess)
pub struct PidController<T, Backend = CortexM4F> {
    kp: T,
    ki: T,
    kd: T,
    _marker: PhantomData<Backend>,
}

impl<T, Backend> PidController<T, Backend>
where 
    T: Copy,
    Backend: LinearOps<T> // Constraint: PID only needs Add/Sub/Mul (Tier 1/2)
{
    pub fn update(&mut self, error: T) -> T {
        // This code is guaranteed to run in O(1) pipeline throughput
        // because Backend is bound by LinearOps, not TranscendentalOps.
        let p_term = Backend::mul(self.kp, error);
        // ... calculation ...
        p_term
    }
}

4. Architectural Safety Checks

We can use this system to break the build if a developer tries to use expensive operations in critical sections.

// A strict backend for a 100kHz control loop.
pub struct HighFreqLoop;

impl LinearOps<f32> for HighFreqLoop { /* ... */ }

// We deliberately DO NOT implement NonLinearOps or TranscendentalOps.
// If a user tries to normalize a vector (requiring sqrt/div) inside
// this loop, the code will fail to compile.