{"id":1338,"date":"2023-08-13T14:21:24","date_gmt":"2023-08-13T12:21:24","guid":{"rendered":"https:\/\/intsight.com\/?p=1338"},"modified":"2024-03-12T14:34:03","modified_gmt":"2024-03-12T13:34:03","slug":"el-gran-secreto-de-los-complejos","status":"publish","type":"post","link":"https:\/\/intsight.com\/index.php\/2023\/08\/13\/el-gran-secreto-de-los-complejos\/","title":{"rendered":"El Gran Secreto de los Complejos"},"content":{"rendered":"<p><span style=\"font-variant: small-caps; font-size: 107%\">Al grano:<\/span> el Gran Secreto de los N\u00fameros Complejos es que, si quieres utilizar instrucciones AVX para acelerar los c\u00e1lculos, la mejor forma de representarlos no es la que todos imaginamos: la parte real y, a continuaci\u00f3n, la parte imaginaria.<\/p>\n<p>Partamos de una regla b\u00e1sica de las instrucciones vectoriales:<\/p>\n<ul>\n<li>Es mejor manejar estructuras de <em>arrays<\/em> que <em>arrays<\/em> de estructuras.<\/li>\n<\/ul>\n<p>Observe que \u00e9sta es una p\u00edldora dif\u00edcil de tragar en la Programaci\u00f3n Orientada a Objetos. <code>Complex<\/code> es una clase que ya est\u00e1 (bien) definida en <code>System.Numerics<\/code>, pero para simplificar la explicaci\u00f3n, voy a fingir que la definimos nosotros. Con la POO en mente, comenzar\u00edamos definiendo la estructura, junto con sus m\u00e9todos, y la har\u00edamos probablemente implementar algunas interfaces, por completitud, en este plan:<\/p>\n<pre class=\"EnlighterJSRAW\">\npublic readonly struct Complex\n{\n    public double Real { get; }\n    public double Imaginary { get; }\n\n    public Complex(double re, double im) =>\n        (Real, Imaginary) = (re, im);\n\n    \/\/ Y as\u00ed, sucesivamente&hellip;\n}\n<\/pre>\n<p>Si quisi\u00e9ramos entonces un vector de n\u00fameros complejos, har\u00edamos algo parecido a esto:<\/p>\n<pre class=\"EnlighterJSRAW\">\npublic readonly struct ComplexVector\n{\n    private readonly Complex[] values;\n\n    public unsafe ComplexVector(Complex[] values) =>\n        this.values = values;\n\n    \/\/ Y as\u00ed, sucesivamente&hellip;\n}\n<\/pre>\n<p>Pues bien, ahora llego yo (o la sacrosanta Realidad, si lo prefiere hacer menos personal) y le digo que la mejor forma de programar un vector de complejos, al menos si queremos acelerarlo con AVX, es la siguiente:<\/p>\n<pre class=\"EnlighterJSRAW\">\npublic readonly struct ComplexVector\n{\n    private readonly double[] re;\n    private readonly double[] im;\n\n    \/\/ Omito verificaciones de igual longitud\n    \/\/ para simplificar el ejemplo.\n    public ComplexVector(double[] re, double[] im) =>\n        (this.re, this.im) = (re, im);\n\n    public unsafe ComplexVector(Complex[] values)\n    {\n        this.re = new double[values.Length];\n        this.im = new double[values.Length];\n        fixed (double* p = re, q = im)\n            for (int i = 0; i < values.Length; i++)\n                (p[i], q[i]) = values[i];\n    }\n\n    \/\/ Y as\u00ed, sucesivamente&hellip;\n}\n<\/pre>\n<p>Podemos dejar la estructura <code>Complex<\/code> original: nos es \u00fatil. Pero al representar la lista de complejos, es mejor que cada campo vaya en su propia lista. Es cierto tambi\u00e9n que deber\u00edamos utilizar AVX para convertir un <em>array<\/em> tradicional de complejos en un vector: existe la posibilidad, pero no lo voy a mostrar aqu\u00ed, para simplificar. He omitido tambi\u00e9n un m\u00e9todo de extensi\u00f3n que he a\u00f1adido en una clase est\u00e1tica para poder \"deconstruir\" f\u00e1cilmente un complejo en sus componentes. No tiene mucha trascendencia, pero ah\u00ed va, para que no haya tantos espacios en blanco:<\/p>\n<pre class=\"EnlighterJSRAW\">\n[MethodImpl(MethodImplOptions.AggressiveInlining)]\npublic static void Deconstruct(\n    this Complex cmp, out double re, out double im) =>\n    (re, im) = (cmp.Real, cmp.Imaginary);\n<\/pre>\n<h4>\u00bfC\u00f3mo se da cuenta uno?<\/h4>\n<p>\u00bfC\u00f3mo se da cuenta uno de estas cosas? StackOverflow est\u00e1 lleno de consejos de este tipo, escritos por personas que ya lo han sufrido en sus carnes. Pero a uno no se le enciende la bombilla hasta que se pega con su propia pared. En este caso, fue intentando acelerar una Transformada Discreta de Fourier para AUSTRA. El c\u00f3digo sin acelerar usaba los complejos como se han usado casi toda la vida: cada real con su parte imaginaria. A priori, cuando uno no conoce bien AVX, se imagina que ser\u00e1 relativamente sencillo manejar dos n\u00fameros complejos dentro de un vector de cuatro valores reales, y que el conjunto de instrucciones va a estar de tu parte. Craso error.<\/p>\n<p>Al final, tuve que limitarme a acelerar algunas partes que, oportunistamente, \"se dejaban\", casi siempre con c\u00f3digo SSE y vectores de s\u00f3lo 128 bits. Lo normal, cuando he acelerado otros algoritmos, ha sido reducir los tiempos de ejecuciones de cuatro hasta incluso ocho o diez veces. En este caso, la mejora s\u00f3lo ha sido la mitad, en la mayor\u00eda de los casos, y en algunos, la tercera parte. Como resultado, tengo pendiente replantearme todo el asunto de la Transformada Discreta de Fourier, pero utilizando listas paralelas para los reales y los imaginarios.<\/p>\n<p>Quiero que vea, de todas maneras, lo sencillo que es usar la t\u00e9cnica de <em>estructura de listas en vez de listas de estructuras<\/em>. El siguiente m\u00e9todo es el producto escalar de dos vectores complejos, representados como Dios manda:<\/p>\n<pre class=\"EnlighterJSRAW\">\npublic static unsafe Complex operator *(\n    ComplexVector v1, ComplexVector v2)\n{\n    if (v1.Length != v2.Length)\n        throw new VectorLengthException();\n    fixed (double* pr = v1.re, pi = v1.im,\n                   qr = v2.re, qi = v2.im)\n    {\n        double sumRe = 0, sumIm = 0;\n        int i = 0, size = v1.Length;\n        if (Avx.IsSupported)\n        {\n            Vector256&lt;double&gt; accRe = Vector256&lt;double&gt;.Zero;\n            Vector256&lt;double&gt; accIm = Vector256&lt;double&gt;.Zero;\n            for (int top = size & ~3; i < top; i += 4)\n            {\n                var vpr = Avx.LoadVector256(pr + i);\n                var vpi = Avx.LoadVector256(pi + i);\n                var vqr = Avx.LoadVector256(qr + i);\n                var vqi = Avx.LoadVector256(qi + i);\n                accRe = Avx.Add(accRe,\n                    Avx.Multiply(vpr, vqr)\n                       .MultiplyAdd(vpi, vqi));\n                accIm = Avx.Add(accIm,\n                    Avx.Multiply(vpi, vqr)\n                       .MultiplySub(vpr, vqi));\n            }\n            sumRe = accRe.Sum();\n            sumIm = accIm.Sum();\n        }\n        for (; i < size; i++)\n        {\n            sumRe += pr[i] * qr[i] + pi[i] * qi[i];\n            sumIm += pi[i] * qr[i] - pr[i] * qi[i];\n        }\n        return new(sumRe, sumIm);\n    }\n}\n<\/pre>\n<blockquote style=\"font-size:85%\"><p>Si no lo recuerda del \u00e1lgebra, cuando se trata de vectores complejos, el producto escalar usa la conjugada del segundo operando. Esto es: la parte imaginaria del segundo operando invierte su signo. Esto es lo que permite que el producto escalar de un vector consigo mismo sea un valor real.<\/p><\/blockquote>\n<p>El c\u00f3digo vectorial tiene una correspondencia uno a uno con el c\u00f3digo escalar que maneja la parte de los <em>arrays<\/em> que puede sobrar al final. Si repasa la entrada anterior, la del algoritmo de <a href=\"https:\/\/intsight.com\/index.php\/2023\/08\/07\/el-algoritmo-de-welford\/\" rel=\"noopener\" target=\"_blank\">Welford<\/a>, notar\u00e1 el mismo patr\u00f3n: a pesar de que el algoritmo escalar es bastante oscuro, es relativamente sencillo convertir esa parte en c\u00f3digo vectorial. La parte m\u00e1s complicada de la entrada pasada era cuando hab\u00eda que mezclar los cuatro acumuladores individuales. Es el mismo problema que hemos visto ya unas cuantas veces cuando calculamos un producto escalar: acumular es sencillo. Lo complicado es sumar despu\u00e9s los cuatro acumuladores.<\/p>\n<h4>El Misterioso Constructor Postergado<\/h4>\n<p>Para no dejar demasiados cabos sueltos, aqu\u00ed tiene una posible implementaci\u00f3n del c\u00f3digo que reparte partes reales e imaginarias a sus respectivos <em>arrays<\/em>. Mi primer impulso fue utilizar <code>Avx2.GatherVector<\/code>: leer cuatro partes reales salt\u00e1ndome las imaginarias, y luego leer cuatro partes imaginarias. Pero, por desgracia, el tiempo de ejecuci\u00f3n del constructor se dispar\u00f3 al doble. No hay forma humana de predecir estas cosas, que no sea prueba, benchmark y error.<\/p>\n<p>La versi\u00f3n que funciona, y que reduce casi a la mitad el tiempo de la versi\u00f3n de m\u00e1s arriba, lee cuatro complejos en dos vectores de 256 bits. Lo primero que hace es usar <code>Avx.Shuffle<\/code> para \"barajar las cartas\" y juntar todas las partes reales en un mismo vector de 256 bits, y las imaginarias en otro. No me pida que calcule estas cosas de memorias. Cuando me tocan estos marrones, tengo todav\u00eda que pillar un cuaderno y l\u00e1piz, e irme a <a href=\"https:\/\/www.felixcloutier.com\/x86\/shufpd\" rel=\"noopener\" target=\"_blank\">p\u00e1ginas como \u00e9sta<\/a>, para repasar los diagramas. He visto tambi\u00e9n que el <code>Shuffle<\/code> se puede sustituir tambi\u00e9n por llamadas a <code>UnpackHigh<\/code>\/<code>UnpackLoad<\/code>. Es probable que estas llamadas den mejores tiempos, pero no me ha dado tiempo a hacer la prueba.<\/p>\n<p>El problema de <code>Shuffle<\/code> (y de las alternativas mencionadas) es que te deja los n\u00fameros en el orden <code>[1, 3, 2, 4]<\/code>. Si no es importante respetar el orden, se pueden quedar as\u00ed. Pero si hay que reordenar los elementos, hay que usar <code>Avx2.Permute4x64<\/code> para ello. En general, AVX intenta, dentro de lo posible, no pasar valores de una mitad a la otra mitad del vector. Hay que usar cosas introducidas en AVX2 para conseguirlo. Por ese motivo, el constructor verifica si <code>Avx2.IsSupported<\/code> antes de lanzarse al r\u00edo:<\/p>\n<pre class=\"EnlighterJSRAW\">\npublic unsafe ComplexVector(Complex[] values)\n    : this(values.Length)\n{\n    fixed (double* p = re, q = im)\n    fixed (Complex* r = values)\n    {\n        int i = 0;\n        if (Avx2.IsSupported)\n        {\n            for (int top = values.Length & ~7; i < top; i += 4)\n            {\n                var v1 = Avx.LoadVector256((double*)(r+i));\n                var v2 = Avx.LoadVector256((double*)(r+i+2));\n                Avx.Store(p + i, Avx2.Permute4x64(\n                    Avx.Shuffle(v1, v2, 0b0000), 0b11011000));\n                Avx.Store(q + i, Avx2.Permute4x64(\n                    Avx.Shuffle(v1, v2, 0b1111), 0b11011000));\n            }\n        }\n        for (; i < values.Length; i++)\n            (p[i], q[i]) = r[i];\n    }\n}\n<\/pre>\n<p>N\u00fameros: en mi m\u00e1quina, un i9-11900K, crear un <code>ComplexVector<\/code> directamente a partir de un <em>array<\/em> de 1024 complejos, tardaba m\u00e1s o menos un milisegundo. Con las mejoras AVX2, tarda 650 microsegundos. Casi la mitad. Y lo mejor, para mi gusto, es que no he tenido que usar paralelismo con tareas. El usuario de la librer\u00eda ya usar\u00e1 ese paralelismo cuando lo considere necesario, y tendr\u00e1 las manos m\u00e1s libres.<\/p>\n<p>Como regalo, le dejo la conversi\u00f3n inversa: de vector complejo a <em>array<\/em> de complejos:<\/p>\n<pre class=\"EnlighterJSRAW\">\npublic unsafe static explicit\n    operator Complex[](ComplexVector v)\n{\n    Complex[] result = new Complex[v.Length];\n    fixed (double* p = v.re, q = v.im)\n    fixed (Complex* r = result)\n    {\n        int i = 0;\n        if (Avx2.IsSupported)\n        {\n            for (int top = v.Length & ~3; i < top; i += 4)\n            {\n                var vr = Avx.LoadVector256(p + i);\n                var vi = Avx.LoadVector256(q + i);\n                Avx.Store((double*)(r + i),\n                    Avx2.Permute4x64(Avx.Permute2x128(\n                    vr, vi, 0b0010_0000), 0b11_01_10_00));\n                Avx.Store((double*)(r + i + 2),\n                    Avx2.Permute4x64(Avx.Permute2x128(\n                    vr, vi, 0b0011_0001), 0b11_01_10_00));\n                }\n            }\n        for (; i < result.Length; i++)\n            r[i] = new(p[i], q[i]);\n    }\n    return result;\n}\n<\/pre>\n<p>El tiempo de ejecuci\u00f3n se reduce \"s\u00f3lo\" a las tres cuartas partes, pero yo creo que merece la pena.<\/p>\n","protected":false},"excerpt":{"rendered":"<p>Al grano: el Gran Secreto de los N\u00fameros Complejos es que, si quieres utilizar instrucciones AVX para acelerar los c\u00e1lculos, la mejor forma de representarlos no es la que todos imaginamos: la parte real y, a continuaci\u00f3n, la parte imaginaria. Partamos de una regla b\u00e1sica de las instrucciones vectoriales: Es mejor manejar estructuras de arrays [&hellip;]<\/p>\n","protected":false},"author":1,"featured_media":1339,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[4],"tags":[65,15,73,78,77,21,23],"class_list":["post-1338","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-c","tag-net","tag-algorithms","tag-austra","tag-avx","tag-complex-numbers","tag-eigenvalues","tag-simd"],"jetpack_featured_media_url":"https:\/\/i0.wp.com\/intsight.com\/wp-content\/uploads\/2023\/08\/complex.png?fit=490%2C324&ssl=1","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1338","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/comments?post=1338"}],"version-history":[{"count":34,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1338\/revisions"}],"predecessor-version":[{"id":1710,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1338\/revisions\/1710"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media\/1339"}],"wp:attachment":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media?parent=1338"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/categories?post=1338"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/tags?post=1338"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}