test(blas): port the OpenBLAS BLAS-level checks, first batch

NOTICE

@@ -33,12 +33,19 @@ terms and their notices ship alongside the crate:
                                         column-pivoting QR, Schur, Hessenberg,
                                         Schur reordering, tridiagonal QR step,
                                         selected LAPACK EIG/LIN numerical validation tests
+  faer/COPYING.BLAS.NETLIB              credit for the reference BLAS test programs
+                                        whose BLAS-level checks the matmul,
+                                        rank-update and triangular-solve unit tests
+                                        adapt (docs/blas-test-inventory.md)
   faer/COPYING.SUITE_SPARSE.AMD.BSD     AMD ordering
   faer/COPYING.SUITE_SPARSE.COLAMD.BSD  COLAMD ordering
 
 The faer package's `license = "MIT AND BSD-3-Clause"` metadata includes the
 LAPACK numerical test adaptations. It does not relicense the other ported files;
-read the COPYING files for their actual terms.
+read the COPYING files for their actual terms. The netlib reference BLAS test
+programs ask for credit rather than imposing a license condition, and no Fortran
+is translated: `faer/COPYING.BLAS.NETLIB` carries that credit and the quoted
+notice, so the `license` expression stays as it is.
 
 The packages published from this repository are the temporary tensor4all fork
 names `t4a-faer` and `t4a-faer-traits`, so that downstream crates can depend on

PATCHSET.md

