Общее оглавление

Цели и критерии

Цели: изучение работы с двумерными динамическими массивами.

Критерии полноты

1. Должна быть реализована отдельная функция, выполняющая действие, указанное в задании.
2. Данная функция должна принимать произвольную (не обязательно квадратную) матрицу, не копируя её. Элементы матрицы — числа с плавающей запятой.
3. В виде отдельной функции должны быть реализованы по два автоматических теста на “положительный” и “отрицательный” ожидаемые ответы функции. Матрицы в тестах должны быть разных размеров.

Примеры

Данная работа предполагает знание материала по многомерным массивам.

amatrix.hpp

Далее применяется способ 2 с прямым использованием упакованного двумерного массива, хранящегося в памяти как массив строк. Например, для матрицы размеров 3×4, используя нумерацию элементов, начинающуюся с нуля по каждому измерению, имеем (каждая строка матрицы выделена своим цветом)

Ряд вспомогательных определений, которые можно использовать при решении данной самостоятельной работы, собран в файл amatrix.hpp (при копировании файла префикс 0700- надо убрать). Итак, пусть информация о матрице (адрес массива, количество строк и столбцов) хранится в объектах структуры Matrix_info:

/// Описание матрицы.
struct Matrix_info
{
  double *data; ///< Указатель на массив элементов.
  size_t rows,  ///< Количество строк.
         cols;  ///< Размер строки (количество столбцов).
};

Все функции, оперирующие матрицами, принимают указатели на соответствующие объекты структуры Matrix_info. Если функция не изменяет объект матрицы, то соответствующий указатель объявляется как указатель на константу (с помощью ключевого слова const). Например, проверка двух матриц на равенство принимает два таких указателя:

/// Сравнить матрицы на равенство.
inline bool matrices_are_equal(const Matrix_info *a, const Matrix_info *b)
{
  // Если размеры матриц не равны, то и матрицы не равны.
  if (a->rows != b->rows || a->cols != b->cols)
    return false;

  // Размеры совпали, сравним поэлементно.
  const auto elements = matrix_size(a);
  auto ad = a->data, bd = b->data;
  for (size_t i = 0; i < elements; ++i)
    if (ad[i] != bd[i])
      return false;

  // Равны поэлементно => равны.
  return true;
}

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

Определив переменную типа Matrix_info, следует инициализировать её. Пустая матрица без выделенного куска памяти:

Matrix_info m = {}; // инициализировать нулями все поля структуры

Создать матрицу заданного размера, заполненную неизвестными значениями, предполагая, что память ещё не выделялась:

alloc_matrix(&m, 10, 20); // матрица 10x20

Если память ранее выделялась, то можно “пересоздать” матрицу, указав новые размеры:

realloc_matrix(&m, 20, 10); // переформатировать в матрицу 20x10

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

free_matrix(&m);

C++, HTML

Указатели или ссылки

Подробнее об указателях и ссылках см. здесь.

В примере amatrix.hpp можно было все функции сделать принимающими ссылки на Matrix_info, а не указатели. Это было бы как раз “в духе” C++. Однако пока было принято решение показать на примере стиль C (хотя и с C++-кодом).

Модификатор const при применении к указателю или структуре с полем-указателем применяется к самому значению указателя (адресу), а не к значению, на которое он указывает. Поэтому, хотя ряд функций в amatrix.hpp и принимает const Matrix_info, сигнализируя, что они не будут изменять матрицу, формально они имеют право изменять значения элементов. В случае ссылок, const применяется к значению, на которое ссылается ссылка (что и предполагается обычно программистом).

Пример 1

Задача. Для заданной матрицы (a_ij) построить транспонированную матрицу (a_ji).

Решение

Итак, “простой” алгоритм “по определению” выглядит как двойной цикл вида:

1. Выделить память на транспонированную матрицу.
2. Для каждой строки i:
   • для каждого столбца j:
     • элемент скопировать в (j, i) транспонированной матрицы элемент (i, j) исходной матрицы.

Используя amatrix.hpp, это можно записать как

