Metadata-Version: 2.4
Name: blaspy
Version: 0.0.6
Summary: Mini motor de álgebra lineal en C (BLA de BLAS, py de Python) con bindings a Python vía ctypes. El paquete de PyPI se llama 'blaspy', pero el módulo real es 'blapy' -- 'pip install blaspy' + 'import blapy'.
License: MIT
Classifier: Development Status :: 2 - Pre-Alpha
Classifier: Intended Audience :: Developers
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: C
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: MacOS
Classifier: Operating System :: Microsoft :: Windows
Requires-Python: >=3.7
Description-Content-Type: text/markdown
Dynamic: classifier
Dynamic: description
Dynamic: description-content-type
Dynamic: license
Dynamic: requires-python
Dynamic: summary

# BLApy

> ⚠️ **BETA — API inestable, cobertura parcial de BLAS, sin
> autotuning ni multithreading probado en hardware con más de un
> núcleo.** Este proyecto es utilizable y está validado con cientos
> de tests, pero no reemplaza a una librería BLAS madura para uso en
> producción. Ver "Qué falta" para el detalle completo de
> limitaciones conocidas antes de decidir si sirve para tu caso.

`BLApy` (BLA de BLAS, py de Python) es un mini motor de álgebra
lineal en C, con bindings a Python. **No busca reemplazar a BLAS**
(OpenBLAS, MKL, BLIS) — es un proyecto propio para tener un núcleo
matemático chico pero real, usando la misma técnica de fondo que usan
las implementaciones serias: **packing de memoria + microkernel con
acumulador en registros**, no solo un triple-loop con SIMD.

Uso principal: **C directo**. Los bindings de Python (vía `ctypes`,
sin dependencias de compilación) son para prototipar una idea rápido
sin escribir tanto C — una vez que algo funciona en Python, se
confirma y se usa desde C.

## Estado: 0.0.6

**Instalable vía `pip install blaspy`** (el paquete de PyPI se llama
`blaspy`, pero el módulo que se importa es `blapy` — ver "Instalar" 
más abajo) o clonando el repo y compilando con `make`.

Cubre los tres niveles clásicos de BLAS, en `double` y un subconjunto
creciente en `float32`:

- **Nivel 1** (vector-vector): `dot`, `axpy`, `scal`, `nrm2`, `asum`,
  `swap`, `copy`, `iamax`, `rot` (las últimas cuatro, **nuevas en
  0.0.6**) — las nueve, en `double` y `float32`.
- **Nivel 2** (matriz-vector): `gemv`, `ger` (producto exterior /
  rango 1 — añadido en 0.0.2), y **nuevas en 0.0.6**: `symv`
  (matriz simétrica), `trmv`/`trsv` (matriz triangular, producto y
  resolución de sistema) — las cinco, en `double` y `float32`.
- **Nivel 3** (matriz-matriz): `gemm` — la función central de este
  proyecto, ver la sección de diseño más abajo — más `gemm_tn` y
  `gemm_nt` (variantes con A o B transpuesta sin copiarla — añadidas
  en 0.0.2, `float32` en 0.0.6), y **nuevas en 0.0.6**: `symm`
  (matriz simétrica), `trmm`/`trsm` (matriz triangular, producto y
  resolución de sistema para todas las columnas/filas de una matriz a
  la vez), `syrk` (actualización simétrica de rango-k) — las siete,
  en `double` y `float32`. `symm`/`trmm`/`trsm` implementan ambos
  lados de multiplicación (`side`: A por izquierda o por derecha),
  fieles al estándar BLAS completo, no solo el caso más común.

Convención de nombres: `float32` usa el prefijo `s` (`sdot`, `sgemv`,
`ssymv`, ...) sobre el nombre de `double` — misma convención clásica
de BLAS (S=single, D=double).

Todas las funciones vectorizables (todo nivel 1/3, y `gemv`/`ger` de
nivel 2 — no `symv`/`trmv`/`trsv`, ver más abajo) detectan el nivel
SIMD disponible **en runtime** en x86_64 (vía `cpuid`), con una
cascada de fallback: AVX-512 → AVX2 → SSE2 → escalar puro. Esto
significa que un mismo `.so`, compilado una sola vez sin
`-march=native`, corre correcto en cualquier CPU x86_64 y usa lo
mejor que esa máquina en particular tenga disponible — no hace falta
recompilar por CPU de destino. En ARM64, todo lo vectorizable usa
NEON, incluyendo nivel 2 -- aunque `gemv`/`ger` no tienen un
microkernel NEON *dedicado*, SÍ corren sobre NEON en la práctica
porque llaman a `dot`/`axpy` (nivel 1), que sí lo tienen (ver
`bla_level2.c` para el detalle, incluida la investigación de si un
microkernel dedicado ayudaría — sin forma de medirlo en el entorno de
desarrollo de este proyecto). `symv`/`trmv`/`trsv` son escalares
puras en todas las arquitecturas — sin camino SIMD dedicado en
ninguna versión todavía, ver "Qué falta". `bla_simd_info()` devuelve
el nombre del nivel detectado (`avx512`, `avx2`, `sse2`, `neon`, o
`scalar`), para poder confirmarlo sin adivinar.

### Novedades de 0.0.6

- **Nombre de distribución en PyPI cambiado a `blaspy`** (el módulo
  importable sigue siendo `blapy`, sin cambios). `setup.py` separa
  `name=` (lo que ve el índice de paquetes) de `packages=`/
  `package_dir=` (la carpeta real que queda instalada e
  importable) — confirmado instalando de punta a punta: el wheel se
  llama `blaspy-0.0.6-...`, `pip show blaspy` reporta la versión
  correcta, `import blaspy` falla como corresponde (no existe ningún
  módulo con ese nombre), e `import blapy` funciona completo.
