Metadata-Version: 2.4
Name: blaspy
Version: 0.0.5
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

`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.5

**Instalable vía `pip install .`** (o desde un wheel, ver "Compilar
e instalar" más abajo) — antes de 0.0.4 solo se podía usar clonando
el repo y compilando con `make`.

Cubre los tres niveles clásicos de BLAS en `double`, más un
subconjunto creciente en `float32` (ver más abajo):

- **Nivel 1** (vector-vector, `double`): `bla_dot`, `bla_axpy`,
  `bla_scal`, `bla_nrm2`, `bla_asum`.
- **Nivel 2** (matriz-vector, `double`): `bla_gemv`, `bla_ger`
  (producto exterior / rango 1 — añadido en 0.0.2).
- **Nivel 3** (matriz-matriz, `double`): `bla_gemm` — la función
  central de este proyecto, ver la sección de diseño más abajo — más
  `bla_gemm_tn` y `bla_gemm_nt` (variantes con A o B transpuesta sin
  copiarla — añadidas en 0.0.2, ver "Uso desde C").
- **`float32`** (prefijo `s`, convención clásica de BLAS S=single,
  D=double): `bla_sdot` (nivel 1, 0.0.3), `bla_saxpy`/`bla_sscal`
  (nivel 1, 0.0.4), `bla_sgemv` (nivel 2, 0.0.4), `bla_sgemm` (nivel
  3, 0.0.3). `bla_snrm2`/`bla_sasum`/`bla_sger`/`bla_sgemm_tn`/
  `bla_sgemm_nt` siguen sin versión `float32`, ver "Qué falta".

Todas las funciones de `double` 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, nivel 1 y `bla_gemm` usan NEON; nivel 2 todavía no.
`bla_simd_info()` devuelve el nombre del nivel detectado (`avx512`,
`avx2`, `sse2`, `neon`, o `scalar`), para poder confirmarlo sin
adivinar.

### Novedades de 0.0.4

- **Paquete PyPI instalable** (`setup.py` + `pyproject.toml`). No es
  trivial en este caso concreto: BLApy carga `libbla.so` vía `ctypes`
  en vez de exponer una extensión de Python real (sin `Python.h` de
  por medio), y con la configuración por defecto de `setuptools` eso
  genera un wheel etiquetado `py3-none-any` ("puro Python, funciona
  en cualquier plataforma") — **falso**, porque el `.so` es un
  binario compilado específico de la arquitectura/SO donde se armó.
  Instalar ese wheel en otra plataforma rompería recién al intentar
  *usar* la librería, no al instalarla, que es el peor momento para
  descubrirlo. La solución (declarar un `Extension` real para forzar
  que `setuptools` etiquete el wheel por plataforma, con un
  `build_ext` personalizado para que el archivo de salida se llame
  `libbla.so` y no `libbla.cpython-312-x86_64-linux-gnu.so`) tuvo
  tres intentos fallidos antes de andar — cada uno encontrado
  corriendo el build real, no por inspección del código: un
  `AttributeError` por orden de ejecución entre `get_ext_filename` y
  `build_extensions`; un directorio duplicado (`blapy/blapy/`); y la
  causa raíz final, que `get_ext_filename()` se usa en dos lugares
  distintos de `setuptools` con formatos de retorno incompatibles
  entre sí (confirmado leyendo el código fuente real de
  `setuptools==68.1.2`, no adivinando). Validado con `pip install .`
  de punta a punta en un entorno virtual limpio: genera
  `blapy-0.0.4-cp312-cp312-linux_x86_64.whl` (tag de plataforma
  correcto, no `any`), y el paquete instalado funciona correctamente
  importado desde fuera del directorio del proyecto.
- **`bla_saxpy`, `bla_sscal`, `bla_sgemv`** (`float32`): mismo
  patrón que sus equivalentes `double`, con SIMD dedicado en las
  cuatro rutas (AVX-512, AVX2, SSE2, NEON). Validado con las mismas
  combinaciones de tamaños de borde que sus versiones `double`.
- **Microkernels AVX2 (`double` y `float32`) y NEON (`double`)
  optimizados con acumuladores desenrollados a mano** en vez de un
  array `acc[N][M]` — encontrado revisando el ensamblado real
  generado (mismo tipo de análisis que en 0.0.2 confirmó que AVX-512
  ya estaba óptimo, esta vez aplicado por primera vez a AVX2 y NEON).
  El hallazgo: con `array acc[8][2]` (double/AVX2) o `acc[8][4]`
  (NEON), el compilador necesita más registros lógicos de los que
  hay físicos disponibles (19 contra 16 YMM en AVX2; 35 contra 32
  v0-v31 en NEON) — 3 acumuladores terminaban en la pila, con
  lectura+escritura en **cada** iteración del loop sobre `kc` (no
  solo una vez), algo más costoso que el supuesto de 0.0.3 ("los
  spills caen en L1, son baratos") sin haber revisado el ensamblado
  para confirmarlo entonces. Con los acumuladores como variables
  sueltas, el compilador logra otra asignación de registros que evita
  esos spills recurrentes:
  - **`double`/AVX2**: confirmado con evidencia sólida — 20 corridas
    en dos tandas separadas (`bla_gemm` real, `n=512`, AVX2 forzado):
    mediana bajó de 10.545ms a 10.081ms, una mejora de **~4.4%**
    consistente en dirección entre ambas tandas, sin el solapamiento
    de rangos que sí se vio en los intentos descartados de 0.0.3
    (como MR=4).
  - **`float32`/AVX2**: se aplicó el mismo cambio de código por
    consistencia, pero la medición (`bla_sgemm`, `n=512`, dos tandas
    de 10 corridas) dio resultados con **dirección opuesta entre
    tandas** — ruido, no señal. Documentado así en el código en vez
    de exagerar el resultado copiando la justificación de `double`
    sin haberla confirmado ahí.
  - **`double`/NEON**: el mismo patrón de spills se confirmó en el
    ensamblado ARM64 real (bajo el cross-compiler usado para validar
    este proyecto), y el cambio se aplicó con **correctitud**
    confirmada (224 combinaciones bajo QEMU, diferencia exacta 0.0)
    pero **sin poder medir velocidad**, mismo motivo ya documentado
    en 0.0.3 (QEMU no modela hardware real, y el toolchain bare-metal
    no expone `clock_gettime` utilizable) — queda como una
    extrapolación razonada del resultado de AVX2, pendiente de
    confirmar en hardware ARM real (ver
    `tests/VALIDAR_EN_ANDROID.md`).
- **Se investigó también el tamaño de bloque de cache (`KC`/`MC`),
  nunca probado antes, y el resultado fue no aplicar ningún
  cambio:** medido en `n=512` y `n=1024`, en el gemm completo real
  (usando los archivos de producción tal como están separados, no una
  copia de un solo archivo — un primer intento con todo en un solo
  archivo dio números 800x más chicos de lo esperado por optimización
  inter-procedimental no representativa del build real, descubierto y
  descartado antes de sacar conclusiones de ahí). Con el proyecto
  real: todas las combinaciones probadas cayeron dentro de ~3-9% de
  variación entre sí, sin que ninguna ganara de forma consistente en
  varias corridas — ruido, no señal. `BLA_KC`/`BLA_MC` quedan en sus
  valores de 0.0.1 sin cambios. Se corrigió, de paso, un error de
  razonamiento en la documentación anterior sobre qué debe caber en
  cada nivel de cache — ver "Parámetros de bloqueo" más abajo.
- **No se investigó multithreading en esta versión**: este sandbox de
  desarrollo tiene una sola CPU disponible, confirmado con cuatro
  métodos independientes (`nproc`, `/proc/cpuinfo`, afinidad de
  proceso, y una prueba directa de 4 procesos C haciendo trabajo de
  CPU real en paralelo — que tardaron ~4x más que 1 proceso solo con
  4x el trabajo, exactamente el patrón esperado de un solo núcleo, no
  el de paralelismo real). Se prefirió no implementarlo sin poder
  validarlo, mismo criterio que ya se aplicó con la velocidad de NEON.
  Sigue siendo la optimización de mayor impacto potencial pendiente
  para cerrar la brecha contra NumPy/OpenBLAS (que sí es
  multithreaded), pero necesita un entorno con más de un núcleo para
  desarrollarse con la misma disciplina que el resto del proyecto.

### Comparación contra NumPy (objetivo de 0.0.4)

Medido en este entorno de desarrollo (NumPy con backend OpenBLAS,
confirmado con `np.show_config()` — aunque el sandbox de desarrollo
solo expone 1 CPU, así que esta comparación es single-thread de
ambos lados, no representa la ventaja de multithreading que OpenBLAS
tendría en una máquina con más núcleos):

| n    | BLApy `gemm` (`double`) | NumPy (`double`) | ratio | BLApy `sgemm` (`float32`) | NumPy (`float32`) | ratio |
|------|--------------------------|-------------------|-------|------------------------------|----------------------|-------|
| 128  | 0.10 ms                  | 0.08 ms           | 1.39x | 0.07 ms                      | 0.04 ms              | 1.51x |
| 256  | 0.75 ms                  | 0.57 ms           | 1.31x | 0.45 ms                      | 0.33 ms              | 1.37x |
| 512  | 6.09 ms                  | 4.22 ms           | 1.44x | 3.64 ms                      | 2.40 ms              | 1.51x |
| 1024 | 45.80 ms                 | 29.82 ms          | 1.54x | 27.60 ms                     | 18.54 ms             | 1.49x |

**Nota sobre variabilidad:** estos números vienen de la medición más
reciente de esta sesión de desarrollo, pero repetir exactamente esta
misma comparación (`double`, `n=512`) en distintos momentos de la
misma sesión dio ratios entre 1.38x y 3.01x — variación real del
entorno del sandbox compartido (probablemente carga del host, no algo
controlable desde el código), no del código en sí. La tendencia
consistente en todas las mediciones es: NumPy/OpenBLAS más rápido en
todo el rango, con la brecha en `double` creciendo con `n` — coherente
con que OpenBLAS explota mejor la jerarquía de cache en matrices
grandes gracias a autotuning específico por CPU que este proyecto no
tiene (ver "Qué falta"). En `float32` la brecha se mantuvo más
estable en las distintas mediciones (~1.4x-1.5x) — no se investigó
por qué específicamente, es una observación, no una conclusión
establecida. Quien reproduzca esta tabla en su propia máquina
probablemente obtenga números absolutos distintos; lo que debería
sostenerse es la dirección (NumPy más rápido) y la tendencia
(la brecha en `double` crece con `n`).

**Sobre el objetivo de "velocidades parecidas":** esta versión sí
encontró y aplicó una mejora real de microkernel (ver "Novedades de
0.0.4" arriba: acumuladores desenrollados en AVX2/NEON, ~4.4%
confirmado en `double`/AVX2 con evidencia sólida de 20 corridas en
dos tandas), a diferencia de las investigaciones de 0.0.2/0.0.3 que
no habían encontrado nada aplicable. Sigue sin cerrar la brecha
completa contra NumPy/OpenBLAS — el camino de mayor impacto potencial
que queda es multithreading (no investigado en esta versión por falta
de un entorno con más de un núcleo para desarrollarlo con el mismo
rigor que el resto del proyecto).


## 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 de
   `BLA_MR × BLA_NR` (8×8 con los valores por defecto) del resultado
   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 8×8 vectorizado a
mano en AVX-512, AVX2 (0.0.2), y NEON/ARM64 (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"

// 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 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 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}
```

