R1 Linear Algebra (Core Ops)

Release 1 ships a “foundation” linear algebra surface with stable Python endpoints and a single routing boundary (so future CPU/GPU/out-of-core optimizations can land behind the same entry points).

This guide is about what you can call and what guarantees exist, not which backend is currently fastest.

Minimal example

import pycauset as pc

A = pc.matrix(((4.0, 1.0), (2.0, 3.0)))
b = pc.vector((1.0, 0.0))

x = pc.solve(A, b)

Core endpoints in Release 1

The Release 1 linalg surface includes these function families:

Block matrices (Release 1 behavior)

Construction: pycauset.matrix(...) builds a block matrix when given a 2D grid where every element is matrix-like; mixed scalars + matrices raise TypeError.
Once block, always block: @, +, -, *, / return block-matrix results with partition refinement (union of boundaries via SubmatrixView tiling). No silent densify; unsupported views error deterministically.
Semi-lazy orchestration: outputs are thunked per block. Triggers: element access, dense conversion, persistence, and crossing the compute boundary. Non-triggers: shape/partition metadata, repr/str, get_block.
Deterministic matmul: fixed k order and metadata-only dtype fold per output block; accumulation dtype is chosen before evaluation.
Block-aware slicing: slices produce tiled SubmatrixView blocks without copying; errors if a required view cannot be represented.
Device routing: leaf ops still run through AutoSolver; complex blocks route to CPU today on CUDA builds (see Compute Architecture). IO prefetch/discard hooks are best-effort and traced separately.

Shape constraints

Dense matmul follows NxM rules: (m, k) @ (k, n) -> (m, n).
Many routines are square-only by definition (inverse/determinant and most eigen routines).

Block matrices obey the same shape rules at the block grid level; partition refinement enforces compatible dimensions (A.cols == B.rows after refinement).

See NxM Support Status and Square-only Assumptions.

Property-aware behavior (R1_PROPERTIES)

Some linalg endpoints consult A.properties:

solve:
returns b when is_identity=True (square only)
rejects is_zero=True
routes diagonal/triangular claims to solve_triangular
solve_triangular treats triangular/diagonal claims as gospel and does not truth-validate.

See R1 Properties for the rules and failure modes.

Practical example (property-driven solve)

import pycauset as pc

A = pc.identity(3)
A.properties["is_upper_triangular"] = True
b = pc.vector((1.0, 2.0, 3.0))

x = pc.solve(A, b)  # routes to solve_triangular under the hood