SparseBlas

In-place y = alpha · op(A) · x + beta · y, where op(A) is Aᵀ when transpose — the sparse dgemv. Per BLAS convention beta == 0.0 overwrites y without reading it, so it may arrive uninitialized, and alpha == 0.0 reduces to the beta scale.

Both directions walk the columns, which is what CSC stores. The non-transposed form accumulates x_j times column j into y; the transposed form dots column j with x to make y_j. Neither touches a structural zero.

A · x, or Aᵀ · x when transpose, into a fresh result — the restricted gemv with alpha = 1, beta = 0.

Solve op(T) · x = b in place, where T is the lower or upper triangle of the square a and op transposes when transpose — the sparse dtrsv. x holds the right-hand side on entry and the solution on return.

The solve existed inside SparseLu and nowhere else, so a caller holding a triangular SparseMatrix of their own had no way to solve against it. This is that routine, on the seam a host sparse BLAS would replace (MKL has mkl_sparse_d_trsv; it is one of the routines an inspector-executor interface covers).

Every direction walks columns, which is what CSC stores, and only the selected triangle is read — so the other one may hold anything. The two non-transposed forms push a finished unknown down or up its column; the transposed forms pull, dotting the column against the already-solved entries. That is the same four-way structure the dense cores have, for the same reason: a column of T is a row of Tᵀ.

Unlike the dense cores, this validates the diagonal. Triangular.kt leaves a singular dense triangle to produce infinities on the grounds that the caller knows what it passed; here a structurally missing diagonal entry is detectable for free while walking the column, so reporting it costs nothing and a silent NaN would be a worse trade.

Throws