mdspan: Vistas de arrays multidimensionales

La propuesta de estandarización P0009 para C++23 [1], cuya implementación puede encontrarse en la referencia [2], permite la adecuada gestión en C++ moderno de vistas no-propietarias y potencialmente mutables de arrays multidimensionales. Estas vistas son una herramienta de enorme interés en campos tan diversos como la Física, la Matemática o la Ingeniería y constituyen la piedra angular en torno a la cual el comité ISO C++ pretende diseñar una biblioteca de álgebra lineal basada en BLAS [3].

Sea I un espacio de índices multidimensionales de rango R, definido como el producto cartesiano de intervalos semiabiertos [0, N₀) ⨯ [0, N₁) ⨯ ... ⨯ [0, N_R-1), con N_k natural para cada k = 0, ..., R-1. Una vista no propietaria de un array multidimensional de extensión k-ésima igual a N_k asociará a cada R-tupla de índices de acceso i ∈ I una referencia a un elemento accesible a través de un rango contiguo de índices enteros.

En el caso que nos ocupa, tales vistas se obtendrán a través de la plantilla de clase std::experimental::mdspan, que cuenta con los siguientes cuatro parámetros:

   namespace std::experimental {       template<typename Element_type, typename Extents, typename Layout_policy = layout_right,                typename Accessor_policy = default_accessor<Element_type>>       class mdspan;    }

• Element_type: el tipo de dato referenciado.
• Extents: una especialización de la plantilla de clase variádica std::experimental::extents, cuyo número de parámetros coincide con el número de extensiones (rango R) del array multidimensional, permitiendo especificar sus longitudes N_k de forma estática o dinámica (véanse los ejemplos más abajo).
• Layout_policy (opcional): especifica la fórmula (y las propiedades de la fórmula) que mapea un índice multidimensional i ∈ I a un elemento del array. La biblioteca emplea por defecto el orden de fila-principal propio de C y C++ (layout_right, row-major),  si bien proporciona también el orden de columna-principal característico de Fortran y MATLAB (layout_left, column-major), así como una generalización de los órdenes anteriores que permite registrar un paso diferente (potencialmente distinto a la unidad) para cada extensión (layout_stride).
• Accessor_policy (opcional): el descriptor de acceso que gobierna cómo se leen y se escriben los elementos (por ejemplo, de forma atómica).

A modo de ejemplo, consideremos la matriz

con R = 2 y extensiones N₀ = 3 y N₁ = 4. Adoptando un orden de fila-principal, podríamos almacenar las entradas numéricas de la matriz de forma contigua en memoria mediante un contenedor de tipo std::array<int,12>:

   auto data = std::array{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};

Resaltemos la elevada localidad espacial de las entradas así conseguida. Resulta inmediato generar una vista a un array 2-dimensional 3⨯4 a través de data. Pudiendo optar por establecer las extensiones N₀ y N₁ en tiempo de compilación, codificaríamos:

   namespace stde = std::experimental;    auto ms = stde::mdspan{data.begin(), stde::extents<3, 4>{}};

siendo ms de tipo stde::mdspan<int,stde::extents<3ull,4ull>,stde::layout_right, stde::default_accessor<int>>. La figura inferior proporciona un esquema tanto de los datos almacenados en memoria virtual como de la vista para ellos obtenida:

Resaltemos la separación estricta entre el espacio de almacenamiento de los datos y la vista mdspan, esta última no-propietaria.

En el ejemplo considerado, las llamadas a las funciones miembro públicas ms.rank(), ms.extent(0) y ms.extent(1) devolverán, respectivamente, el rango R = 2 y las extensiones N₀ = 3 (número de filas) y N₁ = 4 (número de columnas) del array multidimensional. La matriz podría ser entonces representada en la terminal mediante el bucle siguiente:

   namespace stdv = std::views;        for (auto&& r : stdv::iota(0uz, ms.extent(0))) {       for (auto&& c : stdv::iota(0uz, ms.extent(1))) {          fmt::print("{:>3}", ms(r, c));       }       fmt::print("\n");    }    // output: 1 2 3 4 \n 5 6 7 8 \n 9 10 11 12 \n

donde hemos empleado rangos-factoría std::views::iota para generar secuencias de índices 1-dimensionales para cada extensión (todas ellas empiezan en cero), así como el sufijo uz propio de C++23 para inicializar valores de tipo std::size_t. Si deseáramos, por ejemplo, multiplicar por 2 todos los elementos de la segunda fila de la matriz, codificaríamos simplemente:

   for (auto&& c : stdv::iota(0uz, ms.extent(1)))       ms(1, c) *= 2;

