Categorías
Austra

Destrucción creativa

El compilador de Austra sigue mejorando. Y la librería también, que todo hay que decirlo.

Finalmente, he podido añadir generación de valores aleatorios con AVX2, además de la que ya había implementado con AVX512. La dificultad residía en la necesidad de rotar circularmente un entero largo. AVX512 te da la instrucción. Pero es fácil ver que en AVX2 es posible hacer lo mismo con dos desplazamientos lógicos y una conjunción lógica.

El otro avance importante es la vectorización del algoritmo de Box-Muller para generar distribuciones normales. Ya teníamos un logaritmo natural para vectores de 256 y de 512 bits. Ahora añadimos una función que calcula simultáneamente el seno y el coseno para estos dos tipos de vectores. De hecho, hay una variante especial que es la que usa Box-Muller, en la que todos los valores están entre cero y dos veces pi, y es un poco más rápida. Añadí esta función, la utilicé en el código del generador normal y, finalmente, se ha usado en las clases de vectores, matrices y secuencias.

Operaciones in situ

El compilador intenta ahora optimizar expresiones como la siguiente:

vector1 + vector2 + vector3

¿Qué hay aquí para optimizar? A primera vista, parecería que esa expresión es candidata para reescribirse mediante un combinador lineal, que de hecho existe y se usa para algunas optimizaciones ya:

