Предупреждение

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

Эта страница описывает с помощью примера, как реализовать новый лёгкий тип выражения в Eigen. Это состоит из трёх частей: сам тип выражения, класс свойств, содержащий информацию о времени компиляции для выражения, и класс вычислителя, который используется для вычисления выражения до матрицы.

НЕОБХОДИМО СДЕЛАТЬ: Напишите страницу, описывающую дизайн, со сведениями о векторизации и т. д., и сделайте ссылку на неё здесь.

Настройка

Циркулянтная матрица — это матрица, где каждый столбец такой же, как столбец слева, за исключением того, что он циклически сдвинут вниз. Например, вот 4×4 циркулянтная матрица:

\[ \begin{bmatrix} 1 & 8 & 4 & 2 \\ 2 & 1 & 8 & 4 \\ 4 & 2 & 1 & 8 \\ 8 & 4 & 2 & 1 \end{bmatrix} \]

Циркулянтная матрица однозначно определяется своим первым столбцом. Мы хотим написать функцию makeCirculant которая, получив первый столбец, возвращает выражение, представляющее циркулянтную матрицу.

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

Начало работы

Мы представим файл, реализующий функцию makeCirculant, по частям. Мы начнём с включения соответствующих заголовочных файлов и объявления класса выражения, который мы назовём Circulant. Функция makeCirculant будет возвращать объект этого типа. Класс Circulant фактически является шаблонным классом; шаблонный аргумент ArgType относится к типу вектора, переданного функции makeCirculant.

#include <Eigen/Core>
#include <iostream>
 
template <class ArgType> class Circulant;

Класс свойств

Для каждого класса выражения X, должен быть класс свойств Traits<X> в пространстве имён Eigen::internal содержащий информацию о X известную как информация времени компиляции.

Как объяснено в Настройка, мы спроектировали класс выражения Circulant для ссылок на плоские матрицы. Элементы циркулянтной матрицы имеют тот же тип, что и элементы вектора, переданного функции makeCirculant. Тип, используемый для индексации элементов, также такой же. Опять же, для простоты, мы будем возвращать только матрицы в порядке следования столбцов. Наконец, циркулянтная матрица является квадратной матрицей (количество строк равно количеству столбцов), и количество строк равно количеству строк столбцового вектора, переданного функции makeCirculant . Если это вектор динамического размера, размер циркулянтной матрицы неизвестен во время компиляции.

Это приводит к следующему коду:

namespace Eigen {
  namespace internal {
    template <class ArgType>
    struct traits<Circulant<ArgType> >
    {
      typedef Eigen::Dense StorageKind;
      typedef Eigen::MatrixXpr XprKind;
      typedef typename ArgType::StorageIndex StorageIndex;
      typedef typename ArgType::Scalar Scalar;
      enum { 
        Flags = Eigen::ColMajor,
        RowsAtCompileTime = ArgType::RowsAtCompileTime,
        ColsAtCompileTime = ArgType::RowsAtCompileTime,
        MaxRowsAtCompileTime = ArgType::MaxRowsAtCompileTime,
        MaxColsAtCompileTime = ArgType::MaxRowsAtCompileTime
      };
    };
  }
}

Класс выражения

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

НЕОБХОДИМО СДЕЛАТЬ: А как насчёт Nested typedef? Похоже, он необходим; это только временная мера?

template <class ArgType>
class Circulant : public Eigen::MatrixBase<Circulant<ArgType> >
{
public:
  Circulant(const ArgType& arg)
    : m_arg(arg)
  { 
    EIGEN_STATIC_ASSERT(ArgType::ColsAtCompileTime == 1,
                        YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
  }
 
  typedef typename Eigen::internal::ref_selector<Circulant>::type Nested; 
 
  typedef Eigen::Index Index;
  Index rows() const { return m_arg.rows(); }
  Index cols() const { return m_arg.rows(); }
 
  typedef typename Eigen::internal::ref_selector<ArgType>::type ArgTypeNested;
  ArgTypeNested m_arg;
};

Вычислитель

Последний большой фрагмент реализует вычислитель для выражения Circulant. Вычислитель вычисляет элементы циркулянтной матрицы; это делается в функции-члене .coeff(). Элементы вычисляются путём поиска соответствующего элемента вектора, из которого построена циркулянтная матрица. Получение этого элемента может быть на самом деле нетривиальным, когда циркулянтная матрица построена из вектора, заданного сложным выражением, поэтому мы используем вычислитель, соответствующий вектору.

Константа CoeffReadCost записывает стоимость вычисления элемента циркулянтной матрицы; мы игнорируем вычисление индекса и говорим, что это то же самое, что и стоимость вычисления элемента вектора, из которого построена циркулянтная матрица.

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

namespace Eigen {
  namespace internal {
    template<typename ArgType>
    struct evaluator<Circulant<ArgType> >
      : evaluator_base<Circulant<ArgType> >
    {
      typedef Circulant<ArgType> XprType;
      typedef typename nested_eval<ArgType, XprType::ColsAtCompileTime>::type ArgTypeNested;
      typedef typename remove_all<ArgTypeNested>::type ArgTypeNestedCleaned;
      typedef typename XprType::CoeffReturnType CoeffReturnType;
 
      enum { 
        CoeffReadCost = evaluator<ArgTypeNestedCleaned>::CoeffReadCost,
        Flags = Eigen::ColMajor 
      };
      
      evaluator(const XprType& xpr)
        : m_argImpl(xpr.m_arg), m_rows(xpr.rows())
      { }
 
      CoeffReturnType coeff(Index row, Index col) const
      {
        Index index = row - col;
        if (index < 0) index += m_rows;
        return m_argImpl.coeff(index);
      }
 
      evaluator<ArgTypeNestedCleaned> m_argImpl;
      const Index m_rows;
    };
  }
}

Точка входа

После всего этого функция makeCirculant очень проста. Она просто создаёт объект выражения и возвращает его.

template <class ArgType>
Circulant<ArgType> makeCirculant(const Eigen::MatrixBase<ArgType>& arg)
{
  return Circulant<ArgType>(arg.derived());
}

Простая функция main для тестирования

Наконец, короткая функция main, которая показывает, как можно вызвать функцию makeCirculant.

int main()
{
  Eigen::VectorXd vec(4);
  vec << 1, 2, 4, 8;
  Eigen::MatrixXd mat;
  mat = makeCirculant(vec);
  std::cout << mat << std::endl;
}

Если все фрагменты объединить, будет получен следующий результат, показывающий, что программа работает как ожидается:

1 8 4 2
2 1 8 4
4 2 1 8
8 4 2 1