Compilar y enlazar: `gcc tu_programa.c -L. -lbla -lm` (o enlazar
`libbla.so` directo si no está instalada como librería del sistema).

## Uso desde Python

```python
import blapy

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

blapy.dot([1, 2, 3], [4, 5, 6])  # 32.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)

# float32: mismo patrón, prefijo "s"
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)
```

`blapy.ger`, `blapy.gemm_tn`, y `blapy.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).

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):

```bash
pip install .
```

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).

**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.

## 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 completa 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)
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`)

| Parámetro | Valor | Qué es |
|-----------|-------|--------|
| `BLA_MR`  | 8     | Filas de A en el microkernel (viven en registros) |
| `BLA_NR`  | 8     | Columnas de B en el microkernel (viven en registros) |
| `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 |

Estos son valores fijos razonables para CPUs x86_64 modernos típicos
— **no** se autotunearon para ninguna máquina puntual (ver "Qué
falta"). `BLA_MR`/`BLA_NR`=8 coincide exactamente con el ancho de un
registro AVX-512 (8 doubles) — un registro por fila del bloque. Con
AVX2 (4 doubles por registro YMM) y NEON (2 doubles por registro de
128 bits), el mismo bloque lógico 8×8 se cubre con 2 y 4 registros
por fila respectivamente — más registros en total (16 en AVX2, 32 en
NEON, saturando lo disponible en ambos casos), pero con acumulador
vectorial dedicado igual, no el camino escalar. Solo SSE2 sigue
cayendo al escalar.

### `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)