void make_transpose_naive(const Matrix_info *m, Matrix_info *mt)
{
  // размеры
  const auto rows = m->rows, cols = m->cols;
  
  // подготовить место
  realloc_matrix(mt, cols, rows);
  
  // скопировать
  for (size_t i = 0; i < rows; ++i)
    for (size_t j = 0; j < cols; ++j)
      set_element(mt, j, i, get_element(m, i, j));
}

Можно избавиться от умножений, зашитых в функции set_element и get_element, так как и по исходной и по транспонированной матрице мы проходим с постоянным шагом.

void make_transpose_naive(const Matrix_info *m, Matrix_info *mt)
{
  // размеры
  const auto rows = m->rows, cols = m->cols;
  
  // подготовить место
  realloc_matrix(mt, cols, rows);
  
  // код без умножений
  auto m_row = m->data; // начало текущей строки в m
  auto mt_col = mt->data; // начало текущего столбца в mt
  for (size_t i = 0; i < rows; ++i)
  {
    auto mt_pos = mt_col; // позиция записи
    for (size_t j = 0; j < cols; ++j)
    {
      // элементы строки идут в массиве друг за другом непосредственно (шаг = 1)
      *mt_pos = m_row[j];
      // элементы столбца идут в массиве с шагом в одну строку
      // (шаг = rows -- количество столбцов в mt)
      mt_pos += rows;
    }
    
    ++mt_col; // перейти на следующий столбец в mt
    m_row += cols; // перейти к следующей строке в m
  }
}

C++, HTML

Пример 2

Задача. Определить, являются ли строки или столбцы за пределами наибольшей по размеру квадратной угловой подматрицы нулевыми. В случае, если матрица квадратная, возвращать истину. Данному критерию удовлетворяют n×m-матрицы следующего вида (помимо случая n = m):

Решение

Определить случай в соответствии с соотношением размеров матрицы. Проверить на равенство нулю все дополнительные строки (в случае 1), либо все элементы в позициях (n+1):m каждой строки (“хвостах” строк; в случае 2). При использовании нумерации, начинающейся с нуля, а не с единицы, это позиции n:(m–1).

Итоговый код:

bool are_outer_lines_zeroes(const Matrix_info *mtx)
{
  const auto data = mtx->data;
  const auto rows = mtx->rows, cols = mtx->cols;

  if (rows < cols)
  {
    // Строк меньше, чем столбцов -- проверим "хвосты" строк на равенство нулю.
    // Длина "хвоста".
    const auto diff = cols - rows;
    for (size_t i = 0; i < rows; ++i)
    {
      // Вычислим указатель на первый элемент "хвоста" i-й строки.
      const auto tail = data + i * cols + rows;
      // Проверим элементы "хвоста" на равенство нулю.
      for (size_t j = 0; j < diff; ++j)
        if (tail[j] != 0)
          return false;
    }
  }
  else if (cols < rows)
  {
    // Столбцов меньше, чем строк -- проверим "лишние" строки на равенство нулю.
    // Вычислим указатель на первую "лишнюю" строку,
    // все последующие элементы массива принадлежат "лишним" строкам, проверим их на равенство нулю.
    const auto excess = data + cols * cols;
    const auto size  = (rows - cols) * cols; // сколько оставшихся элементов
    for (size_t i = 0; i < size; ++i)
      if (excess[i] != 0)
        return false;
  }
  
  // Либо матрица является квадратной, либо проверка в ветках if прошла успешно.
  return true;
}

C++, HTML

Пример 3: умножение матриц

Пусть даны n×s-матрица A = (a_ij) и s×m-матрица B = (b_jk). Тогда произведением этих матриц AB называется n×m-матрица C = (c_ik) такая, что

Имеем тройной цикл: указанную сумму (скалярное произведение строки A на столбец B) надо посчитать для всех i и k, всего это (n s m) умножений.

Если реализовать это определение “в лоб” на основе функционала из amatrix.hpp, то получим следующий код (возвращает ложь, если матрицы нельзя перемножить, истину, если операция умножения была выполнена):

