{"id":192,"date":"2020-03-23T13:11:36","date_gmt":"2020-03-23T12:11:36","guid":{"rendered":"https:\/\/intsight.com\/?p=192"},"modified":"2021-01-21T17:12:28","modified_gmt":"2021-01-21T16:12:28","slug":"entran-una-matriz-y-un-vector-en-un-bar","status":"publish","type":"post","link":"https:\/\/intsight.com\/index.php\/2020\/03\/23\/entran-una-matriz-y-un-vector-en-un-bar\/","title":{"rendered":"Entran una matriz y un vector en un bar"},"content":{"rendered":"

… y claro, al rato<\/span> sale un vector \u00abtransformado\u00bb.<\/p>\n

Esta entrada no es, aunque pueda parecerlo, un ripio de la anterior. Algor\u00edtmicamente, transformar un vector con una matriz se parece mucho a una sucesi\u00f3n de productos escalares. Pero resulta que el producto escalar, al menos hasta AVX2, tiene su truco. Vamos a comenzar por la implementaci\u00f3n m\u00e1s tonta:<\/p>\n

public static double[] Mult(double[,] a, double[] x)\n{\n    int m = a.GetLength(0);\n    int n = a.GetLength(1);\n    double[] b = new double[m];\n    for (int = 0; i < m; i++)\n    {\n        double d = 0;\n        for (int j = 0; j < n; j++)\n            d += a[i, j] * x[j];\n        b[i] = d;\n    }\n    return b;\n}<\/pre>\n
Recordemos que tenemos un \u00abhandicap\u00bb 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\u00f1ear. Esta es la implementaci\u00f3n final que necesitamos, con soporte opcional de AVX para cuando est\u00e9 disponible y merezca la pena:<\/p>\n
public static unsafe double[] Mult(double[,] a, double[] x)\n{\n    int m = a.GetLength(0);\n    int n = a.GetLength(1);\n    double[] b = new double[m];\n    int lastBlockIndex = n - (n % 4);\n    fixed (double* pA = a)\n    fixed (double* pX = x)\n    fixed (double* pB = b)\n    {\n        double* pA1 = pA;\n        double* pB1 = pB;\n        if (n >= 12 && Avx2.IsSupported)\n            for (int i = 0; i < m; i++)\n            {\n                int j = 0;\n                var v = Vector256<double>.Zero;\n                while (j < lastBlockIndex)\n                {\n                    v = Avx.Add(\n                        v,\n                        Avx.Multiply(\n                            Avx.LoadVector256(pA1 + j),\n                            Avx.LoadVector256(pX + j)));\n                    j += 4;\n                }\n                v = Avx.HorizontalAdd(v, v);\n                double d = v.ToScalar() + v.GetElement(2);\n                for (; j < n; j++)\n                    d += pA1[j] * pX[j];\n                *pB1 = d;\n                pA1 += n;\n                pB1++;\n            }\n        else\n            for (int i = 0; i < m; i++)\n            {\n                int j = 0;\n                double d = 0;\n                while (j < lastBlockIndex)\n                {\n                    d += (*(pA1 + j) * *(pX + j)) +\n                        (*(pA1 + j + 1) * *(pX + j + 1)) +\n                        (*(pA1 + j + 2) * *(pX + j + 2)) +\n                        (*(pA1 + j + 3) * *(pX + j + 3));\n                    j += 4;\n                }\n                for (; j < n; j++)\n \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 d += pA1[j] * pX[j];\n                *pB1 = d;\n                pA1 += n;\n                pB1++;\n            }\n    }\n    return b;\n}<\/pre>\n
Esta vez, el c\u00f3digo SIMD s\u00f3lo se usa cuando hay doce o m\u00e1s elementos en el vector. La cifra la he elegido experimentando en mi i7-4770. Puede que en otros ordenadores, el umbral sea m\u00e1s bajo incluso.<\/p>\n
Tengo que explicar c\u00f3mo se implementa un producto escalar con SIMD, porque no es muy evidente. Uno dir\u00eda que hay que acumular un escalar en una variable global al bucle… pero no hay ninguna operaci\u00f3n SIMD que calcule directamente la suma de las cuatro multiplicaciones necesarias. La explicaci\u00f3n oficial es que una suma de ese tipo destrozar\u00eda 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.<\/p>\n
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\u00edneas 27 y 28 del listado anterior. Seg\u00fan mis experimentos, es la forma m\u00e1s r\u00e1pida de conseguirlo. HorizontalAdd<\/code>, cuando se trata de Vector256<double><\/code>, 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\u00e9todos de extensi\u00f3n ToScalar()<\/code> y GetElement()<\/code> acceden entonces directamente al primer y tercer elemento y los suma. Mantengo la llamada inicial a HorizontalAdd<\/code> porque, te\u00f3ricamente, puede hacer dos de las sumas en paralelo, pero puedes experimentar a ver qu\u00e9 pasa si accedes directamente a los cuatro elementos y los sumas como toda la vida. A m\u00ed ya se me ha acabado la partida de tiempo libre para este experimento.<\/p>\n
La raz\u00f3n para la controversia es que, en realidad, Internet est\u00e1 lleno de recomendaciones para hacer esta suma final de esta otra manera:<\/p>\n
v = Avx2.Permute4x64(\n    Avx.HorizontalAdd(v, v),\n    0b00_10_01_11);\ndouble d = Avx.HorizontalAdd(v, v).ToScalar();\n\/\/ v = Avx.HorizontalAdd(v, v);\n\/\/ double d = v.ToScalar() + v.GetElement(2);<\/pre>\n
Es decir: se llama dos veces a HorizontalAdd<\/code>, pasando entre medias por una permutaci\u00f3n entre escalares. En la arquitectura Haswell, al menos, esto funciona m\u00e1s lento que mi soluci\u00f3n.<\/p>\n
Si multiplico una matriz aleatoria de 64×64 por un vector de 64 elementos, obtengo estas cifras:<\/p>\n\n\n\n\n\n\n
Method<\/th>\nMean<\/th>\nError<\/th>\nStdDev<\/th>\nMedian<\/th>\n<\/tr>\n<\/thead>\n
MultVector<\/td>\n5.762 \u03bcs<\/td>\n0.1142 \u03bcs<\/td>\n0.2227 \u03bcs<\/td>\n5.646 \u03bcs<\/td>\n<\/tr>\n
FMultVector<\/td>\n1.814 \u03bcs<\/td>\n0.0320 \u03bcs<\/td>\n0.0416 \u03bcs<\/td>\n1.818 \u03bcs<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n
No est\u00e1 mal, aunque no conseguimos tanta ventaja como con la multiplicaci\u00f3n entre matrices. La versi\u00f3n con punteros y sin SIMD tampoco va mal, pero queda muy claro que el SIMD acelera este c\u00f3digo. De paso, ya tenemos un patr\u00f3n de c\u00f3digo para productos escalares (y para cosas m\u00e1s raras como multiplicar un vector de sensibilidad delta-gamma por un escenario hist\u00f3rico: cosas de la valoraci\u00f3n de productos financieros).<\/p>\n
Por cierto, el mejor chiste que conozco sobre gente que entra en un bar tiene que ver con la Mec\u00e1nica Cu\u00e1ntica. Dice as\u00ed: entra el Gato de Schr\u00f6dinger en un bar… y no entra.<\/p>\n","protected":false},"excerpt":{"rendered":"
… y claro, al rato sale un vector \u00abtransformado\u00bb. Esta entrada no es, aunque pueda parecerlo, un ripio de la anterior. Algor\u00edtmicamente, transformar un vector con una matriz se parece mucho a una sucesi\u00f3n de productos escalares. Pero resulta que el producto escalar, al menos hasta AVX2, tiene su truco. Vamos a comenzar por la […]<\/p>\n","protected":false},"author":1,"featured_media":195,"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":[27,20,10,25,23,26],"class_list":["post-192","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-c","tag-dot-product","tag-matrices","tag-optimization","tag-pointers","tag-simd","tag-vector"],"jetpack_featured_media_url":"https:\/\/i0.wp.com\/intsight.com\/wp-content\/uploads\/2020\/03\/catbox.png?fit=360%2C234&ssl=1","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/192","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=192"}],"version-history":[{"count":21,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/192\/revisions"}],"predecessor-version":[{"id":571,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/192\/revisions\/571"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media\/195"}],"wp:attachment":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media?parent=192"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/categories?post=192"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/tags?post=192"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}