@@ -196,6 +196,7 @@ gate has changed and two kinds of adoption are now in the ledger.
 | Explicit sparse column ordering stops paying for COLAMD (#322) | `factorize_symbolic_lu_with_column_ordering` and `factorize_symbolic_qr_with_column_ordering` build their internal `MemBuffer` from `StackReq::or(colamd::order_scratch::<I>(m, n, A_nnz), ..)`, but COLAMD runs only in the `None` branch, so a caller that brings the ordering still allocates a COLAMD-sized workspace on every symbolic analysis (dyn-stack's `MemBuffer` takes it uninitialized, so the cost is the allocation, not a memset). Both functions include `order_scratch` only when `column_ordering.is_none()` now; the `None` request and the analysis request are unchanged. Both entry points are this fork's own (#239, `3889bd1`), so upstream faer does not carry the defect and no `upstream-reports/` entry is added. | KEEP | tensor4all/faer#322 | Allocation only, no numerical change: the two functions' permutations, symbolic structures and every downstream numeric path are untouched, and no benchmark was run. The issue's reproducer (n = 2000, nnz = 7990, identity ordering, `--release`, counting global allocator) goes from 332712 B largest allocation to 239928 B for simplicial LU and QR: the 332712 B was exactly `colamd::order_scratch`, and 239928 B is the analysis request (4 + 5 n-element arrays, the transpose arrays, and the largest simplicial arm), while the supernodal arm stays at 468080 B. It was attributed afterwards (tensor4all/faer#340, measured at `e756f13` and `b989876`, two runs each): it is one `Vec<usize>` of 58510 elements, allocated three times per symbolic analysis inside `cholesky::supernodal::ghost_factorize_supernodal_symbolic`; the 64000 B four-array helper only sizes a `MemBuffer` and never appears as an allocation of its own. The figure itself is #337's reported run of issue #322's counting program and was not re-measured for this change. `faer/tests/sparse_symbolic_colamd_workspace.rs` pins both directions on a 100 x 100 sample: with an explicit ordering the largest allocation stays below `order_scratch`, without one it reaches it; the first fails before the change and both pass after. `cargo test -p t4a-faer -j6` passes except `lq_no_alloc::lq_allocates_nothing_60x300`, which fails identically on `origin/main` (dense LQ, unrelated), plus `no_step_alloc`, which counts heap allocations, and `triangular_solve::aliased_rhs::staged_solve_is_bitwise_the_in_place_one`, which counts staged solve nodes through `stage::STAGED` (`faer/src/linalg/triangular_solve.rs:509,579`) and asserts that a 1 KiB stride was **not** staged (the message names the failure), not that anything allocated. Both pass standalone at this revision and both were seen to fail inside the whole-suite run; what makes them fail there was not identified. `no_step_alloc` is its own test binary, so another test binary cannot perturb its counter directly, which rules out the obvious explanation without supplying a better one. |
 | Complex `abs_max` drops a NaN in the real part (#320) | `ComplexField::simd_abs_max` takes the larger of the two part magnitudes with an unordered `re > im` select, and that select returns `im` when `re` is NaN: a NaN real part read as the finite imaginary magnitude, while a NaN imaginary part did return NaN. The imaginary side already lands on `im`, so the real side takes an explicit `!(re == re)` mask in the generic `Complex<T>` SIMD path, and the scalar fallbacks of `ComplexImpl<f32>` / `ComplexImpl<f64>` test both parts before the comparison. No in-tree type reaches either path (native `c32` / `c64` go through the real coercion in `norm_max`), so no dispatch, allocation or numerical result on an existing input changes. | KEEP | tensor4all/faer#320 | Correctness; the `faer-traits` unit test `complex_abs_max_propagates_nan` drives `SimdCtx::abs_max` on `Complex<f64>` and fails before / passes after (real-part NaN, imaginary-part NaN, both, and three finite controls). No performance claim: two mask ops per SIMD chunk on a path no in-tree type reaches. |
 | Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it normalizes the eigenvector by its norm and `(alpha, beta)` by `max(|alpha|, |beta|)`, which leaves the residual unchanged, and requires the two eigenvector norms, that scale and the residual denominator to be finite and positive; a numerator formed from the values as they come loses the residual instead, and no later division brings it back: `u = v = 1e-200 * e0` with `alpha = 1e-150` underflows it to `0`, and `alpha = f64::MAX / 2` with `beta = 0` sends the denominator to `inf`, whose quotient is `0` as well. It then asserts both residuals finite, with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero eigenvector column would otherwise read as a perfect residual. Measured: with `ul = Mat::<T>::zeros(n, n)` after the `gevd_cplx` call the test passed before these asserts and fails after with `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`, and the same mutation on `ur` fails with `column 0: right eigenvector has norm 0e0`; replacing each column's eigenvectors with `k * e0`, `k = 1.125 * (f64::MAX / den)`, passed with the single-division form (`nu * den = inf`, residual recorded as `0e0` against a staged `2.978e-1`) and fails with `c32 left residual 2.977833357222273e-1 >= 1.430511474609375e-5` after it; and `residual_oracle_rejects_unnormalized_eigenpairs` hands the oracle the review's two columns directly (`alpha = f64::MAX / 2, beta = 0, u = v = 1e-300 * e0` and `alpha = 1e-150, beta = 0, u = v = 1e-200 * e0`): on the oracle without normalization both read `left residual 0e0 < 1.430511474609375e-5` and `right residual 0e0 < 1.430511474609375e-5`, and with it `2.977833357222273e-1` and `2.941381743992797e-1`. The dispatch claim was measured too, with a temporary `eprintln!` at both `hessenberg_to_qz_unblocked` call sites and default `params`: a 20 x 20 pencil prints `n=20 nmin=75 blocking_threshold=75 -> unblocked`, and a 100 x 100 one prints `n=100 nmin=75 -> final unblocked call after the blocked loop`. The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). |
+| BLAS-level checks ported from the OpenBLAS `test/` and `ctest/` programs (#166) | Inventory of the check families in `test/` and `ctest/` at OpenBLAS `791b02a` against their faer counterparts, plus a first batch of bounded unit tests: triangular multiply (`DTRMM`), triangular solve (`DTRSM`/`DTRSV`), the rank-1 (`DGER`/`DSYR`) and symmetric rank-k (`DSYRK`) updates, and `norm_l2` marked covered. The inventory rows state which of the source's checks each port reproduces; it is not a claim that every check of every routine was compared. Each test asserts an independently accumulated reference product and the source's rogue-value properties (an unreferenced triangle and an implicit unit diagonal that must not be read, operands unchanged, padding untouched), with NaN in place of the source's `-1.0D10` and exact integers in place of its random fixture. In the `DTRMM` and `DTRSM` ports the written operand also lives in a padded parent, so the padding check covers the destination and that destination is checked entry by entry against the reference. Test-only: no production path changes. | KEEP | tensor4all/faer#166 | Correctness only, no timing claim. 8 tests over f64/f32/c64/c32, 32 in the default run: `cargo test -p t4a-faer --lib blas` (32 passed); two of them (`trsv`/`trsm` with an imaginary diagonal) return immediately on f32/f64, where `conj` is the identity; dropping the `maybe_conj_lhs` at `triangular_solve.rs:147` or `:156` turns the c64 cases red (got `-0.5i`, want `0.5i`); `linalg::matmul` 34 passed/1 ignored and `linalg::triangular_solve` 42 passed/8 ignored, measured at `cbb46be` (the 16 BLAS tests are inside that 42). Deliberate mutations fail the matching test and are reverted: unit-diagonal read in `diag_kind`, triangular destination widened to `Rectangular`, `rank_update` ignoring `beta`, the unit lower solve reading the stored diagonal, and — for the guarded destinations — the `DTRMM` product overwritten with NaN, the `DTRMM` operand's unreferenced triangle overwritten after the call, the `DTRSM` solution overwritten with NaN, one padding entry of either destination's parent overwritten, one padding entry of the `DTRSM` coefficient matrix's parent overwritten, the `(n, k) = (3, 1)` rank-1 product leaving its off-diagonal entries unwritten, the unconjugated `A A^T` arm handed a conjugation, its conjugated `A A^H` sibling handed none, the strided `norm_l1` fallback reading its input contiguously, and one entry of the rank-1 update's `x` overwritten after the call. Netlib BLAS credit: `faer/COPYING.BLAS.NETLIB`; inventory, source identities and omissions: [docs/blas-test-inventory.md](docs/blas-test-inventory.md). Not ported: banded/packed storage, `SYMM`/`SYMV`/`SYR2K`, 3m variants, mixed-precision comparisons, XERBLA argument validation, the `#[ignore]`d full GEMM sweep. |
 
 ## Ground rules