§ 2. Библиотеки для численных методов

Многомерные массивы и ленивые вычисления

Библиотеки многомерных массивов выполняют две функции:

• автоматическая векторизация кода — генерация кода, который использует векторные регистры процессора для поэлементных операций над многомерными массивами,
• и использование ленивых вычислений для генераций самих поэлементных операций без нагромождения вложенных циклов.

Ленивые вычисления — это прием в программировании, который позволяет откладывать вычисления до тех пор, пока их результат не понадобится использовать в обычных вычислениях. Этот прием широко используется в функциональных языках для следующих целей:

• использование кода как данных для передачи по сети, записи в файл или вывода на экран;
• выполнение кода при необходимости (если результат не понадобился, то и код не будет зря выполняться);
• оптимизация математических и аналогичных выражений, используя набор правил.

В C++ ленивые вычисления применяются, в основном, для оптимизации математических выражений, работающих с многомерными массивами.

Ленивые вычисления в C++ реализуются путем создания классов, которые описывают синтаксическое дерево математических операций, которые необходимо оптмизировать (арифметические операции, тригонометрические, гипергеометрические и специальные функции и другие). Для каждой операции создается отдельный класс, содержащий аргументы операции в виде полей, а сама операция реализуется в методе evaluate для отдельного элемента массива.

template <class E1, class E2>
class Plus: public Expression {
public:
    using value_type =
        typename std::common_type<typename E1::value_type, typename E2::value_type>::type;
private:
    const E1& _a;
    const E2& _b;
public:
    explicit Plus(const E1& a, const E2& b): _a(a), _b(b) {}
    value_type evaluate(int i) { return this->_a.evaluate(i) + this->_b.evaluate(i); }
    int size() const { return this->_a.size(); }
    void display(std::ostream& out) const {
        out << "Plus(" << this->_a << ", " << this->_b << ')';
    }
};

В классе присутствует метод для определения размера массива, вывода операции на экран в виде выражения, а также возвращаемый тип данной операции value_type. В примере показана операция сложения для элементов одномерного массива, но подход легко обобщить для нескольких измерений.

В дополнение к классам, описывающим операции, идут глобальные функции, которые вызывают операции для каждого элемента массива.

template <class E>
Array<typename E::value_type> evaluate(const E& expr) {
    using value_type = typename E::value_type;
    const auto n = expr.size();
    Array<value_type> result(n);
    for (int i=0; i<n; ++i) {
        result[i] = expr.evaluate(i);
    }
    return result;
}

Здесь Array<T> является одномерным массивом. В этой функции сначала определяется размер массива, а затем вызывается метод evaluate для каждого элемента массива. Обычно, в этой функции делают различные оптимизации, позволяющие векторизовать поэлементные операции, т.е. использовать векторные регистры процессора, чтобы применять одну и ту же операцию одновременно к нескольким элементам массива. Это значительно ускоряет код и является одной из причин создания библиотек для эффективной работы с многомерными массивами.

Помимо функции evaluate полезно добавить оператор вывода для отладки выражений.

template <class E>
typename std::enable_if<std::is_base_of<Expression,E>::value,std::ostream&>::type
operator<<(std::ostream& out, const E& expr) {
    expr.display(out);
    return out;
}

Поскольку выражение может иметь любой тип, то мы используем специальный шаблон std::enable_if, который позволяет определить оператор вывода только для тех типов, для которых класс Expression является базовым. Вторым аргументом этого шаблона является возвращаемый тип. Если первый аргумент равен true, тогда тип определяется в виде поля type шаблона std::enable_if, иначе никакого поля в нем нет. А если поля нет, то значит оператор вывода для соответствующего типа не определен. Сам базовый тип Expression является пустым, такие типы иногда называют классами-метками, поскольку они используются только во время компиляции для метапрограммирования.

Наконец, надо определить оператор сложения для наших выражений, чтобы использовать знакомый и интуитивный математический синтаксис.

template <class E1, class E2>
typename std::enable_if<std::is_base_of<Expression,E1>::value &&
                        std::is_base_of<Expression,E2>::value,
                        Plus<E1,E2>>::type
operator+(const E1& a, const E2& b) {
    return Plus<E1,E2>(a,b);
}