- **Ejemplos del README corregidos — no compilaban/corrían tal como
  estaban escritos.** El bloque de C no tenía `main()` (a nivel de
  archivo, C no permite ejecutar sentencias sueltas — el compilador
  interpretaba cada llamada como el comienzo de una declaración de
  función, con errores como "initializer element is not a
  compile-time constant"), y el comando de compilación sugerido le
  faltaba `-Wl,-rpath,.` (compilaba y enlazaba sin error, pero fallaba
  al *correr* el binario con "cannot open shared object file" — el
  link exitoso no confirma que el binario sepa dónde buscar el `.so`
  en tiempo de ejecución). El bloque de Python corría sin error, pero
  ninguna línea imprimía su resultado, así que copiarlo y pegarlo no
  mostraba los valores `# 32.0` que los comentarios prometían. Los
  tres bugs se encontraron recién al copiar el bloque exacto del
  README y ejecutarlo tal cual — no alcanzaba con leerlo o con
  probar variantes propias escritas a mano, que sí tenían el `main()`
  puesto correctamente desde el principio.
- **Nivel 1 completo en ambas precisiones**: `swap`, `copy`, `iamax`,
  `rot` — las cuatro funciones de BLAS Level 1 que quedaban fuera.
  `copy` usa `memcpy` de la librería estándar en vez de un camino
  SIMD escrito a mano (una implementación de `memcpy` madura,
  típicamente con dispatch por CPU real vía `ifunc` en glibc,
  difícilmente se supera escribiendo algo a mano para este caso —
  mismo criterio que ya llevó a implementar `nrm2` como
  `sqrt(dot(x,x))` en vez de un camino nuevo). `iamax` es 0-indexada
  (no 1-indexada como el BLAS de referencia en Fortran) — decisión
  deliberada por coherencia con el resto de las firmas de este
  proyecto, documentada explícitamente en el código para quien
  conozca BLAS estándar y no lo espere. Ninguna de las cuatro tiene
  camino SIMD dedicado para `iamax`/`rot` en esta versión (encontrar
  un índice, o necesitar dos valores por posición simultáneamente,
  es más complejo de vectorizar que las reducciones de suma que ya
  existían — queda como candidato para una versión futura si se mide
  que vale la pena).
- **Nivel 2 simétrico/triangular, en ambas precisiones**: `symv`
  (`y = alpha*A@x + beta*y`, A simétrica), `trmv` (`x = A@x` en el
  lugar, A triangular), `trsv` (resuelve `A@x = b` para `x` en el
  lugar, A triangular). Las tres usan un parámetro `uplo` (0=lower,
  1=upper) para indicar qué triángulo de A está poblado con valores
  reales — el otro triángulo puede contener basura sin inicializar
  (en `symv`, porque nunca se lee — se usa el elemento espejado por
  simetría; en `trmv`/`trsv`, porque se trata matemáticamente como
  cero). Validado con basura *deliberada* en el triángulo no usado en
  cada test, para confirmar que nunca se lee esa memoria, y `trsv`
  además con verificación round-trip (`A @ x_resultado` reconstruye
  el `b` original) en vez de solo comparar contra un valor
  precalculado. **Sin camino SIMD dedicado en ninguna arquitectura**
  — el patrón de acceso (mezcla de lectura por fila y por columna
  según qué mitad de la matriz está poblada, o dependencia
  secuencial estricta en el caso de `trsv`) no se adapta directo al
  mismo enfoque de bloque fijo que usa el resto del proyecto; queda
  como trabajo futuro si se mide que vale la pena.
- **Nivel 3 simétrico/triangular, en ambas precisiones**: `symm`
  (equivalente de nivel 3 a `symv`), `trmm`/`trsm` (equivalentes de
  nivel 3 a `trmv`/`trsv` — resuelven/multiplican para TODAS las
  columnas de `B` a la vez, no un vector), `syrk` (`C = alpha*A@A^T +
  beta*C`, estructuralmente distinta a las otras tres — una sola
  matriz de entrada, no dos). A diferencia de `symv`/`trmv`/`trsv`
  (que solo cubrían el caso implícito de "A multiplica por la
  izquierda"), `symm`/`trmm`/`trsm` implementan un parámetro `side`
  explícito (0=left, 1=right) — ambos lados de BLAS real, no solo el
  más común, decisión de alcance completo confirmada antes de
  escribir código. Cada fórmula (10 combinaciones side×uplo×trans en
  total, contando las cuatro funciones) se calculó y verificó a mano
  con NumPy antes de escribir una sola línea de C — en particular,
  `trsm` con `side=right` no se derivó "por simetría del problema"
  sin más: se confirmó explícitamente que la sustitución recorre las
  columnas de A en sentido opuesto al de `side=left` para el mismo
  `uplo`, un detalle que una derivación apurada podría haber
  invertido sin que ningún test lo detectara si el caso de prueba
  hubiera sido demasiado simétrico. `syrk` solo lee/escribe el
  triángulo de `C` indicado por `uplo` — el otro triángulo no se toca
  en absoluto, ni para leer el `beta*C` existente ni para escribir el
  resultado. Mismas limitaciones que nivel 2 simétrico/triangular:
  sin SIMD dedicado, `A` sigue siendo una matriz completa `n×n` (no
  el formato "packed" real de BLAS). Validado con 160 combinaciones
  por precisión (todas las combinaciones de `side`/`uplo`/`trans`,
  tamaños "feos" no cuadrados, basura deliberada en las regiones no
  usadas, y verificación round-trip para `trsm`) — pero, igual que
  `gemm`/`gemm_tn`/`gemm_nt`/`sgemm`/`sgemm_tn`/`sgemm_nt`, estas
  cuatro funciones viven en `bla_gemm.c`/`bla_gemm_f32.c`, que ahora
  requieren OpenMP obligatorio — así que **no pudieron
  cross-compilarse con el toolchain bare-metal usado para NEON/ARM64
  en este proyecto** (mismo motivo exacto que gemm, ver más abajo),
  solo se validaron en las cuatro rutas SIMD x86. Ver
  `tests/VALIDAR_EN_ANDROID.md`, actualizado con estos dos tests
  nuevos para quien quiera confirmarlos en hardware ARM64 real.
  `sgemm_tn`/`sgemm_nt`, mismo patrón que sus equivalentes `double`
  (añadidos en 0.0.2). Se extrajo `gemm_f32_impl` de lo que antes era
  el cuerpo directo de `bla_sgemm` — mismo refactor que ya se había
  hecho para `double` en 0.0.2, ahora también en `float32`.
- **Microkernel AVX-512 de `double` ensanchado de 8×8 a 8×16** —
  hasta 0.0.5, ese camino usaba solo 8 de los 32 registros ZMM
  disponibles para acumuladores. Medido antes de integrarlo (mismo
  criterio de "no aplicar cambios sin medir" del resto del proyecto):
  un primer intento con pocas repeticiones (2000) dio una varianza
  enorme entre corridas (0.22x a 4.48x) — ruido del entorno de
  desarrollo dominando la medición, no señal real. Con 10x más
  repeticiones (20000) para promediar mejor ese ruido: 7 de 7
  corridas con 8×16 más rápido, en un rango angosto (1.08x-1.23x).
  También se probó 16×16 (el banco completo de 32 registros) contra
  8×16 — resultado: NO ayuda (0.92x-1.00x en 5 de 5 corridas,
  consistentemente igual o levemente peor), probable presión de
  registros al no dejar margen para los operandos de la operación en
  sí. AVX2 y escalar NO se tocaron — ya estaban en su óptimo medido
  desde 0.0.4 (banco de registros YMM lleno, con evidencia de
  spilling ya documentada entonces). Esto requirió que `mr`/`nr`
  pasaran a ser parámetros en runtime en el microkernel de `double`
  (mismo patrón que `float32` ya tenía desde 0.0.3) — antes,
  `BLA_MR`/`BLA_NR`=8 servía igual para las tres rutas x86; desde
  esta versión, ya no. Ganancia medida de punta a punta en `bla_gemm`
  completo: mayor que la del microkernel aislado (ver benchmark
  actualizado más abajo), consistente con que el microkernel domina
  una fracción todavía mayor del tiempo total en `gemm` real que en
  la medición aislada.
- **Multithreading con OpenMP**, en `gemm`/`gemm_tn`/`gemm_nt` y sus
  equivalentes `float32` — la limitación que 0.0.5 documentaba como
  "no investigada por falta de un entorno con más de un núcleo". Se
  paraleliza el loop `ic` (bloques de `BLA_MC` filas de `A`/`C`): es
  el que tiene más iteraciones reales en los tamaños donde este
  proyecto se probó (a diferencia de `jc`, que con `BLA_NC=2048`
  suele tener una sola iteración), y cada iteración escribe a un
  rango de filas de `C` que nunca se solapa con el de otra — no hace
  falta ningún lock. El buffer de packing de `A` (`Ap`) pasa a ser
  privado por thread (se escribe en cada iteración de `ic`); el de
  `B` (`Bp`) sigue siendo un único buffer compartido (es de solo
  lectura dentro de ese loop). **OpenMP es obligatorio desde esta
  versión** — no hay build sin `-fopenmp`, una decisión consciente de
  simplicidad sobre compatibilidad con entornos sin `libgomp`.
  Validado con dos niveles de rigor: la suite de correctitud completa
  forzando distintas cantidades de threads (2, 4, 8, 16 — más threads
  que núcleos físicos reales en el entorno de desarrollo, que solo
  tiene 1), y un test de estrés dedicado que corre el mismo cálculo
  20 veces bajo 9 configuraciones de threads distintas, confirmando
  resultado **bit-idéntico** entre corridas (una condición de carrera
  real produciría resultados *distintos* entre corridas de la misma
  entrada, no solo un resultado incorrecto una vez — comparar
  corrida-contra-corrida es más sensible a esto que comparar contra
  una referencia fija una sola vez). **Confirmado además en hardware
  ARM64 real** (un teléfono Android vía Termux, no el sandbox de
  desarrollo) — 80 corridas repetidas (4 configuraciones de threads ×
  20 corridas), 0 divergencias, la primera confirmación de esto bajo
  paralelismo físico genuino (núcleos reales compitiendo por memoria,
  no el time-slicing de un único core que es todo lo que el sandbox
  de desarrollo puede ofrecer) — ver `tests/VALIDAR_EN_ANDROID.md`.
  **Pendiente**: medir si esto da *speedup* real en algún hardware con
  varios núcleos — el sandbox de desarrollo (1 CPU) solo puede
  confirmar correctitud bajo time-slicing, nunca paralelismo
  genuino, y la corrida de Android que sí confirmó correctitud real
  todavía no incluyó esa medición de velocidad.

### Comparación contra NumPy

Medido en este entorno de desarrollo (NumPy con backend OpenBLAS
0.3.31, `MAX_THREADS=64`, confirmado con `np.show_config()` —
aunque el sandbox de desarrollo solo expone 1 CPU física, así que
OpenBLAS tampoco puede usar threads reales acá, mismo límite que
tiene la paralelización OpenMP de BLApy en este entorno; ver
`bench/compare_numpy.py` para reproducir):

| n    | BLApy `gemm` (`double`) | NumPy (`double`) | ratio | BLApy `sgemm` (`float32`) | NumPy (`float32`) | ratio |
|------|--------------------------|-------------------|-------|------------------------------|----------------------|-------|
| 32   | 0.006 ms                 | 0.002 ms          | 0.37x | 0.005 ms                     | 0.002 ms             | 0.35x |
| 64   | 0.018 ms                 | 0.008 ms          | 0.46x | 0.012 ms                     | 0.006 ms             | 0.50x |
| 128  | 0.091 ms                 | 0.077 ms          | 0.85x | 0.052 ms                     | 0.040 ms             | 0.78x |
| 256  | 0.622 ms                 | 0.555 ms          | 0.89x | 0.324 ms                     | 0.248 ms             | 0.76x |
| 384  | 2.00 ms                  | 1.84 ms           | 0.92x | 0.90 ms                      | 0.84 ms              | 0.93x |
| 512  | 4.77 ms                  | 4.32 ms           | 0.91x | 2.78 ms                      | 1.91 ms              | 0.69x |
| 768  | 16.17 ms                 | 13.57 ms          | 0.84x | 6.88 ms                      | 6.25 ms              | 0.91x |
| 1024 | 38.64 ms                 | 32.39 ms          | 0.84x | 18.05 ms                     | 15.52 ms             | 0.86x |

**Tendencia real**: NumPy/OpenBLAS más rápido en absolutamente todo
el rango medido — nunca hay un tamaño donde `bla_gemm`/`bla_sgemm`
ganen en esta máquina. La brecha es más grande en tamaños chicos
(0.35x-0.50x en `n=32`/`n=64` — el trabajo total es tan poco que el
costo fijo de armar los buffers de packing pesa proporcionalmente
más) y se achica notablemente hacia el medio del rango (0.84x-0.93x
en `n=256` a `n=1024` — la brecha *no* sigue una tendencia monótona
simple con `n`, a diferencia de lo que sugerían mediciones de
versiones anteriores de este proyecto). Quien reproduzca esta tabla
en su propia máquina probablemente obtenga números absolutos y
ratios distintos — la comparación depende fuerte de qué CPU concreta
se use (la disponibilidad de AVX-512, el tamaño real de cada nivel
de cache) y de si hay más de un núcleo físico disponible (donde
OpenBLAS tiene una ventaja estructural que BLApy no puede igualar en
este entorno de desarrollo).

**Sobre el microkernel 8×16 (0.0.6):** antes de ese cambio, `n=512`
medía 0.72x; después, 0.91x — una mejora real de punta a punta, más
grande que el ~10-23% medido en el microkernel aislado (ver
"Novedades de 0.0.6"), consistente con que el microkernel domina una
fracción todavía mayor del tiempo total en `gemm` completo que en la
medición aislada. Sigue sin cerrar la brecha — el camino de mayor
impacto potencial que queda es medir si el multithreading (ya
implementado, ver "Novedades de 0.0.6") da speedup real en un
hardware con varios núcleos físicos, algo que este entorno de
desarrollo no puede confirmar.


## Por qué `bla_gemm` es más rápido que un triple-loop con SIMD

Un `matmul` con buen tiling de cache y SIMD por fila (reordenar loops
a i-k-j, bloquear en `(kb, jb)`, vectorizar el loop interno) ya
resuelve el problema de *localidad* de memoria: evita relecturas
innecesarias de cache. Pero el acumulador de salida sigue viviendo en
memoria/cache durante todo el barrido de `k` — cada paso de `k` es
una lectura y escritura real a `C[i,j]`, aunque esa lectura/escritura
casi siempre sea un cache-hit rápido.

`bla_gemm` va un paso más allá con dos técnicas combinadas (la
técnica de Goto, popularizada por
[BLIS](https://github.com/flame/blis) — ver cita más abajo):

1. **Packing**: antes de multiplicar, se copian bloques de `A` y `B`
   a buffers temporales, reescritos en el orden EXACTO en que el
   microkernel los va a leer. El microkernel nunca sufre un stride
   raro ni un salto de fila — siempre lee memoria contigua.
2. **Microkernel con acumulador en registros**: un bloque fijo del
   resultado (8×16 en AVX-512 desde 0.0.6, 8×8 en AVX2/escalar/NEON —
   ver "Novedades de 0.0.6" para por qué AVX-512 es distinto) se
   mantiene **en registros SIMD del CPU**, no en cache, durante TODO
   el barrido de `k`. Cada paso de `k` es una carga de `A` empaquetada
   + una carga de `B` empaquetada + FMA contra el acumulador en
   registro — cero tráfico de memoria para el acumulador hasta el
   final del barrido.

> BLIS distilled the de-facto Goto algorithm for high performance
> matrix-matrix multiplication into a single architecture-specific
> computation microkernel and two packing routines. By providing
> custom implementations of the microkernel and/or packing routines,
> [...] an entire BLAS library can be generated for a specific
> architecture.
> — [SMaLL: A Software Framework for portable Machine Learning
> Libraries](https://arxiv.org/pdf/2303.04769), citando el diseño de
> BLIS

`bla_gemm` usa exactamente esa idea: un microkernel vectorizado a
mano en AVX-512 (8×16 desde 0.0.6, antes 8×8), AVX2 (8×8, 0.0.2), y
NEON/ARM64 (8×8, 0.0.2) — con fallback escalar correcto para SSE2 y
cualquier otro caso —, y packing con padding a cero en los bordes
para que el microkernel nunca necesite un caso especial para bloques
parciales.

### Benchmark: `bla_gemm` vs `nx_matmul` (de la librería `numerx`)

Medido en este entorno de desarrollo, matrices cuadradas aleatorias,
mediana de 7 repeticiones, gemm/matmul puro (sin incluir conversión
de tipos ni construcción de arrays en el tiempo medido), ambas
librerías compiladas con `-O3` **sin** `-march=native` (para que la
comparación sea justa: mismas condiciones de portabilidad en ambos
lados, no una ventaja artificial de una sobre otra). Esta CPU tiene
AVX-512, así que BLApy usa ese camino por defecto; `numerx` se queda
en AVX2 (es su techo, ver su propio README):

| n   | BLApy `gemm` AVX-512 (ms) | numerx `matmul` AVX2 (ms) | speedup |
|-----|----------------------------|-----------------------------|---------|
| 32  | 0.006                      | 0.009                       | 1.6x    |
| 64  | 0.019                      | 0.057                       | 2.9x    |
| 128 | 0.110                      | 0.513                       | 4.7x    |
| 256 | 0.686                      | 4.192                       | 6.1x    |
| 384 | 2.511                      | 13.538                      | 5.4x    |
| 512 | 6.184                      | 36.479                      | 5.9x    |

`numerx.matmul` ya usa tiling de cache + SIMD por fila (ver su propio
README) — no es una comparación contra un triple-loop ingenuo. La
ventaja de `bla_gemm` viene específicamente del packing + microkernel
con acumulador en registros, no de tener SIMD que el otro no tenga —
aunque tener AVX-512 disponible en este CPU puntual también ayuda,
ver el punto siguiente.

**El resultado más importante de 0.0.2, comparado con 0.0.1:** en
0.0.1, si se forzaba el camino AVX2 de BLApy (simulando una máquina
sin AVX-512 — ver `make bench-forced`), `bla_gemm` caía al fallback
escalar y perdía contra `numerx.matmul` (~1.8x más lento en
`n=256`). Con el microkernel AVX2 vectorizado a mano de 0.0.2, el
mismo experimento en `n=256` da `bla_gemm`(AVX2 forzado) ≈ 1.20 ms
contra `numerx.matmul`(AVX2) ≈ 4.19 ms — **BLApy ahora es ~3.5x más
rápido que numerx incluso sin AVX-512 disponible**, revirtiendo por
completo la pérdida que tenía 0.0.1. Reproducir con
`make bench-forced` (que fuerza el runtime-dispatch a cada nivel sin
necesitar hardware distinto) y comparar contra `python3
bench/compare_numerx.py` corriendo en la misma máquina.

En `n=32`, la ventaja es menor porque el trabajo total es tan chico
que el costo fijo de armar los buffers de packing pesa
proporcionalmente más — para matrices muy chicas, el overhead de
packing puede no compensarse (ver sección "Qué falta" más abajo, caso
concreto de matrices diminutas en un loop). Reproducir la tabla
completa: `make bench` (requiere tener numerx compilado como
extensión de Python en el mismo entorno; ver
`bench/compare_numerx.py`).

## Uso desde C

```c
#include "bla.h"
#include <stdio.h>

int main(void) {
    // Nivel 1: producto punto
    double x[] = {1.0, 2.0, 3.0};
    double y[] = {4.0, 5.0, 6.0};
    double result = bla_dot(x, y, 3);  // 32.0

    // Nivel 1, nuevo en 0.0.6: swap, intercambia en el lugar
    double xa[] = {1.0, 2.0};
    double ya[] = {9.0, 8.0};
    bla_swap(xa, ya, 2);  // xa = {9,8}, ya = {1,2}

    // Nivel 3: C = A @ B (A es 2x3, B es 3x2, C es 2x2)
    double A[] = {1, 2, 3, 4, 5, 6};
    double B[] = {7, 8, 9, 10, 11, 12};
    double C[4];  // el LLAMADOR reserva C -- bla_gemm no reserva memoria
    bla_gemm(1.0, A, 2, 3, B, 2, 0.0, C);

    // Nivel 2: A = x @ y^T + A (producto exterior, in-place sobre A)
    double x2[] = {1.0, 2.0};
    double y2[] = {3.0, 4.0, 5.0};
    double A2[6] = {0};  // A es 2x3, mutada en el lugar
    bla_ger(1.0, x2, 2, y2, 3, A2);  // A2 = {3,4,5, 6,8,10}

    // Nivel 2, nuevo en 0.0.6: symv, con A simétrica pasada como
    // "solo el triángulo inferior" (uplo=0) -- el resto de A puede
    // ser cualquier cosa, nunca se lee.
    double As[] = {
        2, 0, 0,   // triángulo superior: no se lee, puede ser basura
        1, 5, 0,
        3, 4, 6
    };
    double xs[] = {1.0, 2.0, 3.0};
    double ys[] = {10.0, 20.0, 30.0};
    bla_symv(1.5, As, 3, 0 /* lower */, xs, 0.5, ys);  // ys = {24.5, 44.5, 58.5}

    // Nivel 3: C = A_stored^T @ B, sin materializar la transpuesta.
    // A_stored está almacenada 3x2 (k x m) -- representa A^T.
    double A_stored[] = {1, 2, 3, 4, 5, 6};  // A logica (2x3) es {{1,3,5},{2,4,6}}
    double B2[] = {1, 0, 0, 1, 1, 1};        // 3x2
    double C2[4];
    bla_gemm_tn(1.0, A_stored, 2, 3, B2, 2, 0.0, C2);  // C2 = {6,8, 8,10}

    // float32: mismas firmas de forma, con prefijo "s" y tipo float en
    // vez de double. Mismos valores que el ejemplo de bla_gemm arriba,
    // para poder comparar el resultado directo.
    float xf[] = {1.0f, 2.0f, 3.0f};
    float yf[] = {4.0f, 5.0f, 6.0f};
    float resultf = bla_sdot(xf, yf, 3);  // 32.0f

    float Af[] = {1, 2, 3, 4, 5, 6};
    float Bf[] = {7, 8, 9, 10, 11, 12};
    float Cf[4];
    bla_sgemm(1.0f, Af, 2, 3, Bf, 2, 0.0f, Cf);  // Cf = {58,64, 139,154}

    printf("dot=%.1f swap=%.1f,%.1f gemm=%.1f,%.1f,%.1f,%.1f symv=%.1f,%.1f,%.1f sdot=%.1f\n",
           result, xa[0], xa[1], C[0], C[1], C[2], C[3], ys[0], ys[1], ys[2], resultf);
    return 0;
}
```

Compilar y enlazar (`-I.` es necesario para que `#include "bla.h"`
encuentre el header si tu_programa.c no está en el mismo directorio
que bla.h -- si estás compilando desde DENTRO del repo clonado y
bla.h está ahí mismo, `-I.` es redundante pero no hace daño):

```sh
gcc tu_programa.c -I. -L. -lbla -lm -Wl,-rpath,. -o tu_programa
```

(`-L.` asume que `libbla.so` está en el directorio actual -- ajustá
la ruta si no es así. `-Wl,-rpath,.` es IMPRESCINDIBLE, no opcional:
sin esto, el binario compila y enlaza sin error, pero falla al
CORRERLO con `error while loading shared libraries: libbla.so:
cannot open shared object file` -- el link exitoso solo confirma que
el compilador encontró los símbolos en tiempo de compilación, no que
el binario sepa dónde buscar el .so en tiempo de ejecución; sin
rpath, el sistema solo busca en las rutas estándar del sistema, no en
el directorio actual. Alternativa sin rpath:
`LD_LIBRARY_PATH=. ./tu_programa` al correrlo, en vez de en la
compilación. Con el paquete instalado vía pip, ver "Compilar e
instalar" más abajo para las rutas exactas de bla.h y libbla.so
dentro del paquete.)

## Uso desde Python

```python
import blapy

print(blapy.simd_info())  # 'avx512', 'avx2', 'sse2', 'neon', o 'scalar'

print(blapy.dot([1, 2, 3], [4, 5, 6]))  # 32.0

x, y = blapy.swap([1.0, 2.0], [9.0, 8.0])  # nuevo en 0.0.6
print(x, y)  # [9.0, 8.0] [1.0, 2.0]

# gemm acepta listas planas en row-major, o matrices 2D de NumPy
# (se aplanan automáticamente si tenés numpy instalado)
C = blapy.gemm(A=[1,2,3,4,5,6], m=2, k=3, B=[7,8,9,10,11,12], n=2)
print(C)  # [58.0, 64.0, 139.0, 154.0]

# symv, nuevo en 0.0.6: A simétrica, solo el triángulo inferior
# poblado (uplo=0) -- el resto puede ser cualquier cosa
A_lower = [2,0,0, 1,5,0, 3,4,6]
y_symv = blapy.symv(1.5, A_lower, 3, 0, [1.0, 2.0, 3.0], 0.5, [10.0, 20.0, 30.0])
print(y_symv)  # [24.5, 44.5, 58.5]

# float32: mismo patrón, prefijo "s"
print(blapy.sdot([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]))  # 32.0
Cf = blapy.sgemm(A=[1,2,3,4,5,6], m=2, k=3, B=[7,8,9,10,11,12], n=2)
print(Cf)  # [58.0, 64.0, 139.0, 154.0]
```

Todas las funciones de C tienen su binding de Python equivalente
(mismo nombre sin el prefijo `bla_`) — incluyendo las nuevas de
0.0.6 (`swap`, `copy`, `iamax`, `rot`, `symv`, `trmv`, `trsv`, y sus
variantes `float32`). `ger`, `gemm_tn`, y `gemm_nt` también están
disponibles (existían en C desde 0.0.2, pero no tenían binding de
Python hasta 0.0.3 — omisión corregida entonces).

Los bindings de Python cargan `libbla.so` vía `ctypes` — no hace
falta compilar una extensión de Python aparte, solo tener el `.so` de
BLApy compilado y en el mismo directorio que `blapy.py` (o en el
directorio de trabajo actual). Ver `python/blapy.py` para el detalle
de cómo se busca la librería.

## Instalar (vía pip) o compilar (desarrollo local)

**Instalar como paquete** (nuevo en 0.0.4; nombre de distribución
`blaspy` desde 0.0.6 — ver más abajo):

```bash
pip install blaspy
```

El paquete de PyPI se llama `blaspy`, pero el módulo que se importa
sigue siendo `blapy` (`import blapy`, no `import blaspy`) — mismo
proyecto, dos nombres distintos: uno para el índice de paquetes, otro
para el código. `setup.py` separa `name=` (lo que ve PyPI) de
`packages=`/`package_dir=` (la carpeta real instalada). Instalando
desde el repo clonado en vez de PyPI, `pip install .` funciona igual
(usa el `name=` de `setup.py`, que ya es `blaspy`).

Esto compila `libbla.so` y lo coloca dentro del paquete `blapy`
instalado — después, `import blapy` funciona desde cualquier
directorio, sin necesitar clonar el repo ni tener `libbla.so` a mano.
Ver "Novedades de 0.0.4" más arriba para el detalle de por qué esto
no fue trivial (BLApy no es una extensión de Python real, es una
librería C plana cargada vía `ctypes`, lo cual `setuptools` no maneja
bien por defecto).

**Requiere OpenMP (`libgomp`) desde 0.0.6** — `bla_gemm`/`bla_sgemm`
y sus variantes transpuestas paralelizan con `#pragma omp`, y esto es
obligatorio, no opcional (ver "Novedades de 0.0.6"). En la gran
mayoría de los sistemas Linux/macOS con `gcc`/`clang` esto ya está
disponible sin instalar nada aparte; si el build falla por no
encontrar `omp.h` o `libgomp`, instalar el paquete de desarrollo de
OpenMP de tu distribución (por ejemplo `libomp-dev` en Debian/Ubuntu)
antes de reintentar.

**Desarrollo local** (compilar sin instalar, para trabajar sobre el
propio código fuente):

```bash
make          # compila libbla.so + corre todos los tests
```

El código fuente de los bindings de Python vive en un solo lugar,
`src/blapy/__init__.py` — es lo que `setup.py` empaqueta. `python/`
(el directorio que usa `make lib` para dejar `libbla.so` accesible
durante desarrollo local, sin necesitar `pip install`) contiene
`blapy.py` como un **symlink** hacia `src/blapy/__init__.py`, no una
copia — antes de 0.0.4 eran dos archivos separados que casi se
desincronizaron sin que nadie lo notara (se encontraron idénticos por
casualidad al revisar antes de empaquetar 0.0.4, no por diseño), así
que se unificaron en una sola fuente de verdad. `bla.h` sigue el mismo
patrón desde 0.0.5: `src/blapy/bla.h` es un symlink hacia el `bla.h`
real de la raíz del repo, para que `package_data` pueda incluirlo en
el paquete instalado sin duplicar la fuente de verdad (confirmado
compilando C contra el paquete `pip` instalado, desde un directorio
sin ningún archivo del repo — funciona).

## Compilar

```bash
make          # compila libbla.so + corre todos los tests
make lib      # solo compila libbla.so
make test     # solo corre los tests (requiere lib ya compilada... make se encarga)
make bench    # compara contra numerx.matmul (requiere numerx instalado)
make bench-forced  # mide bla_gemm (double) forzando cada nivel SIMD x86
                    # (avx512/avx2/sse2/none), sin necesitar hardware
                    # distinto -- ver bench/bla_simd_forced.c. Mide
                    # bla_gemm, NO bla_sgemm (float32) -- no hay
                    # bench-forced-f32 todavía, ver "Qué falta".
make test-arm64 ARM_TOOLCHAIN=<ruta> ARM_QEMU=<ruta>
              # corre la suite de tests en un binario ARM64 real, vía
              # cross-compiler + QEMU (ver tests/run_arm64_tests.sh y
              # tests/VALIDAR_EN_ANDROID.md para medir velocidad real).
              # LIMITACIÓN DESDE 0.0.6: el toolchain bare-metal usado
              # para esto NO soporta OpenMP (falta -pthread, normal en
              # un target sin sistema operativo detrás) -- los tests
              # que dependen de gemm/gemm_tn/gemm_nt/sgemm/sgemm_tn/
              # sgemm_nt (que ahora tienen #pragma omp obligatorio) NO
              # pueden cross-compilarse con este toolchain, solo los
              # que no tocan esos archivos. Ver
              # tests/VALIDAR_EN_ANDROID.md para la forma real de
              # validar esos tests en ARM64 desde 0.0.6 (clang nativo
              # en Termux/Android, que sí corre sobre un kernel Linux
              # real con pthreads completos).
make native   # compila con -march=native -- MÁS RÁPIDO en esta CPU puntual,
              # pero el .so resultante NO es portable a otras máquinas.
              # No usar para distribuir, solo para benchmarking local.
```

Por defecto (`make`, `make lib`), se compila con `-O3` pero **sin**
`-march=native` — mismo criterio que documenta el propio `numerx` en
su `setup.py`: la detección de SIMD ya es en runtime, así que
`-march=native` solo arriesgaría generar instrucciones incompatibles
con la CPU de quien instale el paquete después, sin ninguna ganancia
que ese runtime-dispatch no dé ya. Las funciones que usan intrínsecos
AVX2/AVX-512 llevan `__attribute__((target(...)))` por función (mismo
patrón que usa `numerx.c` para su propio código AVX2) para que
compilen correctamente sin necesitar `-march=native` a nivel de
archivo completo.

## Parámetros de bloqueo

Definidos en `bla_internal.h` (no son API pública, pueden cambiar).

### `bla_gemm` (`double`)

`BLA_MR`/`BLA_NR` fijos (8×8), definidos en `bla_internal.h`, ya NO
son el tamaño real del bloque en las tres rutas x86 desde 0.0.6 —
solo AVX2/escalar/NEON usan ese tamaño; AVX-512 usa 8×16 (ver
"Novedades de 0.0.6"), calculado en runtime por `gemm_block_size`
(`bla_gemm.c`), mismo mecanismo que `bla_sgemm` ya tenía desde 0.0.3
(ver más abajo):

| Parámetro | Valor | Qué es |
|-----------|-------|--------|
| `BLA_MR`  | 8     | Filas de A en el microkernel, para AVX2/escalar/NEON (viven en registros) |
| `BLA_NR`  | 8     | Columnas de B, para AVX2/escalar/NEON (ídem) — AVX-512 usa `nr=16` en su lugar, no esta macro |
| `BLA_KC`  | 256   | Bloque de la dimensión `k` para el buffer de packing |
| `BLA_MC`  | 96    | Bloque de la dimensión `m` para el buffer de packing |
| `BLA_NC`  | 2048  | Bloque de la dimensión `n` para el buffer de packing |

| Arquitectura   | `mr` | `nr` | Acumuladores | Registros que satura |
|----------------|------|------|---------------|------------------------|
| AVX-512 (0.0.6)| 8    | 16   | 16            | 16 de 32 ZMM (`__m512d`, 8 doubles c/u) |
| AVX2           | 8    | 8    | 16            | 16 YMM (`__m256d`, 4 doubles c/u) — óptimo medido, ver "Novedades de 0.0.6" |
| NEON           | 8    | 8    | 32            | 32 v0-v31 (`float64x2_t`, 2 doubles c/u) |
| SSE2/escalar   | 8    | 8    | —             | Sin microkernel dedicado, camino escalar |

`BLA_KC`/`BLA_MC`/`BLA_NC` son valores fijos razonables para CPUs
x86_64 modernos típicos — **no** se autotunearon para ninguna máquina
puntual (ver "Qué falta"). El 8×16 de AVX-512 usa 16 de los 32
registros ZMM disponibles, no el banco completo — se midió que 16×16
(el banco completo) no ayuda, ver "Novedades de 0.0.6" para el
detalle.

### `bla_sgemm` (`float32`, añadido en 0.0.3)

A diferencia de `bla_gemm`, acá el bloque lógico **cambia según la
arquitectura detectada en runtime**, no es el mismo MRxNR para las
tres rutas x86 — porque un registro `float32` de un ancho en bits
dado carga el DOBLE de elementos que uno `double` del mismo ancho, así
que el bloque óptimo también es más ancho, y ese "más ancho" no
escala igual en las tres arquitecturas (ver `bla_internal.h` para el
razonamiento completo):

| Arquitectura | `MR` | `NR` | Acumuladores | Registros que satura |
|--------------|------|------|---------------|------------------------|
| AVX-512      | 16   | 32   | 32            | 32 ZMM (`__m512`, 16 floats c/u) |
| AVX2         | 8    | 16   | 16            | 16 YMM (`__m256`, 8 floats c/u) |
| NEON         | 8    | 16   | 32            | 32 v0-v31 (`float32x4_t`, 4 floats c/u) |

Estos valores salieron de **medir varias combinaciones en el
contexto de un gemm completo** antes de fijarlos (no de una fórmula
ni de una suposición) — ver "Novedades de 0.0.3" más arriba para el
resultado de esa investigación en AVX-512 y AVX2. **NEON es la
excepción**: no se pudo medir tiempo de ejecución con el entorno de
desarrollo usado en este proyecto (ver "Qué falta"), así que su
`8×16` es una extrapolación razonada del mismo patrón que sí se
confirmó en las otras dos arquitecturas ("saturar el banco de
registros rinde bien"), no una medición directa — pendiente de
confirmar o corregir con hardware real.

`BLA_KC_F32`/`BLA_MC_F32`/`BLA_NC_F32` siguen el mismo razonamiento de
bloqueo de cache que sus equivalentes de `double`, ajustados
proporcionalmente (un `float` pesa la mitad que un `double`, así que
`KC` puede ser el doble para el mismo tamaño en bytes del buffer de
packing) — tampoco autotuneados.

## Qué falta (honesto, no escondido)

Este proyecto es **beta**: utilizable y validado con cientos de
tests, pero con cobertura parcial de BLAS y varias limitaciones
reales que conviene conocer antes de decidir si sirve para tu caso.
Lista actualizada a 0.0.6:

- **`symv`/`trmv`/`trsv`/`symm`/`trmm`/`trsm`/`syrk` (y sus variantes
  `float32`) son escalares puras, sin SIMD dedicado en ninguna
  arquitectura.** El patrón de acceso de `symv`/`symm` (decidir por
  elemento si leer `A[i][j]` o el espejado `A[j][i]`, según en qué
  triángulo cae) y la dependencia secuencial estricta de `trsv`/
  `trsm` (cada elemento necesita el resultado del anterior ya
  calculado) no se adaptan directo al mismo enfoque de bloque fijo +
  acumulador en registros que usa el resto del proyecto —
  vectorizarlas bien necesitaría separar el cálculo en una parte
  vectorizable + una parte con máscara/dependencia, que no se
  intentó en esta versión.
- **`symm`/`trmm`/`trsm`/`syrk` (y sus variantes `float32`) no
  pudieron validarse bajo NEON/QEMU en este proyecto** — viven en
  `bla_gemm.c`/`bla_gemm_f32.c`, que desde 0.0.6 requieren OpenMP
  obligatorio, y el toolchain bare-metal usado para validar NEON no
  soporta OpenMP (mismo motivo exacto que afecta a `gemm`/`gemm_tn`/
  `gemm_nt` desde esa versión, ver "Novedades de 0.0.6"). Solo se
  validaron en las cuatro rutas SIMD x86 en este sandbox; ver
  `tests/VALIDAR_EN_ANDROID.md` para validarlas en hardware ARM64
  real vía Termux.
- **`iamax`/`rot` (y sus variantes `float32`) tampoco tienen SIMD
  dedicado.** Encontrar un índice (no solo un valor agregado) o
  necesitar dos valores por posición simultáneamente es más complejo
  de vectorizar que las reducciones de suma que ya existían
  (`dot`/`asum`/`nrm2`) — candidatas para una versión futura si se
  mide que vale la pena.
- **`float32` ya cubre lo mismo que `double` en todos los niveles
  desde 0.0.6.**
- **Microkernel vectorizado a mano: AVX-512 (8×16 desde 0.0.6, antes
  8×8), AVX2 (8×8), y NEON/ARM64 (8×8) — para `gemm`/`gemm_tn`/
  `gemm_nt` y sus equivalentes `float32`. SSE2 sigue cayendo al
  escalar** en esas mismas funciones — impacto no medido todavía (es
  un CPU cada vez más raro de encontrar sin AVX2 también disponible,
  así que la prioridad de escribirlo es baja, pero la limitación es
  real, no se afirma que esté cubierto).
- **NEON (ARM64) cubre nivel 1, nivel 3, y (indirectamente, vía
  `dot`/`axpy`) `gemv`/`ger` de nivel 2, en `double` y `float32`.**
  `symv`/`trmv`/`trsv` son escalares puras en todas las
  arquitecturas, ARM64 incluido (ver punto de arriba) — no es una
  limitación específica de NEON, es que ninguna arquitectura las
  tiene vectorizadas todavía. `gemv`/`ger` no tienen un microkernel
  NEON *dedicado* (se investigó un prototipo de 2 filas para `gemv`,
  correcto pero sin poder confirmar si ayuda por el mismo motivo de
  velocidad no medible en ARM que se explica en el punto siguiente) —
  pero SÍ corren sobre NEON en la práctica, porque llaman a
  `dot`/`axpy`, que sí lo tienen.
- **La velocidad de NEON en hardware ARM real está parcialmente
  medida desde 0.0.6.** El camino NEON se desarrolló y validó por
  *correctitud* con un cross-compiler + QEMU (emulación de
  instrucción por instrucción, sin modelar ningún chip real) — miles
  de combinaciones de test en total, 0 fallos. QEMU no puede dar un
  número de velocidad confiable en ningún chip real, pero **el
  multithreading de 0.0.6 sí se confirmó por correctitud en hardware
  ARM64 real** (un teléfono Android vía Termux — 80 corridas
  repetidas, 0 divergencias, ver "Novedades de 0.0.6") — la primera
  vez que este proyecto corrió algo nativo en un chip ARM real, no
  solo bajo QEMU. Sigue pendiente: medir *velocidad* real de NEON
  (correctitud confirmada no es lo mismo que rendimiento confirmado),
  y confirmar en hardware real si un microkernel NEON dedicado para
  `gemv`/`ger` (o el ensanchamiento del microkernel de `gemm`, hecho
  para AVX-512 en 0.0.6) también ayudaría en ARM. Ver
  `tests/VALIDAR_EN_ANDROID.md`, actualizado en 0.0.6 con
  instrucciones para todo esto.
- **Multithreading limitado a `gemm`/`gemm_tn`/`gemm_nt` (ambas
  precisiones).** El resto de las funciones (todo nivel 1 y nivel 2)
  son de un solo hilo — no se paralelizaron, porque su costo es
  `O(n)`/`O(n²)`, no `O(n³)` como `gemm`, así que el margen para que
  paralelizar compense el overhead de crear threads es mucho menor.
  Además, **el speedup real del multithreading en `gemm` no está
  medido en ningún hardware con varios núcleos físicos todavía** — el
  entorno de desarrollo de este proyecto tiene 1 sola CPU (confirmado
  con múltiples métodos independientes), así que solo pudo validar
  *correctitud* bajo paralelismo simulado (más threads que núcleos,
  vía time-slicing del sistema operativo) y, para ARM64, en un
  teléfono real (también solo correctitud, sin medir tiempo). Sigue
  siendo la pregunta abierta de mayor impacto potencial para cerrar
  la brecha contra NumPy/OpenBLAS (que sí es multithreaded de forma
  confirmada).
- **OpenMP es una dependencia obligatoria desde 0.0.6, sin
  fallback.** Si `libgomp`/`omp.h` no están disponibles, el proyecto
  entero no compila — no hay una versión sin OpenMP a la que caer.
  Decisión consciente de simplicidad sobre compatibilidad máxima; ver
  "Instalar" para el paquete a instalar si esto falla.
- **Sin autotuning de `BLA_MC`/`BLA_NC`/`BLA_KC`.** Los valores son
  fijos, elegidos por razonamiento general sobre tamaños típicos de
  cache L1/L2 (con una investigación real en 0.0.4 que confirmó que
  variarlos no daba señal clara, ver historial de versiones
  anteriores), no autotuneados por script para ningún hardware
  específico.
- **Matrices muy chicas en un loop repetido:** el benchmark actual
  (ver "Comparación contra NumPy") muestra que justamente en matrices
  chicas (`n=32`, `n=64`) es donde BLApy pierde más contra NumPy
  (0.35x-0.50x) — el costo fijo de armar los buffers de packing pesa
  proporcionalmente más cuando hay poco trabajo real que hacer. No
  hay un camino especial "sin packing para matrices chicas" en esta
  versión.
- **`nrm2`/`snrm2` sin escalado para extremos.** BLAS de referencia
  (Netlib) escala internamente para evitar overflow/underflow
  prematuro en valores extremos; acá sigue siendo `sqrt(dot(x,x))`
  directo — un vector con elementos muy grandes puede dar `Inf` de
  forma prematura por overflow en el cuadrado intermedio, aunque el
  resultado final fuera representable. El límite es proporcionalmente
  más estrecho en `float32` (rango ~1e38) que en `double` (~1e308).
- **Manejo de errores mínimo.** Si `malloc` falla dentro de `gemm`
  (buffers de packing, en cualquier precisión), la función
  simplemente vuelve sin escribir nada en `C` — no hay un código de
  error ni un mecanismo de excepción. Aceptable en esta etapa
  temprana, no para una librería que alguien más vaya a integrar a
  ciegas.
- **`trsv`/`strsv` no validan si A es singular.** Un elemento cero en
  la diagonal produce división por cero (Inf/NaN de IEEE 754,
  propagándose de forma visible en el resultado) en vez de un error
  explícito — comportamiento aceptado por ahora, no ideal para uso a
  ciegas.
- **`symv`/`trmv`/`trsv`/`symm`/`trmm`/`trsm`/`syrk` no usan el
  formato "packed" real de BLAS.** `A` sigue siendo una matriz `n×n`
  completa (row-major); `uplo` controla qué mitad se *lee*, no cuánta
  memoria ocupa el array — simplificación deliberada para reducir
  riesgo de bugs de indexado, a costa de no tener el ahorro de
  memoria real de BLAS empaquetado (`n(n+1)/2` en vez de `n²`).
- **No hay column-major.** Todo asume row-major (`C` order); no existe
  un parámetro tipo `CBLAS_ORDER` para elegir layout.
- **`bla_sgemm` no tiene benchmark de velocidad forzado por SIMD
  todavía.** El benchmark contra NumPy (ver más arriba) sí cubre
  `float32`, pero no hay un `bench-forced-f32` análogo a
  `bench-forced` (que fuerza cada nivel SIMD x86 sin necesitar
  hardware distinto) — solo existe para `double`.

## Licencia

Ver `LICENSE`.