Cosas que 0.0.3 deliberadamente no cubre, para que quede claro qué
esperar y qué no:

- **`float32` cubre solo `bla_sdot` y `bla_sgemm`.** El resto de
  nivel 1 (`axpy`/`scal`/`nrm2`/`asum`), nivel 2 completo
  (`gemv`/`ger`), y las variantes transpuestas de nivel 3
  (`gemm_tn`/`gemm_nt`) no tienen versión `float32` todavía —
  quedaron fuera de 0.0.3 a propósito (se cubrió lo que ya tenía SIMD
  dedicado en `double`, ver "Novedades de 0.0.3"), no es un olvido.
- **Microkernel vectorizado a mano: AVX-512, AVX2, y NEON (ARM64).
  SSE2 sigue cayendo al escalar.** Resuelto en 0.0.2 para AVX2 y NEON
  — en 0.0.1, ambos caían al fallback escalar genérico (~6x más
  lento que AVX-512, y perdía contra `numerx.matmul` en esa
  situación, ver benchmark arriba). SSE2 es el único nivel x86 que
  sigue sin camino dedicado — 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 y gemm (nivel 3) en `double`, NO nivel
  2, y NO nivel 1/3 en `float32`.** Resuelto en 0.0.3 para nivel 1
  (`bla_dot`/`bla_axpy`/`bla_scal`/`bla_asum`) — en 0.0.2 corrían con
  el camino escalar de siempre en ARM64. `bla_gemv`/`bla_ger` (nivel
  2) siguen sin camino NEON dedicado; y `bla_sdot`/`bla_sgemm`
  (`float32`) usan NEON solo en `bla_sgemm` (nivel 3), no hay
  `float32` de nivel 1 en ninguna arquitectura todavía (ver punto de
  arriba).