Здесь тоже используется std::enable_if, но с более сложным условием. Сам оператор возвращает выражение, а не результат вычислений, т.е. является функцией-конструктором для выражения Plus.

Результатом наших усилий является набор классов и функций, которые позволяют генерировать эффективные и легко векторизуемые поэлементные циклы, используя знакомую математическую нотацию. Например, код

Array<T> a{1,2,3}, b{4,5,6}, c{7,8,9};
Array<T> d = a+b+c;

сгенерирует цикл следующего вида.

for (int i=0; i<n; ++i) {
    d[i] = a[i] + b[i] + c[i];
}

Данный подход работает для поэлементных операций, но аналогично можно реализовать любые другие операции над массивами. Поэлементные реализуют в первую очередь из-за оптимизаций, которые можно таким образом получить. В заданиях мы добавим в этот код больше выражений.

Наиболее популярными на данный момент библиотеками многомерных массивов являются Eigen и Armadillo. Первой библиотекой, которая использовала описанные в данном разделе приемы программирования, является Blitz.

Матрицы и векторы

Кроме многомерных массивов в программах используются их частные случаи: матрицы и векторы. Матрицы и векторы в программе — это двух и трехмерные массивы, для которых определены арифметические операции, а также функции, реализующие методы решения систем линейных уравнений, методы нахождения определителя матрицы и все остальные. В программе различают матрицы и векторы динамического размера (как правило, это большие матрицы, описывающие системы линейных уравнений) и статического размера (как правило, это матрицы, описывающие геометрические преобразования). Для матриц и векторов статического размера используются отдельные типы данных, поскольку геометрические операции можно оптимизировать с помощью векторных регистров, в которые маленький вектор помещается целиком, и специальных инструкций процессора, которые работают с такими регистрами. Таким образом, мы имеем большие матрицы и векторы, для которых определены операции и методы линейной алгебры, а также маленькие матрицы и векторы, для которых определены геометрические операции и методы.

Наиболее популярными на данный момент библиотеками для матриц и векторов являются Eigen и GLM. Blitz содержит классы для матриц и векторов статических размеров, но не содержит отдельные классы для динамического размера, а также большого количества операций для этих типов данных.

Наиболее популярными библиотеками матричных и векторных операций являются BLAS и LAPACK. BLAS содержит операции линейной алгебры, а LAPACK более высокоуровневые методы раложения матриц и решения СЛАУ. Оригинальные версии этих библитек написаны на Фортране, однако существуют аналоги этих библиотек (с тем же самым программным интерфейсом), оптимизированные под различные архитектуры процессоров, а также для графических ускорителей. Среди них можно выделить открытый OpenBLAS (бывший GotoBLAS), проприетарный MKL (Math Kernel Library) и опять же Eigen.

Библиотеки с остальными методами перечислены ниже.

Структура программы

Несмотря на то что существует множество библиотек с уже готовыми программными реализациями численных методов, они не всегда подходят для решения реальных задач. Также для реальных задач часто нужен уникальный численный метод, которого просто нет в библиотеках. В таком случае его реализуют в основной программе.

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

Класс является наиболее удачной абстракцией не только из-за большого количества параметров и входных данных, которые неудобно передавать в функцию в виде аргументов, но и из-за того что у класса-солвера может быть несколько реализаций: одна для многопроцессорной машины, другая для видеокарты, третья для кластера. В случае класса несколько реализаций можно сделать путем объявления основных методов виртуальными и создания отдельных дочерних классов, каждый из которых реализует данные методы уникальным способом. Затем в основной программе вместо объекта базового класса используется указатель на базовый класс, которому присваивается адрес объекта одного из дочерних классов (по выбору пользователя). Такая структура позволяет эффективно разрабатывать несколько реализаций одного и того же метода под различные архитектуры вычислительных машин без нагромождения условных переходов и макросов, которые включают или выключают реализации во время компиляции.

Верификация и тестрование программ

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