bool multiply_v1(const Matrix_info *a, const Matrix_info *b, Matrix_info *c)
{
  const auto n = a->rows, s = a->cols, m = b->cols;
  if (n == 0 || s == 0 || m == 0 || s != b->rows)
    return false;

  realloc_matrix(c, n, m);
  for (size_t i = 0; i < n; ++i)
  {
    for (size_t k = 0; k < m; ++k)
    {
      double sum = 0.;
      for (size_t j = 0; j < s; ++j)
        sum += get_element(a, i, j) * get_element(b, j, k);
      set_element(c, i, k, sum);
    }
  }

  return true;
}

Насколько он хорош или плох? Как и в примере с транспонированием можно заметить, что функции get_element и set_element выполняют лишнюю работу, так как в действительности не требуется умножать каждый раз старший индекс на размер, ведь во всех циклах шаг по матрицам постоянный. Убрав их, получим вариант 2:

bool multiply_v2(const Matrix_info *a, const Matrix_info *b, Matrix_info *c)
{
  const auto n = a->rows, s = a->cols, m = b->cols;
  const auto b_data = b->data; // Указатель на массив элементов матрицы b.
  if (n == 0 || s == 0 || m == 0 || s != b->rows)
    return false;

  realloc_matrix(c, n, m);
  auto a_row = a->data,  // Указатель на текущую строку матрицы a.
       c_item = c->data; // Указатель на текущий элемент матрицы c.

  for (size_t i = 0; i < n; ++i, a_row += s)
  {
    for (size_t k = 0; k < m; ++k)
    {
      double sum = 0.;
      auto b_item = b_data + k; // k-й столбец в j-й строке матрицы b.
      for (size_t j = 0; j < s; ++j, b_item += m)
        sum += a_row[j] * (*b_item);
      *c_item++ = sum;
    }
  }

  return true;
}

Однако, скорее всего, мы не увидим заметного увеличения быстродействия, так как на современных процессорах в случае больших матриц даже вариант 1 будет упираться не в скорость выполнения арифметических операций, а в скорость работы подсистемы памяти. Главный недостаток что варианта 1, что варианта 2 в том, что внутренний цикл пробегает по памяти с большим шагом m*sizeof(элемента матрицы) байт. Поэтому в хорошей реализации умножения матриц “по определению” внутренний цикл будет двигаться подряд как по первой матрице, так и по второй матрице. Было бы очень просто “исправить” код, если бы мы вычисляли не AB, а AB^T. Впрочем, здесь мы тоже можем поменять циклы, немного пожертвовав эффективностью: вся текущая строка матрицы C будет набором сумм, обновляемых за один проход внутреннего цикла. Естественно, теперь матрицу C необходимо предварительно обнулить.

bool multiply_v3(const Matrix_info *a, const Matrix_info *b, Matrix_info *c)
{
  const auto n = a->rows, s = a->cols, m = b->cols;
  if (n == 0 || s == 0 || m == 0 || s != b->rows)
    return false;

  realloc_matrix(c, n, m);
  fill_matrix(c, 0.);
  auto a_row = a->data, // Указатель на текущую строку матрицы a.
       c_row = c->data; // Указатель на текущую строку матрицы c.

  for (size_t i = 0; i < n; ++i, a_row += s, c_row += m)
  {
    auto b_row = b->data; // Указатель на текущую строку матрицы b.
    for (size_t j = 0; j < s; ++j, b_row += m)
    {
      const auto aij = a_row[j];
      for (size_t k = 0; k < m; ++k)
        c_row[k] += aij * b_row[k];
    }
  }

  return true;
}

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

Итак, результаты оценки времени работы (в секундах на одно умножение матриц) представленных выше трёх вариантов на системе с процессором AMD K10 3.6ГГц и 16Гб DDR3-1600 под управлением ОС Windows 7 x64 SP1 для двух компиляторов: Microsoft Visual C++ 2015 x64 (MSVC; стандартная конфигурация Release для платформы x64 с добавлением опции /Qpar) и TDM-GCC 5.1 (g++; опции компиляции: -march=amdfam10 -fomit-frame-pointer -fexpensive-optimizations -O3). Оба этих компилятора поддерживают векторизацию циклов (использовались команды набора SSE2). Оценка погрешности измерения времени получена путём трёхкратного повторения теста.