- **La velocidad de NEON en hardware ARM real no está medida desde
  este proyecto.** 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) — mismas
  224+176+60 combinaciones de test que x86_64, 0 fallos. Pero QEMU no
  puede dar un número de velocidad confiable en ningún chip real
  (Snapdragon, Apple Silicon, etc.). Ver `tests/VALIDAR_EN_ANDROID.md`
  para correr el mismo proyecto nativo en un teléfono real y obtener
  ese número — pendiente de que alguien con ese hardware lo corra.
- **Se investigó optimizar los microkernels de AVX-512 (0.0.2), AVX2
  y NEON (0.0.3), y en las tres arquitecturas el resultado fue no
  aplicar ningún cambio.** AVX-512 (0.0.2): prefetch manual empeoró
  ~110% en aislamiento, desenrollado sin mejora medible — revisando
  el ensamblado se confirmó que GCC ya fusiona el broadcast de A en
  `vfmadd231pd` (modo `{1to8}`), sin margen para mejorar a mano. AVX2
  (0.0.3): prefetch, desenrollado, y MR=4 mostraron mejora en el
  microkernel **aislado** (hasta ~16%), pero se revirtieron casi por
  completo o se invirtieron al medir en un **gemm completo** (MR=4:
  de ~15% a ruido entre -6%/+3%; prefetch: de ~2% a **-20%**
  consistente) — ver "Novedades de 0.0.3" para el detalle. NEON
  (0.0.3): no se pudo medir tiempo de ejecución en absoluto con el
  entorno de desarrollo de este proyecto (el cross-compiler
  bare-metal no expone `clock_gettime`/`clock()` de forma utilizable,
  y QEMU no modela hardware real de todos modos), así que no se
  aplicó ningún cambio sin poder confirmar que ayuda — la lección de
  AVX2 (que una mejora aislada se revierte en el sistema real) hace
  que extrapolar sin medir sea especialmente arriesgado acá.
- **Sin multithreading.** `bla_gemm` es de un solo hilo. BLIS/OpenBLAS
  paralelizan varios de los bucles de bloqueo con OpenMP; acá no hay
  ningún `#pragma omp` todavía.
- **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, no medidos exhaustivamente en ningún hardware
  específico ni ajustados por script.
- **Matrices muy chicas en un loop repetido:** para una sola llamada,
  el overhead de armar los buffers de packing es marginal incluso en
  matrices chicas (ver benchmark, n=32 ya gana). Pero si el caso de
  uso es MUCHAS matrices diminutas (por ejemplo 4x4) en un loop
  ajustado, ese overhead por-llamada podría no compensarse contra un
  camino sin packing en absoluto — no hay un camino especial
  "sin packing para matrices chicas" en esta versión, ni se midió ese
  patrón de uso específico.
- **`bla_nrm2` sin escalado para extremos.** BLAS de referencia
  (Netlib) escala internamente para evitar overflow/underflow
  prematuro en valores extremos; `bla_nrm2` sigue siendo
  `sqrt(dot(x,x))` directo — un vector con elementos mayores a
  ~1e150 puede dar `Inf` de forma prematura por overflow en el
  cuadrado intermedio, aunque el resultado final fuera representable
  en `double`.
- **Manejo de errores mínimo.** Si `malloc` falla dentro de
  `bla_gemm` o `bla_sgemm` (buffers de packing), 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.
- **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 todavía.** La tabla
  de benchmark contra `numerx.matmul` (más arriba) es solo para
  `bla_gemm` (`double`) — no hay una comparación equivalente para
  `bla_sgemm` (`float32`) contra ninguna otra librería, ni un
  `bench-forced-f32` análogo a `bench-forced`. La correctitud de
  `bla_sgemm` está validada con las mismas 224 combinaciones en las
  cuatro rutas SIMD (ver "Novedades de 0.0.3"), pero su velocidad
  real no se comparó contra nada externo.

## Licencia

Ver `LICENSE`.
