Categorías
C#

La la, how their life goes on

Si insisto tanto en estos temas de instrucciones vectoriales, es parar preparar a mis lectores potenciales para que usen estas técnicas cuando lo crean necesario. Naturalmente, las primeras explicaciones van a ser del tipo heroico: todo hay que hacerlo a mano, razonando desde los primeros principios. Pero, como decía Molly Jones, la mujer de Desmond Jones, la la, how their life goes on: la vida sigue su curso, y lo normal en la vida de un programador es usar las técnicas que todos conocemos y apreciamos para ahorrarnos esfuerzo.

Lo primero para que todos nos ahorremos trabajo, por supuesto, consiste en usar estas cosas a través de una librería bien probada. Ya existen unas cuantas de estas, pero estoy creando Austra porque hay un nicho muy concreto, y porque la librería tiene méritos propios. La idea es hacer de Austra un proyecto open source, por supuesto. En este momento está en GitHub, pero no es por ahora un repositorio público porque la aplicación de «pruebas» está basada en WPF y DevExpress. O me compro personalmente una licencia comercial, o cambio los controles por algo gratuito. No es sencillo. Hay cosas como Avalonia o Uno Platform, que además, son multiplataformas, pero son bastante pobres en componentes. Y no quiero ni oír hablar de .NET MAUI: es un fracaso, de momento. Antes que usar .NET MAUI, prefiero «degradar» la aplicación a Windows Forms (tengo un editor de código bastante bueno) y currarme los gráficos y las cosas que falten. Pero a eso vamos: a que Austra esté disponible gratuitamente como open source, que cualquier que quiera pueda colaborar, que haya un conjunto de tests y benchmarks exhaustivo, etc, etc.

Pero no se escribe un post para explicar lo obvio. Mi verdadero objetivo ahora es explicar un par de técnicas muy tontas que bajan el listón del heroísmo al escribir este tipo de código, y que todos conocemos, pero que nos puede dar miedo usar a primeras, porque no sabes cómo va a interferir en la generación de código de .NET 7 y .NET 8 (que tiene sus cagadas, como todo).

Primer ejemplo: resulta que mucho de estos algoritmos algebraicos calculan el producto escalar de dos zonas de memoria en muchos casos. Hay una clase Vector con su correspondiente producto escalar, e incluso un ComplexVector, pero es mejor tener una rutina que maneje directamente zonas de memoria con un tamaño predeterminado. Austra tiene una clase estática CommonMatrix que precisamente implementa estas cosas. Ésta es la última versión del método en cuestión:

/// <summary>
/// Computes the dot product of two double arrays.
/// </summary>
/// <param name="p">Pointer to the first array.</param>
/// <param name="q">Pointer to the second array.</param>
/// <param name="size">Number of items in each array.</param>
/// <returns>A sum of products.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public unsafe static double DotProduct(
    double* p, double* q, int size)
{
    double sum;
    int i = 0;
    if (Avx.IsSupported)
    {
        Vector256<double> acc = Vector256<double>.Zero;
        for (int top = size & AVX_MASK; i < top; i += 4)
            acc = acc.MultiplyAdd(p + i, q + i);
        sum = acc.Sum();
    }
    else
        sum = 0;
    for (; i < size; i++)
        sum += p[i] * q[i];
    return sum;
}

A primera vista, lo único raro aquí es el atributo que fuerza un inlining agresivo. Pero hay un par de cosillas más que no son métodos habituales de AVX. Los muestro aparte aquí, porque son métodos de extensión de la misma clase CommonMatrix:

/// <summary>Sums all the elements in a vector.</summary>
/// <param name="v">
/// A intrinsics vector with four doubles.
/// </param>
/// <returns>The total value.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static double Sum(this Vector256<double> v)
{ 
    v = Avx.HorizontalAdd(v, v);
    return v.ToScalar() + v.GetElement(2);
}

Resulta que sumar los cuatro elementos de un vector de cuatro dobles no es moco de pavo. Hay teorías que pululan por la Internet sobre cuál es la forma eficiente de lograrlo. Yo no estoy seguro aún de que la mía sea la más eficiente, pero gracias al encapsulamiento en un método separado, si descubro algo mejor, tendré que tocar el código en un único lugar. Obvio, ¿no? Pero hasta comprobar que el JIT de .NET Core no se volvía loco, no me atreví a dar este paso. Ya he comprobado que no pasan cosas raras.

