C++ многомерный массив - есть что-то новое стандартное?

Проблема с std::array не в том, что он на стеке, а в том что размер фиксированный во время компиляции. Поместить его в кучу легко: auto arr = std::make_unique<std::array<...>>();.

Еще иногда любят делать std::vector<std::vector<...>> - но это плохо, потому что в куче создается несколько блоков памяти, а не один (в вашем примере та же проблема). Лучше - один длинный одномерный std::vector, и вручную вычислять индекс в нем.

В C++23 добавили std::mdspan, чтобы не считать индексы в одномерных контейнерах руками: https://gcc.godbolt.org/z/7oqhs4K8e

#include <cstddef>
#include <iostream>
#include <mdspan>
#include <vector>

int main()
{
    std::size_t w = 3, h = 2;
    std::vector<int> storage(w * h);
    auto span = std::mdspan(storage.data(), h, w);

    // Если компилятор - MSVC, пока работает только так: span[std::array{0,0}] = 1;
    span[0,0] = 1;
    span[0,1] = 2;
    span[0,2] = 3;
    span[1,0] = 4;
    span[1,1] = 5;
    span[1,2] = 6;

    for (std::size_t y = 0; y < span.extent(0); y++) // Или `y < h`
    {
        for (std::size_t x = 0; x < span.extent(1); x++) // Или `x < w`
            std::cout << span[y,x] << ' ';
        
        std::cout << '\n';
    }
    // Печатает:
    // 1 2 3
    // 4 5 6
}

Но это пока не поддержали в libstdc++ (стандартной библиотеке компилятора GCC). Зато в libc++ (стандартная библиотека и Clang-а) и MSVC STL работает.

Он поддерживает любое число измерений, не обязательно 2.

А вот так можно передавать его в функцию:

#include <cstddef>
#include <iostream>
#include <mdspan>
#include <vector>

// Здесь 2 - число измерений.
void PrintArray(std::mdspan<int, std::dextents<std::size_t, 2>> span)
{
    for (std::size_t y = 0; y < span.extent(0); y++) // Или `y < h`
    {
        for (std::size_t x = 0; x < span.extent(1); x++) // Или `x < w`
            std::cout << span[y,x] << ' ';
        
        std::cout << '\n';
    }
}

int main()
{
    std::size_t w = 3, h = 2;
    std::vector<int> storage(w * h);
    auto span = std::mdspan(storage.data(), h, w);

    span[0,0] = 1;
    span[0,1] = 2;
    span[0,2] = 3;
    span[1,0] = 4;
    span[1,1] = 5;
    span[1,2] = 6;

    PrintArray(span);
}

Для векторов (в математическом смысле) и матириц, задумывался std::valarray и slice / gslice для пересчета индексов многомерных массивов в одномерные.

В сравнении с vector и array, valarray (в большинстве реализаций) размещается в памяти выровнено, что упрощает применение SSE/AVX инструкций.

пример отсюда: https://en.cppreference.com/w/cpp/numeric/valarray/slice_array

#include <iostream>
#include <valarray>
 
class Matrix
{
    int dim;
    std::valarray<int> data;
public:
    explicit Matrix(int dim, int init = 0)
        : dim{dim}, data(init, dim*dim) {}
    void clear(int value = 0) { data = value; }
    void identity() { clear(); diagonal() = 1; }
    int& operator()(int x, int y) { return data[dim * y + x]; }
 
    std::slice_array<int> diagonal()
    {
        return data[std::slice(0, dim, dim + 1)];
    }
    std::slice_array<int> secondary_diagonal()
    {
        return data[std::slice(dim - 1, dim, dim - 1)];
    }
    std::slice_array<int> row(std::size_t row)
    {
        return data[std::slice(dim * row, dim, 1)];
    }
    std::slice_array<int> column(std::size_t col)
    {
        return data[std::slice(col, dim, dim)];
    }
    template<unsigned, unsigned> friend class MatrixStack;
};
 
template<unsigned dim = 3, unsigned max = 8> class MatrixStack
{
    std::valarray<int> stack;
    unsigned count = 0;
public:
    MatrixStack() : stack(dim * dim * max) {}
    void print_all() const
    { 
        std::valarray<int> row(dim*count);
        for (unsigned r = 0; r != dim; ++r) // screen row
        {
            row = stack[std::gslice(r * dim, {count, dim}, {dim * dim, 1})];
            for (unsigned i = 0; i != row.size(); ++i)
                std::cout << row[i] << ((i + 1) % dim ? " " : " │ ");
            std::cout << '\n';
        }
    }
    void push_back(Matrix const& m)
    {
        if (count < max)
        {
            stack[std::slice(count * dim * dim, dim * dim, 1)]
                = m.data[std::slice(0, dim * dim, 1)];
            ++count;
        }
    }
};
 