vec([1, 1, 1, vector1, vector2, vector3)

Se podría utilizar, pero este constructor añadiría tres multiplicaciones escalares innecesarias. Lo que sí nos ahorraríamos sería la creación de dos búferes temporales como resultado de las dos sumas de la expresión original. Lo que ha aprendido a hacer el compilador es a modificar la expresión original a este otro patrón, que es más general:

(vector1 + vector2).InplaceAdd(vector3)

InplaceAdd no se puede llamar directamente desde el lenguaje: es una instrucción peligrosa, porque sobrescribe el búfer del primer operando. El compilador tiene que detectar que el primer operando es de tipo temporal, y no va a utilizarse en el resto de la fórmula, o en el peor de los casos, en el resto de la sesión. Pero en el ejemplo mostrado, la primera suma es evidentemente una construcción temporal, independientemente del sitio de donde estemos sacando las tres variables usadas.

Naturalmente, podemos hacerlo mejor si ya sabemos qué son vector1 y sus dos amigos. Si son variables de sesión, no hay nada más que podamos hacer. Lo mismo ocurre si son parámetros de una función definida por el usuario o variables locales, introducidas por una cláusula let. Pero supongamos que la expresión original fuese la siguiente:

vec::random(10) + vector2 + vector3

En ese caso, el primer operando se crea en un constructor y sólo se usa en esa expresión. El compilador puede entonces crear el siguiente código equivalente:

vec::random(10).InplaceAdd(vector2).InplaceAdd(vector3)

No importa de dónde salgan vector2 y vector3. El búfer del primer operando puede ser reutilizado, y eso es lo que hacemos. El resultado de InplaceAdd, dicho sea de paso, también se puede sobrescribir con total seguridad.

Este tipo de optimizaciones es mejor que las haga un compilador, y no la librería. Es el compilador quien tiene toda la información sobre el uso de un operando. De momento, nos movemos sobre terreno seguro. Hay algunas operaciones que podrían optimizarse con este truco, pero de momento no lo hacemos. Por ejemplo, si una variable local sólo se utiliza una vez, no hay memoria compartida que tengamos que respetar. De momento, no contamos el número de usos de cada variable local.

Finalmente, hemos ampliado las optimizaciones que hacíamos sobre vectores reales a vectores complejos. Es un poco más de código, pero no se hace más lenta la compilación.

Se supone que la señorita de la imagen es Freya bailando sobre calaveras. Eran Kali y Shiva quienes practicaban este tipo de danza, pero estaba seguro de que la tía que la IA iba a generar para Freya iba a ser más guapa. ¿Y las calaveras? Pues están debajo. Pero los generadores de Inteligencia Artificial siguen teniendo problemas para generar manos humanas. De hecho, la imagen que utilizo para Austra la llamo en mi cabeza «la chica del eccema». Llevo meses haciendo retoques a sus manos, y todavía tienen problemas. Pero me gusta el tono del pelo, el óvalo facial y la postura. No se puede tener todo.

Categorías
C# FinTech

Secuencias

Esto es una perogrullada: en un lenguaje de programación funcional, las cosas se suelen hacer de forma diferente. Y probablemente escribir «perogrullada» sea una pedantería. Dejo el párrafo vivo por pereza, pero no me lo tenga en cuenta.

Factorial

Austra no pretende ser un lenguaje de programación funcional. Ni siquiera pretende ser «Turing-completo». En algún punto del camino, he pensado en diseñar un lenguaje «de verdad», e incluso en saltarme la maquinaria de .NET para compilarlo a código nativo directamente. Pero falta un poco para eso.

De todas maneras, hay que intentar que se puedan hacer todas las cosas posibles en Austra. Por ejemplo, ¿cómo calculas un factorial, si no tienes todavía funciones recursivas, y no quieres introducir bucles? Hasta ahora, la solución era parecida a ésta:

vec(10, i => i + 1).prod

No estamos definiendo una función, sino que estamos usando un parámetro a dedo, pero funciona. Creamos un vector de 10 entradas, lo llenamos con los números del 1 al 10, y multiplicamos todos los elementos del vector. Incluso podemos «optimizar» un poco la expresión:

vec(8, i => i + 2).prod * 2

Es absurdo dejar el 1 en el vector, y ya que tenemos diez elementos, en vez de dejarlos nueve, quito también el 2. Eso lo hago porque sé que internamente el producto mete los ocho elementos en dos registros AVX y no tengo que manejar el elemento nuevo aparte. El problema es que, de todas maneras, estamos pidiendo memoria para el vector.

En Austra 2.0 (que ya tiene soporte para .NET 8 y AVX512, dicho sea de paso), ya están implementadas las secuencias de valores reales, y puedes hacer esto otro, que es más natural:

seq(2, 10).prod

La nueva clase se llama seq, y es una copia descarada del diseño de enumerables y LINQ (o de los streams de Java, o de las «list comprehensions» de tantos otros lenguajes funcionales). Ahora tenemos un seq, pero luego vendrán iseq y cseq, para secuencias de enteros y complejos. Observe también que he acortado vector a vec, para usar cvec en vez de complexvector, y para que un futuro vector de enteros se pueda llamar ivec a secas.

Las secuencias pueden hacer casi todas las cosas que hacen los vectores, y en la mayoría de los casos, generarlas es más sencillo: con el vector, usamos una función lambda. También permiten ahorrar código porque aplican los mismos trucos que LINQ para objetos. Por ejemplo:

-- Esta es una secuencia que divide el intervalo
-- de 0 a 2*pi en 360 trozos.
-- Multiplicamos cada número por dos...
seq(0, τ, 360) * 2
-- .. pero internamente, lo transformamos en esto:
seq(0, 2τ, 360)

Recursividad eficiente

Y esto último todavía no está implementado, pero es bueno pensar en estas cosas antes de lanzarse a programarlas. ¿Cómo definiríamos una función para el factorial que fuese medianamente eficiente, pero usando recursividad? Lo que se suele hacer en un lenguaje funcional, pero con una notación al estilo Austra:

def fact(n: int) =
  let f = (n, acc: int) =>
    if n = 1 then acc else f(n - 1, acc * n) in
      f(n, 1)

El truco está en usar una función auxiliar, que al estilo Austra sería una lambda en una cláusula let. La función auxiliar es recursiva por la cola, por lo que el propio JIT de .NET la puede convertir en un bucle. La recursividad en lambdas va a ser problemática, porque estoy usando directamente Linq.Expressions para generar código, pero hay un truco sencillo, que es el que voy a usar en las primeras pruebas: declarar una variable de tipo lambda, asignarle un puntero vacío, y a continuación, asignarle una expresión lambda que ya podrá usar la variable recursivamente. La alternativa es usar un combinador Y. Es un truco que viene, precisamente, del cálculo lambda, sobre el que escribiré en algún momento. No es que sea cómodo o eficiente, pero es interesante.

De todas maneras, me he dado cuenta, revisando la documentación, que en F# hay que anunciar cuando una función va a ser recursiva. Una alternativa que tengo que pensar es generar un «método dinámico» en estos casos, con el esfuerzo adicional de generar directamente el código IL. Pero tengo ya experiencia de generar IL con Freya, y puedo reutilizar código.

Categorías
Quantum

Estado cuántico

Este blog se llama Quantum Insights porque mi intención inicial era dedicarlo a la Computación Cuántica. Me he distraído un poco con los preliminares, pero es hora ya de saltar a la materia que nos interesa. No nos hará falta saber mucha Física para aprender, pero un poco no nos vendrá mal tampoco.

Mi plan para las próximas entradas es el siguiente: primero, voy a explicar qué es el «estado cuántico», a grandes rasgos. Segunda entrada: el proceso de medición. En la tercera entrada ya veremos, entonces, cómo se define un ordenador cuántico. Entiendo que, por mucho que simplifique, siempre hay temas que necesitarán aclaraciones. Utilice los comentarios ad libitum, y si lo cree necesario, envíeme un correo electrónico. Mi cuenta es mi nombre, Ian, más el dominio de este blog.

Tres axiomas

La primera mitad de la Mecánica Cuántica se explica con estos tres axiomas:

  • El estado de un sistema cuántico se describe mediante un vector en un espacio vectorial complejo $\cal H$, dotado de un producto interior hermitiano.
  • Los observables del sistema se corresponden con operadores lineales auto-adjuntos en $\cal H$.
  • La evolución en el tiempo del estado cuántico está determinada por la ecuación de Schrödinger.

Hay un montón de términos técnicos, y por ello vamos a dedicar varias secciones a cada axioma.

Espacios de Hilbert (I)

Lo primero es ver qué es un espacio vectorial complejo. Llevamos unas cuantas entradas hablando de vectores a secas, por lo que imagino que el concepto es más o menos intuitivo. Los vectores más populares son los vectores euclidianos: tres valores reales, como en $(1.2,\,2,\,-33)$. No obstante, los vectores que nos interesan para la Mecánica Cuántica tienen dos diferencias importantes:

  1. Cada componente de estos vectores va a ser un número complejo, en vez de un número real.
  2. En el espacio euclidiano hay tres dimensiones. El estado cuántico puede tener, dependiendo del sistema que se estudie, un número diferente de dimensiones. Puede ser un número finito o infinito de dimensiones. Y cuando hay infinitas dimensiones, puede tratarse de infinito numerable o infinito no numerable. Sí: manda…

No se asuste: los estados cuánticos que se manejan en Computación Cuántica son espacios finitos, y el número de dimensiones es una potencia de dos, como en $2^N$, donde $N$ es el número de qubits del sistema. Por ejemplo, dedicaremos algún tiempo a estudiar el sistema de 1 qubit, por ser el más sencillo posible. El estado de un qubit, por lo antes dicho, se puede representar como un vector de dos dimensiones complejas, como estos ejemplos:
$$
\eqalign{(0,&\,1)\cr
(0.7071 + 0.7071i,&\, 0.7071 – 0.7071i)}
$$Estos espacios vectoriales complejos, una vez que definimos un producto interior hermitiano (no lo hemos hecho todavía) se conocen como espacios de Hilbert, si tienen una propiedad adicional: que sean espacios métricos completos. Los espacios que se utilizan en computación cuántica tienen dimensiones finitas, y cumplen automáticamente con esta regla. De ahí, la $\cal H$ caligráfica que se menciona en los axiomas.

Producto interior hermitiano (I)

Dos vectores se pueden sumar y restar entre sí, y es fácil imaginar cómo se hace: componente a componente. También se puede multiplicar un vector por un escalar, casi como con los vectores euclidianos. La diferencia está en que el escalar ahora es un número complejo. Esta multiplicación, naturalmente, también se calcula componente a componente. Por ejemplo:
$$\eqalign{
(0,\, 1)+(i,\, 0)=&(i,\, 1)\cr
i\cdot(i,\, 1)=&(-1,\, i)}
$$Los vectores euclidianos tienen una operación de multiplicación escalar entre vectores, que recibe dos vectores y devuelve un número real. Cuando el espacio vectorial es complejo, sin embargo, esta operación se complica. El problema está en cómo se define la longitud de un vector. En un espacio euclidiano, acostumbramos a calcular el producto interior del vector consigo mismo, y aplicarle entonces la raíz cuadrada:
$$
\vert\vert v \vert\vert = \sqrt{v \cdot v}
$$Si queremos que las longitudes de los vectores complejos sean reales y positivas, entonces tenemos que ajustar la definición del producto interior para que, al menos, el producto de un vector consigo mismo sea no solamente real, sino además positivo. La generalización necesaria es sencilla. Supongamos que tenemos un par de vectores, $x$ e $y$, con componentes reales. En este caso, el producto interior se define clásicamente así:
$$
\eqalign{
x =& [x_1,x_2,x_3\cdots x_n]\cr
y =& [y_1,y_2,y_3\cdots y_n]\cr
x \cdot y =& \sum_{i=1}^n{x_i y_i}
}$$Estos son los productos escalares que hemos estado calculando a gogó en entradas anteriores.

Ahora supongamos que los componentes de estos vectores son complejos. La conjugada de un número complejo se define como un segundo número complejo con la misma parte real y la negación de la parte imaginaria del número original:
$$
(a + b i)^* = a – b i
$$El producto interior de dos vectores con componentes complejos es, entonces:
$$
x \cdot y = \sum_{i=1}^n{x_i^* y_i}
$$En el caso de los vectores reales, el producto interior es simétrico. Esto es, $x \cdot y = y \cdot x$. Pero para vectores complejos, la propiedad que se cumple es la siguiente:
$$
x \cdot y = (y \cdot x)^*
$$La igualdad anterior es fácil de demostrar si nos vamos a la definición de producto interior en espacios complejos. Lo que nos importa ahora es lo que ocurre cuando se toma el producto interior de un vector complejo consigo mismo:
$$
x \cdot x = (x \cdot x)^*
$$Esto significa que el valor del producto interior de un vector consigo mismo es, a la vez, igual a su valor conjugado. Pero esto sólo puede ocurrir cuando el valor es real, o sea, cuando la parte imaginaria es cero. Podemos incluso ir más lejos, y demostrar que el producto interior complejo que acabamos de definir es siempre no negativo cuando se multiplica cualquier vector consigo mismo. Por lo tanto, podemos seguir definiendo la longitud de un vector en un espacio complejo como antes:
$$
\vert\vert x \vert\vert = \sqrt{x \cdot x}
$$

La notación de Dirac (I)

Hagamos una pequeña pausa: resulta que para calcular un producto escalar, necesitamos una versión modificada de uno de los vectores… y no del otro. Pero el operador que hemos utilizado para el producto interior (el punto) «sugiere» que se trata de un operador simétrico (y no lo es). A Paul Adrien Maurice Dirac se le ocurrió una idea: digamos que a todo espacio vectorial $\mathbb{C}^n$ le corresponde automáticamente un espacio vectorial dual, y que entre ambos espacios hay una transformación biunívoca. El dual de un vector es un vector con los mismos componentes, pero conjugados. A los elementos del espacio vectorial original los llamaremos «kets» y los escribiremos de esta manera:
$$
\vert \psi \rangle
$$A los vectores del espacio dual los llamaremos «bras», y el dual del «ket» anterior se representa así:
$$
\langle \psi \vert
$$Lo que haremos a continuación es definir el producto interior como una operación que siempre toma un primer operando de tipo «bra» y un segundo operando de tipo «ket». Si yo quiero calcular el producto interior de $\vert \phi \rangle$ con $\vert\psi \rangle$, tengo que convertir antes el primer operando en un «bra», y sólo entonces puedo obtener el producto interior:
$$
\langle \phi \vert \psi \rangle
$$La expresión anterior es un número complejo, un escalar, y se cumple la antisimetría conjugada:
$$
\langle \phi \vert \psi \rangle = \langle \psi \vert \phi \rangle ^*
$$Naturalmente, la longitud de un vector, no importa si es bra o ket, puede definirse así:
$$
\vert\vert \psi \vert\vert = \sqrt{\langle\psi\vert\psi\rangle}
$$Siendo rigurosos, esto es un poco de abuso de notación. Pero no nos supondrá problema alguno.

Observables (II)

Toca explicar qué narices es el operador lineal auto-adjunto del segundo axioma. No es complicado: si el espacio vectorial es un espacio de dimensiones finitas, como en Computación Cuántica, un operador es simplemente una matriz con componentes complejos. Si tuviésemos que tratar con espacios de infinitas dimensiones, tendríamos que hilar un poco más fino, pero no es necesario en nuestro caso.

¿Recuerda que nuestro producto interior no es simétrico? Este detalle provoca que, en general, aplicar un operador al bra y al ket tengan resultados diferentes. En términos generales:
$$
\langle Ax\vert y\rangle\neq\langle x\vert Ay\rangle
$$Si queremos mover el operador $A$ al otro lado de la barra central, tenemos que transformar el operador $A$ en su adjunto $A^{*}$:
$$
\langle Ax\vert y\rangle = \langle x\vert A^{*}y\rangle
$$A nivel de celdas, la adjunta de una matriz es una matriz transpuesta creada a partir de los valores conjugados de sus componentes. Por ejemplo:
$$
\pmatrix{1&-2-i\cr 1+i&i}^* = \pmatrix{1&1-i\cr -2+i&-i}
$$Un operador auto-adjunto, simplemente, es un operador que no cambia al calcular su adjunto: $A = A^*$. Por ejemplo:
$$
\pmatrix{0&-i\cr i&0}^* = \pmatrix{0&-i\cr i&0}
$$Por lo tanto, si volvemos al producto interior, si el operador $A$ es auto-adjunto, se cumple que:
$$
\langle Ax\vert y\rangle = \langle x\vert Ay\rangle
$$¿Qué importancia tienen los operadores auto-adjuntos? Pues que los valores propios, o eigenvalues, de un operador auto-adjunto en un espacio vectorial complejo, son siempre valores reales. Este es el dato técnico. Pasemos a la interpretación física:

  1. Un «observable» de un sistema cuántico es simplemente una propiedad física del sistema que podemos medir. Ejemplos: la posición de una partícula, la velocidad de una partícula, la orientación del espín respecto a una dirección predeterminada, etc, etc.
  2. Nosotros no vamos a medir la velocidad de una partícula en Computación Cuántica. Lo advierto para la salud mental de todos nosotros. El «observable» que vamos a manejar prácticamente siempre en un ordenador cuántico es el estado binario de sus qubits. Esto lo veremos con más detalles en el momento adecuado.
  3. Como bien dice el axioma, cada «observable» se asocia a un operador auto-adjunto.
  4. El valor del observable se obtiene mediante un proceso de «medición», que veremos en la próxima entrada.
  5. Matemáticamente, la medición consiste en escoger probabilísticamente uno de los valores propios del operador asociado al observable.
  6. Como las cantidades físicas suelen ser magnitudes reales (las anoréxicas tienen una masa compleja, con una parte real y otra imaginaria), tenemos que exigir que los operadores observables tengan esta conveniente propiedad de ser auto-adjuntos.

La ecuación de Schrödinger (III)

Ya llegamos al tercer axioma, que es donde se menciona por primera vez la ecuación de Schrödinger:$$
i\hbar{d \over dt}\vert \psi(t) \rangle=H\vert \psi(t) \rangle
$$Tengo una buena noticia: ¡no necesitaremos resolver la ecuación de Schrödinger! Al menos, mientras no tengamos que enredarnos con el hardware a muy bajo nivel, claro. Sin embargo, menciono la dichosa ecuación porque conocerla nos va a ayudar a comprender mejor las reglas de la Computación Cuántica.

A la derecha de la igualdad tenemos la derivada temporal de la función de onda. Y a la izquierda, el operador de Hamilton del sistema. En pocas palabras: se trata de una ecuación lineal. Si $\vert \psi_0\rangle$ es una solución de la ecuación, y también lo es $\vert \psi_1\rangle$, entonces cualquier combinación
$$\alpha\vert \psi_0\rangle + \beta\vert \psi_1\rangle
$$donde $\alpha$ y $\beta$ son números reales, vale también como solución.

En la próxima entrada de este blog, trataremos el modelo de medición sobre un sistema cuántico.

Categorías
C#

Entran una matriz y un vector en un bar

… y claro, al rato sale un vector «transformado».

Esta entrada no es, aunque pueda parecerlo, un ripio de la anterior. Algorítmicamente, transformar un vector con una matriz se parece mucho a una sucesión de productos escalares. Pero resulta que el producto escalar, al menos hasta AVX2, tiene su truco. Vamos a comenzar por la implementación más tonta:

public static double[] Mult(double[,] a, double[] x)
{
    int m = a.GetLength(0);
    int n = a.GetLength(1);
    double[] b = new double[m];
    for (int = 0; i < m; i++)
    {
        double d = 0;
        for (int j = 0; j < n; j++)
            d += a[i, j] * x[j];
        b[i] = d;
    }
    return b;
}

Recordemos que tenemos un «handicap» autoimpuesto por representar las matrices como arrays bidimensionales de C#. Pero esta vez no voy a dar la brasa con los punteros, que ya sabemos que resuelven este problema sin pestañear. Esta es la implementación final que necesitamos, con soporte opcional de AVX para cuando esté disponible y merezca la pena:

public static unsafe double[] Mult(double[,] a, double[] x)
{
    int m = a.GetLength(0);
    int n = a.GetLength(1);
    double[] b = new double[m];
    int lastBlockIndex = n - (n % 4);
    fixed (double* pA = a)
    fixed (double* pX = x)
    fixed (double* pB = b)
    {
        double* pA1 = pA;
        double* pB1 = pB;
        if (n >= 12 && Avx2.IsSupported)
            for (int i = 0; i < m; i++)
            {
                int j = 0;
                var v = Vector256<double>.Zero;
                while (j < lastBlockIndex)
                {
                    v = Avx.Add(
                        v,
                        Avx.Multiply(
                            Avx.LoadVector256(pA1 + j),
                            Avx.LoadVector256(pX + j)));
                    j += 4;
                }
                v = Avx.HorizontalAdd(v, v);
                double d = v.ToScalar() + v.GetElement(2);
                for (; j < n; j++)
                    d += pA1[j] * pX[j];
                *pB1 = d;
                pA1 += n;
                pB1++;
            }
        else
            for (int i = 0; i < m; i++)
            {
                int j = 0;
                double d = 0;
                while (j < lastBlockIndex)
                {
                    d += (*(pA1 + j) * *(pX + j)) +
                        (*(pA1 + j + 1) * *(pX + j + 1)) +
                        (*(pA1 + j + 2) * *(pX + j + 2)) +
                        (*(pA1 + j + 3) * *(pX + j + 3));
                    j += 4;
                }
                for (; j < n; j++)
                     d += pA1[j] * pX[j];
                *pB1 = d;
                pA1 += n;
                pB1++;
            }
    }
    return b;
}

Esta vez, el código SIMD sólo se usa cuando hay doce o más elementos en el vector. La cifra la he elegido experimentando en mi i7-4770. Puede que en otros ordenadores, el umbral sea más bajo incluso.

Tengo que explicar cómo se implementa un producto escalar con SIMD, porque no es muy evidente. Uno diría que hay que acumular un escalar en una variable global al bucle… pero no hay ninguna operación SIMD que calcule directamente la suma de las cuatro multiplicaciones necesarias. La explicación oficial es que una suma de ese tipo destrozaría el paralelismo de la CPU. Y yo me lo creo, de veras. La consecuencia es que necesitamos acumular las multiplicaciones en cuatro variables; es decir, en un vector que hace de acumulador.

Las cosas se ponen de color hormiga cuando terminamos el bucle y tenemos entonces que sumar los cuatro elementos del vector acumulador. Analicemos las líneas 27 y 28 del listado anterior. Según mis experimentos, es la forma más rápida de conseguirlo. HorizontalAdd, cuando se trata de Vector256<double>, suma el primer elemento con el segundo, y lo almacena por partida doble en el primer y segundo elemento. A la vez, suma el tercero y el cuarto y hace lo mismo para guardar el resultado. Los métodos de extensión ToScalar() y GetElement() acceden entonces directamente al primer y tercer elemento y los suma. Mantengo la llamada inicial a HorizontalAdd porque, teóricamente, puede hacer dos de las sumas en paralelo, pero puedes experimentar a ver qué pasa si accedes directamente a los cuatro elementos y los sumas como toda la vida. A mí ya se me ha acabado la partida de tiempo libre para este experimento.

La razón para la controversia es que, en realidad, Internet está lleno de recomendaciones para hacer esta suma final de esta otra manera:

v = Avx2.Permute4x64(
    Avx.HorizontalAdd(v, v),
    0b00_10_01_11);
double d = Avx.HorizontalAdd(v, v).ToScalar();
// v = Avx.HorizontalAdd(v, v);
// double d = v.ToScalar() + v.GetElement(2);

Es decir: se llama dos veces a HorizontalAdd, pasando entre medias por una permutación entre escalares. En la arquitectura Haswell, al menos, esto funciona más lento que mi solución.

Si multiplico una matriz aleatoria de 64×64 por un vector de 64 elementos, obtengo estas cifras:

Method Mean Error StdDev Median
MultVector 5.762 μs 0.1142 μs 0.2227 μs 5.646 μs
FMultVector 1.814 μs 0.0320 μs 0.0416 μs 1.818 μs

No está mal, aunque no conseguimos tanta ventaja como con la multiplicación entre matrices. La versión con punteros y sin SIMD tampoco va mal, pero queda muy claro que el SIMD acelera este código. De paso, ya tenemos un patrón de código para productos escalares (y para cosas más raras como multiplicar un vector de sensibilidad delta-gamma por un escenario histórico: cosas de la valoración de productos financieros).

Por cierto, el mejor chiste que conozco sobre gente que entra en un bar tiene que ver con la Mecánica Cuántica. Dice así: entra el Gato de Schrödinger en un bar… y no entra.