Es curioso, sin embargo, que el mínimo y el máximo de un vector se calcule más eficientemente con un método como el siguiente:

[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static double Max(this Vector256<double> v)
{
    var x = Sse2.Max(v.GetLower(), v.GetUpper());
    return Math.Max(x.ToScalar(), x.GetElement(1));
}

En este caso, no parece que la transición temporal a una instrucción SSE haga mucho daño al código que la rodea. De todas maneras, observe que para la reducción final al escalar hay que ir a por el segundo elemento directamente. Cosas de Intel.

Éste es otro método que usa el código original, omitiendo esta vez los comentarios XML:

internal unsafe static Vector256<double> MultiplyAdd(
    this Vector256<double> summand,
    double* multiplicand,
    double* multiplier) =>
    Fma.IsSupported
        ? Fma.MultiplyAdd(
            Avx.LoadVector256(multiplicand),
            Avx.LoadVector256(multiplier),
            summand)
        : Avx.Add(summand, Avx.Multiply(
            Avx.LoadVector256(multiplicand),
            Avx.LoadVector256(multiplier)));

Hay un par de variantes más que usan vectores directamente en vez de direcciones, y variantes adicionales para MultiplySubtract y MultiplyAddNegated, que sons métodos FMA. Estas aparentes tonterías ahorran líneas y líneas de código. En mi caso, como tengo memoria fotográfica, prefiero que el código fuente sea lo más pequeño posible para poder recordarlo más adelante.

Categorías
C#

Reuniendo piezas de un vector

Cuando hay que utilizar instrucciones vectoriales, es siempre aconsejable que se cargue en memoria un trozo consecutivo de memoria para crear vectores. Para un programador ya entrenado en AVX/AVX2/AVX512, por ejemplo, es fácil darse cuenta de que el siguiente bucle (un producto escalar de dos vectores) es un candidato ideal para usar instrucciones SIMD:

double sum = 0.0;
for (int i = 0; i < size; i++)
    sum += p[i] * q[i];

Hay dos zonas de memoria, p y q, y el índice se mueve ascendentemente de uno en uno, en ambos casos. Por lo tanto, para convertir este código en instrucciones vectoriales de AVX, hacemos que cada paso del código cargue cuatro entradas de cada vector, en cada paso del bucle:

double sum = 0;
int i = 0;
if (Fma.IsSupported)
{
    var acc = Vector256<double>.Zero;
    for (; i < size - 4; i += 4)
        acc = Fma.MultiplyAdd(
            Avx.LoadVector256(p + i),
            Avx.LoadVector256(q + i),
            acc);
    acc = Avx.HorizontalAdd(acc, acc);
    sum = acc.ToScalar() + acc.GetElement(2);
}
for (; i < size; i++)
    sum += p[i] * q[i];

Gracias a este sencillo truco, se aprovecha mejor la caché de memoria del procesador. Es cierto que p y q pueden estar en direcciones distantes, pero al menos la carga de cada trozo de cuatro valores no provoca invalidaciones de la caché. En cambio, este otro fragmento de código, extraído del código de Austra que realiza la descomposición de una matriz según la factorización LU, no cumple con esta regla:

double m = c[k];
for (int i = k + 1; i < size; i++)
    c[i] -= m * a[i * size + k];

El almacenamiento en c es consecutivo, pero si reunimos cuatro pasos del bucle, las cargas desde el vector a van a venir de zonas dispersas de la memoria: esto es porque la variable de bucle en el indexador está multiplicada por size. Sería una pena no poder utilizar AVX aquí, porque podríamos reducir por cuatro las restas y multiplicaciones del bucle.

Por fortuna, Intel agregó al conjunto de instrucciones de AVX2 una instrucción que, precisamente, sirve para reunir en una misma carga fragmentos de memoria no consecutivos. Se trata de la instrucción vgatherdpd, que .NET implementa mediante el método Avx2.GatherVector:

public static Vector256<double> GatherVector256(
    double* baseAddress,
    Vector128<int> index,
    byte scale);

A primera vista, parece arameo antiguo. A segunda y tercera vistas, se asemeja más al sumerio de chabolas. Pero no es tan complicado explicar cómo funciona, sobre todo si lo acompañamos con un ejemplo. A este método le pasamos un puntero a un array de números de doble precisión, esa es la parte obvia. Además, le pasamos cuatro índices de tipo entero, en el parámetro de tipo Vector128<int>. Como este tipo contiene 128 bits, es evidente que ahí podemos pasar cuatro enteros de 32 bits; ya veremos cómo. ¿Y el factor de escala? Aquí es un poco redundante, y Microsoft podría habernos ahorrado el trabajo: cuando se trabaja con vectores de double, tenemos que pasar un 8, o lo que es igual, sizeof(double).

Juntando las piezas, he aquí como quedaría el código transformado del fragmento original:

double m = c[k];
int i = k + 1;
if (Avx2.IsSupported)
{
    var vm = Vector256.Create(m);
    var vx = Vector128.Create(0, size, 2 * size, 3 * size);
    for (double* p = a + i * size + k;
        i < size - 4;
        i += 4, p += size * 4)
        Avx.Store(c + i,
            Avx.Subtract(Avx.LoadVector256(c + i),
            Avx.Multiply(Avx2.GatherVector256(p, vx, 8), vm)));
}
for (; i < size; i++)
    c[i] -= m * a[i * size + k];

Expliquémonos: antes de entrar en el bucle, creamos un vector de números reales a partir del valor de la variable m. Lo que sigue es lo que nos interesa: creamos por adelantado, un Vector128<int> con cuatro offsets o desplazamientos diferentes, múltiplos del tamaño de la fila de la matriz, que viene en size. Y ahora viene el detalle artesanal delicado: al iniciar el bucle, creo una variable adicional de puntero, p, que apunte i * size + k entradas a partir de la dirección inicial, que está en a. Recuerde que leíamos el valor de a[i * size + k] en cada paso del bucle original, y que i es la variable de bucle. Si hubiese querido ser más claro, debería haber escrito:

double* p = a + (k + 1) * size + k

Pero me puede la pereza, y al fin y al cabo, en ese preciso instante, i es equivalente a k + 1, ¿no?

El detalle es que, en cada paso del bucle AVX, la variable p se incrementa en 4 * size entradas (es decir, en cuatro filas enteras de size entradas). Digo que es un detallazo elegante porque la alternativa obvia sería incrementar directamente los cuatro índices enteros que guardamos en el vector vx. ¿Por qué no lo he hecho así? Pues porque, al trabajar con vectores de 128 bits, la instrucción de suma pertenecería al conjunto de instrucciones SSE2… y resulta que el cambio entre AVX y SSE tiene un coste pequeño, pero no desdeñable.

… y este es el momento más apropiado para exclamar ¡pero coño! (con perdón), ¿esto no tendría que hacerlo el compilador, o el JIT, o quien narices sea?. Pues sí. De hecho, teóricamente, Hotspot, el JIT de Java, tendría que ocuparse de estas cosas. En la práctica, la vectorización automática en Hotspot es un truño como un puño. Y el famoso escape analysis tampoco hace lo que se supone que tendría que hacer. O si lo hace, no se nota. Al final, Java ha añadido una librería de vectores basados en instrucciones SIMD, al estilo de .NET. Los compiladores de C/C++ hacen mejor papel en estos casos, pero no son omnipotentes, y a veces tienes que echarles una mano.

NOTA: AVX2, creo recordar, se introdujo en la cuarta generación de Intel, Haswell, en 2013. La primera implementación en hardware de vgatherdpd no era para tirar cohetes, y todavía hay páginas y comentarios en Internet que lo recuerdan. La implementación mejoró notablemente en Skylake, la sexta generación (2015). Yo tuve un Haswell durante muchos años, y empecé a probar AVX2 con esa máquina. Mi penúltimo portátil fue un Skylake. Ahora mismo tengo un PC con CPU de oncena generación, Tiger Lake, que es un pepino. Lo menciono porque, por casualidad, Tiger Lake es la última generación de procesadores de «consumo» que incluye las instrucciones AVX512: con la llegada de las CPU con núcleos heterogéneos, estas instrucciones se deshabilitan porque los núcleos «eficientes» no las implementan. Ahora mismo, si quieres AVX512, tienes que irte a Xeon, o a AMD.