int main()
{
    constexpr int dim = 3;
    Matrix m{dim};
    MatrixStack<dim,12> stack;
 
    m.identity();
    stack.push_back(m);
 
    m.clear(1);
    m.secondary_diagonal() = 3;
    stack.push_back(m);
 
    for (int i = 0; i != dim; ++i)
    {
        m.clear();
        m.row(i) = i + 1;
        stack.push_back(m);
    }
 
    for (int i = 0; i != dim; ++i)
    {
        m.clear();
        m.column(i) = i + 1;
        stack.push_back(m);
    }
 
    m.clear();
    m.row(1) = std::valarray<int>{4, 5, 6};
    stack.push_back(m);
 
    m.clear();
    m.column(1) = std::valarray<int>{7, 8, 9};
    stack.push_back(m);
 
    stack.print_all();
}

Перечитал в сети много уверенных утверждений, что одномерный массив с матричным доступом быстрее классического подхода (отдельный массив для каждой строки матрицы). Сделал пробный класс матрицы, накладывающий адресацию на линейный массив и сравнил производительность с классической матрицей (отдельные массивы на каждую строку) при использовании метода Гаусса решения систем линейных уравнений.

Выводы:

1. На MSVC2022 классика значительно (у меня в 3 раза!) быстрее, чем через класс matrix_le. Причина - MSVC2022 не умеет применять SSE/AVX (по ключу /arch:) к таким матрицам (использует там только MMX даже с ключом /arch:AVX2), а к классическим умеет.
2. clang умеет применять SSE/AVX к вычислениям с matrix_le, и на нём оба варианта дают одинаково высокую (как у MS на классике) производительность.
3. Мало тестировать отдельно разные способы доступа к памяти - их надо тестировать в связке с реальным вычислительным алгоритмом, использующим этот доступ.При этом возможны неожиданности вызванные этой связкой.

Дополнительно:

4. Оказалось что любая инкапсуляция матрицы в класс мешает msvc применять SSE / AVX. Если инкапсулировать в (i, j) классическую [i][j] то получается то же самое замедление.
5. Если вычислять [i * nCol + j] непосредственно в методе Гаусса без инкапсуляции, то это тоже мешает msvc применять SSE / AVX, но оптимизация улучшается и производительность получается примерно средней между вариантами с инкапсуляцией матрицы в класс и классическим [i][j].

Дополнение 2:

6. Забавно - продолжил добавлять в программу новые варианты тестов и с какого-то момента clang разучился использовать SSE / AVX в варианте [i][j] - стало что в clang быстро [i * nCol + j], а в msvc быстро [i][j]. При этом я код этих функций не изменял, только добавил в программу новые версии функций с альтернативными вариантами метода Гаусса. Получается умение компилятора использовать SSE / AVX зависит ещё и от "фазы луны".

Код класса матрицы:

//=======================================================================================
// Класс матрицы для системы линейных уравнений 
// Квадратная матрица + столбец свободных членов размером n x (n + 1)
template<class Type> requires(std::is_floating_point_v<Type>)
class matrix_le 
{
public:
    matrix_le(size_t n);
    matrix_le(matrix_le& A);
    ~matrix_le();
    Type& operator()(size_t i, size_t j) {
        return data[i * mCol + j];
    };
    size_t n() { return N; };
    void fill_rand();
    matrix_le& operator=(const matrix_le&);
private:
    Type* data;
    size_t N;
    size_t mCol;    // длина строки (пригодится при экспериментах с выравниванием, поэтому отдельное поле, а не просто n+1)
};

// Методы класса matrix_le
template<class Type> requires(std::is_floating_point_v<Type>)
matrix_le<Type>::matrix_le(size_t n)
{
    N = n;
    mCol = n + 1;
    data = new Type[n * mCol];
};

template<class Type> requires(std::is_floating_point_v<Type>)
matrix_le<Type>::matrix_le(matrix_le& A)
{
    N = A.N;
    mCol = A.mCol;
    data = new Type[N * mCol];
    memcpy(data, A.data, sizeof(Type) * N * mCol); // копирование содержимого матрицы
};

template<class Type> requires(std::is_floating_point_v<Type>)
matrix_le<Type>::~matrix_le()
{
    delete[] data;
};


template<class Type> requires(std::is_floating_point_v<Type>)
void matrix_le<Type>::fill_rand()
{
    for (size_t i = 0; i < N * mCol; i++)
        data[i] = static_cast<Type>(rand()) / static_cast <Type> (RAND_MAX / 10000) - 5000;
};

