{"id":1214,"date":"2023-07-11T18:41:10","date_gmt":"2023-07-11T16:41:10","guid":{"rendered":"https:\/\/intsight.com\/?p=1214"},"modified":"2024-03-12T14:31:14","modified_gmt":"2024-03-12T13:31:14","slug":"reuniendo-piezas-de-un-vector","status":"publish","type":"post","link":"https:\/\/intsight.com\/index.php\/2023\/07\/11\/reuniendo-piezas-de-un-vector\/","title":{"rendered":"Reuniendo piezas de un vector"},"content":{"rendered":"

Cuando hay que utilizar<\/span> 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\u00e1cil darse cuenta de que el siguiente bucle (un producto escalar de dos vectores) es un candidato ideal para usar instrucciones SIMD:<\/p>\n

\ndouble sum = 0.0;\nfor (int i = 0; i < size; i++)\n    sum += p[i] * q[i];<\/pre>\n
Hay dos zonas de memoria, p<\/code> y q<\/code>, y el \u00edndice se mueve ascendentemente de uno en uno, en ambos casos. Por lo tanto, para convertir este c\u00f3digo en instrucciones vectoriales de AVX, hacemos que cada paso del c\u00f3digo cargue cuatro entradas de cada vector, en cada paso del bucle:<\/p>\n
\ndouble sum = 0;\nint i = 0;\nif (Fma.IsSupported)\n{\n    var acc = Vector256<double>.Zero;\n    for (; i < size - 4; i += 4)\n        acc = Fma.MultiplyAdd(\n            Avx.LoadVector256(p + i),\n            Avx.LoadVector256(q + i),\n            acc);\n    acc = Avx.HorizontalAdd(acc, acc);\n    sum = acc.ToScalar() + acc.GetElement(2);\n}\nfor (; i < size; i++)\n    sum += p[i] * q[i];<\/pre>\n
Gracias a este sencillo truco, se aprovecha mejor la cach\u00e9 de memoria del procesador. Es cierto que p<\/code> y q<\/code> pueden estar en direcciones distantes, pero al menos la carga de cada trozo de cuatro valores no provoca invalidaciones de la cach\u00e9. En cambio, este otro fragmento de c\u00f3digo, extra\u00eddo del c\u00f3digo de Austra<\/a> que realiza la descomposici\u00f3n de una matriz seg\u00fan la factorizaci\u00f3n LU<\/a>, no cumple con esta regla:<\/p>\n
\ndouble m = c[k];\nfor (int i = k + 1; i < size; i++)\n    c[i] -= m * a[i * size + k];<\/pre>\n
El almacenamiento en c<\/code> es consecutivo, pero si reunimos cuatro pasos del bucle, las cargas desde el vector a<\/code> van a venir de zonas dispersas de la memoria: esto es porque la variable de bucle en el indexador est\u00e1 multiplicada por size<\/code>. Ser\u00eda una pena no poder utilizar AVX aqu\u00ed, porque podr\u00edamos reducir por cuatro las restas y multiplicaciones del bucle.<\/p>\n
Por fortuna, Intel agreg\u00f3 al conjunto de instrucciones de AVX2 una instrucci\u00f3n que, precisamente, sirve para reunir en una misma carga fragmentos de memoria no consecutivos. Se trata de la instrucci\u00f3n vgatherdpd<\/a>, que .NET implementa<\/a> mediante el m\u00e9todo Avx2.GatherVector<\/code>:<\/p>\n
\npublic static Vector256<double> GatherVector256(\n    double* baseAddress,\n    Vector128<int> index,\n    byte scale);<\/pre>\n
A primera vista, parece arameo antiguo. A segunda y tercera vistas, se asemeja m\u00e1s al sumerio de chabolas. Pero no es tan complicado explicar c\u00f3mo funciona, sobre todo si lo acompa\u00f1amos con un ejemplo. A este m\u00e9todo le pasamos un puntero a un array<\/em> de n\u00fameros de doble precisi\u00f3n, esa es la parte obvia. Adem\u00e1s, le pasamos cuatro \u00edndices de tipo entero, en el par\u00e1metro de tipo Vector128<int><\/code>. Como este tipo contiene 128 bits, es evidente que ah\u00ed podemos pasar cuatro enteros de 32 bits; ya veremos c\u00f3mo. \u00bfY el factor de escala? Aqu\u00ed es un poco redundante, y Microsoft podr\u00eda habernos ahorrado el trabajo: cuando se trabaja con vectores de double<\/code>, tenemos que pasar un 8, o lo que es igual, sizeof(double)<\/code>.<\/p>\n
Juntando las piezas, he aqu\u00ed como quedar\u00eda el c\u00f3digo transformado del fragmento original:<\/p>\n
\ndouble m = c[k];\nint i = k + 1;\nif (Avx2.IsSupported)\n{\n    var vm = Vector256.Create(m);\n    var vx = Vector128.Create(0, size, 2 * size, 3 * size);\n    for (double* p = a + i * size + k;\n        i < size - 4;\n        i += 4, p += size * 4)\n        Avx.Store(c + i,\n            Avx.Subtract(Avx.LoadVector256(c + i),\n            Avx.Multiply(Avx2.GatherVector256(p, vx, 8), vm)));\n}\nfor (; i < size; i++)\n    c[i] -= m * a[i * size + k];<\/pre>\n
Expliqu\u00e9monos: antes de entrar en el bucle, creamos un vector de n\u00fameros reales a partir del valor de la variable m<\/code>. Lo que sigue es lo que nos interesa: creamos por adelantado, un Vector128<int><\/code> con cuatro offsets o desplazamientos diferentes, m\u00faltiplos del tama\u00f1o de la fila de la matriz, que viene en size<\/code>. Y ahora viene el detalle artesanal delicado: al iniciar el bucle, creo una variable adicional de puntero, p<\/code>, que apunte i * size + k<\/code> entradas a partir de la direcci\u00f3n inicial, que est\u00e1 en a<\/code>. Recuerde que le\u00edamos el valor de a[i * size + k]<\/code> en cada paso del bucle original, y que i<\/code> es la variable de bucle. Si hubiese querido ser m\u00e1s claro, deber\u00eda haber escrito:<\/p>\n
\ndouble* p = a + (k + 1) * size + k<\/pre>\n
Pero me puede la pereza, y al fin y al cabo, en ese preciso instante, i<\/code> es equivalente a k + 1<\/code>, \u00bfno?<\/p>\n
El detalle es que, en cada paso del bucle AVX, la variable p<\/code> se incrementa en 4 * size<\/code> entradas (es decir, en cuatro filas enteras de size<\/code> entradas). Digo que es un detallazo elegante porque la alternativa obvia ser\u00eda incrementar directamente los cuatro \u00edndices enteros que guardamos en el vector vx<\/code>. \u00bfPor qu\u00e9 no lo he hecho as\u00ed? Pues porque, al trabajar con vectores de 128 bits, la instrucci\u00f3n de suma pertenecer\u00eda al conjunto de instrucciones SSE2… y resulta que el cambio entre AVX y SSE tiene un coste peque\u00f1o, pero no desde\u00f1able.<\/p>\n
… y este es el momento m\u00e1s apropiado para exclamar \u00a1pero co\u00f1o! (con perd\u00f3n), \u00bfesto no tendr\u00eda que hacerlo el compilador, o el JIT, o quien narices sea?<\/em>. Pues s\u00ed. De hecho, te\u00f3ricamente, Hotspot, el JIT de Java, tendr\u00eda que ocuparse de estas cosas. En la pr\u00e1ctica, la vectorizaci\u00f3n autom\u00e1tica en Hotspot es un tru\u00f1o como un pu\u00f1o. Y el famoso escape analysis<\/em> tampoco hace lo que se supone que tendr\u00eda que hacer. O si lo hace, no se nota. Al final, Java ha a\u00f1adido una librer\u00eda 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.<\/p>\n
NOTA: AVX2, creo recordar, se introdujo en la cuarta generaci\u00f3n de Intel, Haswell, en 2013. La primera implementaci\u00f3n en hardware de vgatherdpd<\/strong> no era para tirar cohetes, y todav\u00eda hay p\u00e1ginas y comentarios en Internet que lo recuerdan. La implementaci\u00f3n mejor\u00f3 notablemente en Skylake, la sexta generaci\u00f3n (2015). Yo tuve un Haswell durante muchos a\u00f1os, y empec\u00e9 a probar AVX2 con esa m\u00e1quina. Mi pen\u00faltimo port\u00e1til fue un Skylake. Ahora mismo tengo un PC con CPU de oncena generaci\u00f3n, Tiger Lake, que es un pepino. Lo menciono porque, por casualidad, Tiger Lake es la \u00faltima generaci\u00f3n de procesadores de \u00abconsumo\u00bb que incluye las instrucciones AVX512: con la llegada de las CPU con n\u00facleos heterog\u00e9neos, estas instrucciones se deshabilitan porque los n\u00facleos \u00abeficientes\u00bb no las implementan. Ahora mismo, si quieres AVX512, tienes que irte a Xeon, o a AMD.<\/p><\/blockquote>\n","protected":false},"excerpt":{"rendered":"
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\u00e1cil darse cuenta de que el siguiente bucle (un producto escalar de dos vectores) es un candidato ideal para usar instrucciones SIMD: double […]<\/p>\n","protected":false},"author":1,"featured_media":1213,"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":[73,19,74,20,23],"class_list":["post-1214","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-c","tag-austra","tag-linear-algebra","tag-lu-decomposition","tag-matrices","tag-simd"],"jetpack_featured_media_url":"https:\/\/i0.wp.com\/intsight.com\/wp-content\/uploads\/2023\/07\/eostre.png?fit=450%2C450&ssl=1","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1214","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=1214"}],"version-history":[{"count":28,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1214\/revisions"}],"predecessor-version":[{"id":1709,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/posts\/1214\/revisions\/1709"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media\/1213"}],"wp:attachment":[{"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/media?parent=1214"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/categories?post=1214"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/intsight.com\/index.php\/wp-json\/wp\/v2\/tags?post=1214"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}