Podemos también obtener una vista para un subconjunto de un objeto mdspan ya existente. Ello es posible gracias a la función de slicing submdspan:

   // filas [1, 3) y columnas [1, 4) de la matriz ms tras multiplicar su segunda fila por 2:    auto sbs = stde::submdspan(ms, std::pair{1, 3}, std::pair{2, 4});        for (auto&& i : stdv::iota(0uz, sbs.extent(0))) {       for (auto&& j : stdv::iota(0uz, sbs.extent(1))) {          fmt::print("{:>3}", sbs(i, j));       }       fmt::print("\n");    }    // output: 14 16 \n 11 12 \n

donde sbs posee un memory layout de tipo layout_stride. El especificador de slice empleado para cada extensión en este caso viene dado por un rango contiguo de índices limitado por una pareja de valores de tipo std::pair (alternativamente, podría emplearse std::tuple). El rango de la vista no se ve reducido en ninguno de ambos casos. La biblioteca permite también que un especificador se reduzca a un único índice entero (en cuyo caso el rango de la vista se ve reducido en una unidad) o bien que éste sea el objeto inline y constexpr stde::full_extent (todos los índices de la extensión son entonces seleccionados, sin reducción del rango de la vista) [4]:

Especificador de slice	Argumento de submdspan	Reduce el rango
Único índice	Entero	Sí
Rango de índices	std::pair o std::tuple con dos enteros	No
Todos los índices	std::experimental::full_extent	No

---

Una de las principales características de la propuesta P0009 consiste en la posibilidad de emplear, e incluso combinar, tanto extensiones dinámicas como estáticas. Siguiendo el ejemplo proporcionado en el seminario [4] impartido por Bryce Adelstein Lelbach (Nvidia), consideremos el cálculo paralelizado de la matriz transpuesta de una matriz A cuyas extensiones m⨯n sean establecidas en runtime. A efectos puramente ilustrativos, rellenaremos las entradas de la matriz con valores reales pseudo-aleatorios:

   auto rng = std::mt19937{std::random_device{}()};    auto dst = std::uniform_real_distribution{1.0, 10.0};    auto rand = [&]{ return dst(rng); };    // extensiones del array 2-dimensional:    auto const m = /* establecida en runtime */,               n = /* establecida en runtime */;    // espacio de almacenamiento de la matriz A:    auto a_storage = std::make_unique_for_overwrite<double[]>(m*n);    // establecemos entradas aleatorias para A:    for (auto&& e : stdv::iota(0uz, m*n))       a_storage[e] = rand();    // vista para A (extensiones de longitudes establecidas dinámicamente):    auto a = stde::mdspan{a_storage.get(), m, n};

donde la vista generada es de tipo stde::mdspan<double,stde::dextents<2>, stde::layout_right,stde::default_accessor<double>>, siendo stde::dextents<2> un mero alias para stde::extents<stde::dynamic_extent, stde::dynamic_extent>. La matriz transpuesta B = A^t se calcularía entonces como:

   // espacio de almacenamiento de la matriz transpuesta B = A^t:    auto b_storage = std::make_unique_for_overwrite<double[]>(n*m);    // vista asociada:    auto b = stde::mdspan{b_storage.get(), n, m};    // espacio de índices multidimensionales de acceso a la matriz B:    auto b_indices = tl::views::cartesian_product(stdv::iota(0uz, n), stdv::iota(0uz, m));    // cálculo paralelizado de las entradas de la matriz transpuesta:    std::for_each(       std::execution::par_unseq,       std::begin(b_indices), std::end(b_indices),       [=](auto&& ind){          auto&& [r, c] = ind;          b(r, c) = a(c, r);       }    );

donde hemos empleado la implementación de la vista de producto cartesiano proporcionada por la biblioteca de rangos [5] con el fin de definir el espacio de índices multidimensionales de la matriz B. 

Observemos por último que las vistas mdspan pueden ser capturadas por valor sin incurrir apenas en coste de espacio (especialmente cuando sus extensiones sean establecidas estáticamente), dado que éstas implementan una semántica de referencia.

---

Referencias bibliográficas:

1. Programming Language C++ LEWG - P0009r14 - http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2021/p0009r14.html
2. Reference mdspan implementation - https://github.com/kokkos/mdspan/
3. D1673R4: A free function linear algebra interface based on the BLAS - http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2021/p1673r4.html
4. C++ Standard Parallelism - Bryce Adelstein Lelbach - CppCon 2021 - https://youtu.be/LW_T2RGXego
5. tl libraries - Ranges - https://github.com/TartanLlama/ranges