Последняя строчка таблицы показывает отношение скорости (в умножениях чисел в секунду), достигнутых вариантом 3 (отдельно для каждого компилятора) к средним геометрическим скоростей каждого варианта (отдельно для каждого компилятора). Хорошо видно, что варианты 1 и 2 отличаются незначительно, а вариант 3 быстрее их 7–10 раз.

Вариант “asm” — это вариант-3, вручную переведённый автором на ассемблер x86-64 с минимальной оптимизацией (векторизация и конвейеризация вложенного цикла — обрабатывает по четыре числа за итерацию). Данный вариант демонстрирует наилучшие результаты на матрицах небольшого размера.

Использованный для тестирования процессор снабжён кэшами трёх уровней: 64 кб L1, 512 кб L2 (эксклюзивный), 6 Мб L3. Рассмотрим влияние размера на относительную скорость на перемножении квадратных матриц такого размера, что все три матрицы (операнды и результат) помещаются в кэше соответствующего уровня. Колонки “mulps” (multiplications per second) содержат скорость в “умножениях в секунду” (чисел в плавающей точке) для каждого варианта, при этом одно ядро процессора имеет теоретическую производительность 14.4 млрд операций в плавающей точке в секунду и теоретический максимум mulps равен 4800 млн.

Интересно, что скорость “быстрых” вариантов растёт при увеличении размеров матриц до тех пор, пока преобладает влияние кэшей L1 и L2 (суммарно 576 кб). При выходе за пределы L3 скорость падает почти вдвое, что соответствует соотношению пропускных способностей внутренних шин. Впрочем, до теоретического максимума производительности везде далеко: даже в лучшем случае от него достигается лишь 40%.

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

C++, HTML

Варианты

1

Определить, является ли произвольная матрица диагональной.

2

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

3

Определить, является ли произвольная матрица верхнетреугольной.

4

Определить, является ли произвольная матрица кососимметричной.

5

Определить, является ли произвольная матрица нижнетреугольной.

6

Определить, является ли произвольная матрица трёхдиагональной.

7

Определить, обладает ли произвольная матрица свойством диагонального преобладания.

8

Строки (или столбцы) ортогональной матрицы составляют ортонормированный базис (являются попарно перпендикулярными друг другу векторами единичной длины). Транспонирование ортогональной матрицы даёт обратную матрицу.

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

9

Определить, сколько нулевых столбцов содержит матрица. Функция должна возвращать целое число (тип size_t).

Замечание. Данный вариант отличается от прочих тем, что требует не проверки некоего условия, а вычисления значения. Реализовать три теста: с матрицей без нулевых столбцов и с двумя матрицами с несколькими нулевыми столбцами (разное количество нулевых столбцов).

10

Определить, является ли произвольная матрица антидиагональной.

11

Пример

1 2 3 4
8 0 5 3
1 6 0 2
7 1 8 1

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

12

Определить, является ли произвольная матрица L-матрицей.

13

Пример

5 3 7 1
1 5 3 7
7 1 5 3
3 7 1 5

Определить, является ли произвольная матрица циркулянтом.

14

Пример

1 2 3 4
2 3 4 1
3 4 1 2
4 1 2 3

Определить, является ли произвольная матрица антициркулянтом.

15

Пример

1 5 0 1 1
0 0 1 2 2
0 0 0 1 1
0 0 0 0 0
0 0 0 0 0

Определить, является ли произвольная матрица ступенчатой.

16

Пример

1 2 3 4 5 6 7
8 1 2 3 4 5 6
9 8 1 2 3 4 5

Определить, является ли произвольная матрица матрицей Тёплица.

17

Пример

1 2 3 4 5
2 3 4 5 6
3 4 5 6 7
4 5 6 7 8
5 6 7 8 9

Определить, является ли произвольная матрица матрицей Ганкеля.

Общее оглавление

Кувшинов Д.Р. © 2015