// Оператор присваивания
template<class Type> requires(std::is_floating_point_v<Type>)
matrix_le<Type>& matrix_le<Type>::operator=(const matrix_le& A)
{
    if (this != &A) // проверка на самоприсваивание
        memcpy(data, A.data, sizeof(Type) * N * mCol); // копирование содержимого матрицы
    
    return *this;
};
// Конец шаблона класса для матрицы линейных уравнений

Доступ через (i, j), потому что MSVC2022 в 2024 году всё ещё не умеет [i, j] из С++23.

Добавил код вычислений

Классическая матрица на указателях:

// Решение системы линейных уравнений оптимизированным методом Гаусса
template<typename Type> requires(std::is_floating_point_v<Type>)
void Gauss_le_classic(Type** A, Type* x, size_t n)
{
    // Начало метода Гаусса (треугольная матрица)
    for (size_t i = 0; i < n; i++)      // столбцы
    {
        for (size_t j = i + 1; j < n; j++)  // строки - вычитание строки i из всех нижележащих
        {
            Type K = A[j][i] / A[i][i];
            for (size_t m = i; m <= n; m++) // столбцы - почленное вычитание от i элемента, включая свободный член
                A[j][m] -= A[i][m] * K;
        };
    };
    // Вычисление x[i] = (b[i] - sum(A[i][j+1])*x[j+1])) / A[i][i]
    for (ptrdiff_t i = n - 1; i >= 0; i--)
    {
        Type sum = 0;
        for (size_t j = i; j < n-1; j++)
            sum += A[i][j+1] * x[j+1];
        x[i] = (A[i][n] - sum) / A[i][i];
    };
};

Вариант с class Matrix

template<typename Type> requires(std::is_floating_point_v<Type>)
void Gauss_le_classic(matrix_le<Type>& A, std::vector<Type> &x)    
{
    size_t n = A.n();
    // Начало метода Гаусса (треугольная матрица)
    for (size_t i = 0; i < n ; i++)     // столбцы (без свободного члена)
    {
        for (size_t j = i + 1; j < n; j++)  // строки - вычитание строки i из всех нижележащих
        {
            Type K = A(j, i) / A(i, i);
            for (size_t m = i; m <= n; m++) // столбцы - почленное вычитание от i элемента, включая свободный член
                A(j, m) -= A(i, m) * K;
        };
    };

    // Вычисление x[i] = (b[i] - sum(A1[i][j+1])*x[j+1])) / A1[i][i]
    for (ptrdiff_t i = n - 1; i >= 0; i--)
    {
        Type sum = 0;
        for (size_t j = i; j < n - 1; j++)
            sum += A(i, j + 1) * x[j + 1];
        x[i] = (A(i, n) - sum) / A(i, i);
    };
};

Тест производительности:

#define fp_type double
size_t n = 1000;
matrix_le<fp_type> A(n);
A.fill_rand();
matrix_le B(A);
std::vector<fp_type> x(n);

perfomance.reset();
SetPriorityClass(hInstance, REALTIME_PRIORITY_CLASS);
for (UINT i = 0; i < n_tests; i++)
{
    B = A;  // восстановить исходную матрицу перед повторным тестом
    perfomance.begin();
        Gauss_le_classic(B, x);
    perfomance.end();
};
SetPriorityClass(hInstance, NORMAL_PRIORITY_CLASS);

infoString << L"n = " << n << L", тестов: " << n_tests << L"\n"
    << L"Гаус на CPU (class matrix_le): " << perfomance.min() << L" мс\n";

// здесь опущен аналогичный замер варианта Type** A

MessageBox(0, infoString.str().c_str(), L"Результат", MB_OK);

класс perfomance - мой велосипед скрывающий QueryPerformanceCounter и запоминающий минимальное значение времени из серии тестов.

Но и даже без замеров можно поставить брейки на функции и в дизассемблированном коде увидеть что MSVC в классическом (Type** A) варианте смог использовать AVX2 (компиляция с ключом /arch:AVX2), а в варианте с class Matrix нет, хотя всё компилировалось совместно.

PS: функционал сокращен чтобы не загромождать пример - нет проверки и корректировки нулей на главной диагонали и проверки что размер вектора x достаточен для размещения в нём результата. В конструкторе копирования и операторе присваивания матрицы нет проверки и обработки размеров матрицы - предполагается что копируются/присваиваются только матрицы одинакового размера. И даже так пример получился громоздким.

Тестовая матрица однократно заполняется случайными числами в диапазоне -5000 -:- +5000. Перед повторными замерами матрица восстанавливается из резервной копии чтобы данные были одинаковыми. Это позволяет отделить время вычислений от времени копирования данных.