Как же писать модульные тесты для численных методов? Наиболее подходящим вариантом на данный момент являются так называемые библиотеки тестирования на основе свойств. В этих библиотеках тесты пишутся таким образом, что для параметров и аргументов численного метода перебираются всевозможные значения и их комбинации, и проверяется, работает ли метод при этих значениях. Тесты, написанные таким образом наиболее востребованы для кода, который реализует арифметические или аналогичные операции, для которых сущестуют математические свойства. Например, мы знаем, что если матрицу домножить на обратную ей, то получится единичная матрица. Это свойство выполняется для всех матриц, для которых существует обратная. Это свойство можно проверить с помощью такой библиотеки.

Как же работают эти библиотеки? В них параметры и аргументы метода генерируются псевдослучайным образом в виде целых или вещественных чисел. Затем эти числа программист преобразует к фактическим аргументам метода, вызывает этот метод и выполняет проверки свойств. Если хотя бы одно свойство не выполняется, то библиотека выдает комбинацию аргументов, при которой произошла ошибка. Перед выдачей аргументов библиотека пытается найти наименьшие их значения, для которых свойства не выполняются. Это позволяет избавиться от трудно читаемых больших чисел в выводе, а также в некоторых случаях сразу понять, какой аргумент привел к такой ошибке (если все остальные аргументы равны нулю). Этот процесс очень похож на поиск минимума функции, только функция у нас возвращает булевы значения, а для уточнения точки минимума все аргументы уменьшаются до минимальлно возможных значений, приводящих к ошибке.

Следующий код проверяет корректность реализации кодировки BASE64.

void test_base64_symmetry_encode_decode() {
    using namespace sys::test;
    auto prng = current_test->prng();
    falsify(
        [&prng] (const Argument_array<1>& params) {
            auto size = params[0];
            std::string input;
            input.reserve(size);
            std::uniform_int_distribution<char> dist_char(
                std::numeric_limits<char>::min(), std::numeric_limits<char>::max());
            for (size_t i=0; i<size; ++i) { input += dist_char(prng); }
            std::string encoded(base64_encoded_size(input.size()), '_');
            sys::base64_encode(input.data(), input.size(), &encoded[0]);
            std::string decoded(base64_max_decoded_size(encoded.size()), '_');
            auto actual_size = sys::base64_decode(encoded.data(), encoded.size(), &decoded[0]);
            decoded.resize(actual_size);
            return expect(value(input) == value(decoded));
        },
        make_parameter<size_t>(0,4097));
}

Задания

Перед выполнением заданий скопируйте исходный код программы, используя следующие команды.

git clone https://courses.igankevich.com/numerical-methods/linalg.git
cd linalg
meson build
cd build
ninja
./src/linalg/linalg

Перед вами программа, которая реализует ленивые вычисления на C++. Вычисления производятся поэлементно над одномерными массивами.

Операторы сравнения1 балл

• Реализуйте поэлементный бинарный оператор сравнения <. Оператор сравнивает соответствующие элементы переданных в него выражений.
• Реализуйте оператор all, который возвращает истину, если все элементы выражения являются истиной. Этот оператор непоэлементный, поэтому для него придется написать отдельную глобальную функцию evaluate. Реализация должна вычислять минимально возможное число элементов.
• Реализуйте оператор any, который возвращает истину, если хотя бы один элемент выражения является истиной. (Реализация аналогична all.)
• Реализуйте тернарный оператор ?:. Поскольку переопределить в C++ его нельзя, то роль функции-конструктора будет играть функция where. Оператор where(a,b,c) должен вычислять выражение b или c в зависимости от предиката a, который также является выражением.

Несколько измерений2 балла

• Реализуйте многомерные массивы и выражения. Количество измерений известно во время компиляции. Для этого добавьте параметр шаблона N к классу Array и измените индексацию массивов и выражений, используя формулы для row-major из статьи.
• Реализуйте операцию транспонирования для многомерного массива. Эта операция переставляет индексы в указанном в качестве аргументов порядке. Операция определена для массивов с количеством измерений больше одного.

Определитель3 балла

Реализуйте вычисление определителя для квадратной матрицы произвольного размера. Для этого воспользуйтесь формулой для символа Леви—Чивиты и формулой для детерминанта. Проверьте себя либо с помощью детерминантов для небольших матриц, либо с помощью Wolfram Alpha.

Видео

2023

2022