Титульный лист

Разделим матрицы A, B и C на равные по размеру блочные матрицы:
, , , где .
Тогда,

(1)


Однако с помощью этой процедуры не удастся уменьшить количество умножений. Как и в обычном методе, потребуется 8 умножений.
Теперь определим новые элементы:






.
которые затем используются для выражения . Таким образом, понадобиться всего 7 умножений на каждом этапе рекурсии. Элементы матрицы C выражаются из по формулам:




Рекурсивный процесс продолжается n раз, до тех пор пока размер матриц не станет достаточно малым, далее используют обычный метод умножения матриц. Это делают из-за того, что алгоритм Штрассена теряет эффективность по сравнению с обычным на малых матрицах в силу большего числа сложений. Оптимальный размер матриц для перехода к обычному умножению зависит от характеристик процессора и языка программирования и на практике лежит в пределах от 32 до 128.
Для реализации в курсовой работе выбран модифицированный метод Штрассена позволяющий выполнять перемножения матриц не дополняя их до квадратных.
Пусть даны две матрицы размерностями NxM и MxK. Тогда берется максимальная степень числа 2 меньшая, чем N,M и K. Обозначим ее P. Тогда первая матрица делится на 4 размерностями PxP, Px(M-P), (N-P)xP, (M-P)x(M-P). Аналогично, вторая матрица разделяется на 4. Далее вычисления выполняются по формулам (1).
2. Реализация программы

Для описания матрицы в программе описана структура:

typedef struct Matrix
{
int rows; //количество строк
int columns; //количество столбцов
int realrows; //Реально выделенных строк - для аккуратной работы с диспетчером памяти
double ** data; //данные
} Matrix;
Рассмотрим основные функции программы.
double get(Matrix * M, int row, int col)
Функция get используется для доступа к отдельному элементу матрицы. В качестве параметров передаются указатель на структуру описывающую матрицу и номера строки и столбца элемента. Функция возвращает значение элемента.
void set(Matrix * M, int row, int col, double value)
Функция set используется для записи значения в матрицу. В качестве параметров передаются указатель на структуру описывающую матрицу, номера строки и столбца элемента и значение для записи. Результатом работы функции является изменение матрицы переданной в качестве параметра.
void RegRAM(int delta)
Функция RegRAM является вспомогательной. Используется для учета общего количества задействованной для вычислений памяти. Необходима для кооретной обработки программой ситуации нехватки оперативной памяти.
Matrix * Create(int rows,int columns)
Функция Create используется для создания новой матрицы указанного размера. В качестве параметров передается количество строк и столбцов матрицы. Функция возвращает указатель на созданную структуру. Если создание матрицы не возможно, то возвращается NULL.
Matrix * Generate(int rows, int columns)
Функция Generate создает матрицу и заполняет ее случайными значениями. В качестве параметров передается количество строк и столбцов матрицы. Функция возвращает указатель на созданную структуру. Если создание матрицы не возможно, то возвращается NULL. Для физического размещения матрицы в памяти функция обращается к функции Create.
void Delete(Matrix ** M)
Функция Delete освобождает память заданную структурой Matrix. В качестве параметра функция принимает указатель на область память содержащий структуру. После выполнения функции, переданная матрица перестает существовать, указатель принимает значение NULL.
void Print(char * Head, Matrix * M, int width)
Функция Print используется для вывода матрицы на консоль и/или файл. Параметр Head содержит заголовок для матрицы, он используется для генерации имени текстового файла для записи данных матрицы. Кроме этого в функцию передается указатель на структуру содержащую матрицу и ширину одного столбца при выводе. В том случае, если суммарная ширина строки матрицы превосходит ширину стандартной консоли, выполняется только запись в файл.
void mulStandart(Matrix *A, Matrix *B, Matrix * Result, int size, bool Fast)
Функция mulStandart выполняет стандартное перемножение квадратных фрагментов матриц. На вход функции подаются два указателя на структуры содержащие матрицы для перемножения. Третий указатель используется для размещения результата. size – содержит размер фрагментов матриц. Логический флаг Fast содержит признак, использовать ли более быстрый доступ к элементам матрицы.
Matrix * MulStandart(Matrix *A, Matrix *B)
Функция MulStandart выполняет стандартное перемножение матриц. На вход функции подаются два указателя на структуры содержащие матрицы для перемножения. В качестве результата возвращается указатель на структуру содержащую результат перемножения матриц.
bool Equals(Matrix *A, Matrix *B, double epsilon)
Функция Equals выполняет сравнение двух матриц. На вход функции подаются да указателя на структуры содержащие матрицы для сравнения и точность сравнения. Необходимость точности обусловлена тем, что матрицы содержат вещественные значения, которые не могут сравниваться обычным способом из-за особенности представления в памяти компьютера. Функция возвращает истину, если матрицы идентичны и ложь в противном случае.
int expand(int k)
Функция expand определяет ближайшую степень двойки меньшую k (передается как параметр). Найденное значение возвращается в качестве результата работы.
void Copy(Matrix * From, Matrix * To, int size, int drfrom, int dcfrom, int drto, int dcto)
Функция Copy создает копию квадратной части матрицы в другой. В качестве параметров передается указатель на исходную матрицу, указатель на матрицу где необходимо разместить скопированный фрагмент, размер фрагмента, координаты начала копировки в исходной матрице (строка и столбец), координаты начала размещения в целевой матрице (строка и столбец).
void CopyTo(Matrix * From, Matrix * To, int drfrom, int dcfrom)
Функция CopyTo выполняет копирование части матрицы From в матрицу To. В функцию передаются строка и столбец до которых необходимо выполнить копирование.
void CopyFrom(Matrix * From, Matrix * To, int drto, int dcto)
Функция CopyFrom выполняет копирование матрицу From в часть матрицы To. В функцию передаются строка и столбец до которых необходимо выполнить копирование.
void Add(Matrix * M, Matrix * To, bool Fast)
Функция Add выполняет сложения двух матриц. Результат сохраняется в матрице, переданной через указатель To . На вход функции подаются два указателя на структуры содержащие матрицы для сложения. Логический флаг Fast содержит признак, использовать ли более быстрый доступ к элементам матрицы.
void Sub(Matrix * M, Matrix * From, bool Fast)
Функция Sub выполняет вычитание одной матрицы из другой. Результат сохраняется в матрице, переданной через указатель From . На вход функции подаются два указателя на структуры содержащие матрицы для вычитания. Логический флаг Fast содержит признак, использовать ли более быстрый доступ к элементам матрицы.
void mulStrassen(Matrix *A, Matrix *B, Matrix * Result, int size, bool Fast)
Функция mulStrassen выполняет перемножение матриц c фиктивным размером size алгоритмом Штрассена. Третий указатель используется для размещения результата. size – содержит размер фрагментов матриц. Логический флаг Fast содержит признак, использовать ли более быстрый доступ к элементам матрицы.
Matrix * MulStrassen(Matrix *A, Matrix *B, bool Fast)
Функция MulStrassen выполняет перемножение матриц классическим алгоритмом Штрассена. На вход функции подаются два указателя на структуры содержащие матрицы для перемножения. В качестве результата возвращается указатель на структуру содержащую результат перемножения матриц.
Matrix * MulMatrix(Matrix *A, Matrix *B)
Функция MulMatrix выполняет перемножение матриц модифицированным алгоритмом Штрассена. На вход функции подаются два указателя на структуры содержащие матрицы для перемножения. В качестве результата возвращается указатель на структуру содержащую результат перемножения матриц.
3. Результаты работы

Для запуска программы необходимо запустить на выполнение файл Shtrassen.exe.

Сразу после запуска программа запрашивает размерность матриц. Матрицы заполняются случайными числами и сохраняются в файлы. Если матрица не большого размера, то ее содержимое также выводится на экран.
После этого программа выполняет перемножение матриц двумя способами:
-классическим алгоритмом;
-модифицированным методом Штрассена.
Результаты сохраняются в текстовый файл.
Для каждого алгоритма выводится:
-максимальный объем задействованной оперативной памяти;
-количество выполненных умножений;
-затраченное время.
На рисунке 3.1 показан пример применения алгоритма для матриц небольшого размера. Использованы матрицы размерность 10 на 10.

Рисунок 3.1. Перемножение матриц небольшого размера
На рисунке 3.2 показан пример вычисления перемножения матриц среднего размера 100 на 100.

Рисунок 3.2. Перемножение матриц среднего размера
На рисунках 3.3-3.4 показан пример вычисления перемножения матриц большого размера 1000 на 1000. В данном случае процесс занимает заметное для пользователя время и компьютеру может не хватить памяти.

Рисунок 3.3. Процесс перемножение матриц большого размера

Рисунок 3.4. Нехватка памяти
Для сравнения эффективности алгоритмов проведем ряд вычислений для матриц различных размеров. Измерения проводились для квадратных матриц.
Затраты вычисляются по количеству умножений так как реальные затраты времени сильно зависят от способа реализации. Например, при динамическом выделении памяти основные расходы времени - на работу ОС по ее выделению/освобождению.
Результаты замеров внесем в Excel для последующей обработки (рисунок 3.5).

Рисунок 3.5. Результаты замеров
Коэффициент K это отношение затрат стандартного алгоритма к затратам модифицированного алгоритма Штрассена. Чем больше K, тем эффективнее работает смешанный алгоритм.
На рисунке 3.6 показана зависимость K от размера матрицы.

Рисунок 3.6. Зависимость K от размера матрицы
Видно, что в среднем с ростом размерности матриц выигрыш увеличивается, вероятно, по логарифмическому закону
На размерностях, равных степеням двойки резко улучшается производительность. Потом медленно уменьшается, достигая локального минимума на размерности 2^n-1.
Заключение

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

Разработанная программа позволяет сравнить между собой быстродействие алгоритма Штрассена и классического алгоритма перемножения матриц. Быстродействие сравнивается между собой по количеству выполненных операций умножения.
Полученные результаты сравнений были обработаны и получен график зависимости коэффициента выигрыша от размерности матрицы. Данный коэффициент растет с увеличением размерности матрицы. Максимальный выигрыш достигается на матрицах имеющих размерность равную степеням числа 2.
Для разработки программы использовался язык С++ и среда разработки Visual Studio.
Список использованной литературы

1. Ахо Дж., Хопкрофт А., Ульман Дж. Структуры данных и алгоритмы. — Москва.: Вильямс, 2016. - 400 c.

2. Кормен Томас Х. Алгоритмы: построение и анализ / Томас Х. Кормен, Чарльз И. Лейзерсон, Рональд Л. Ривест, Клиффорд Штайн – М.: Вильямс, 2016 - 1328с.
3. Левитин А. В. Алгоритмы: введение в разработку и анализ / Ананий Левитин, Игорь Красиков – М.: Вильямс, 2006 - 576с.
4. Макконнел Дж. Основы современных алгоритмов / Дж. Макконнел – М.: Техносфера, 2006 - 368с.
5. Пахомов Б. C/C++ и MS Visual C++ 2012 для начинающих / Борис Пахомов –СПб.: БХВ-Петербург, 2015 - 518с.
6. Хортон А. Visual C++ 2010 Полный курс/ А. Хортон – М.: Вильямс Издательский дом, 2012 - 1216с.
7. Шилдт Г. C++. Полное руководство / Герберт Шилдт – М.: Вильямс, 2017 - 704с.
Приложение А. Исходный код программы

#define _CRT_SECURE_NO_WARNINGS

#include
#include
#include
#include

//Отладочная печать

//#define DEBUG

//Удобные определения, которых нет в C

#define bool unsigned char
#define true 1
#define false 0

//Тип "матрица"

typedef struct Matrix
{
int rows; //строчек
int columns; //столбцов
int realrows; //Реально выделенных строк - для аккуратной работы с диспетчером памяти
double ** data; //данные
} Matrix;

//Методы получения и записи данных

double get(Matrix * M, int row, int col)
{
if (row >= M->rows) return 0;
if (col >= M->columns) return 0;
return M->data[row][col];
}

void set(Matrix * M, int row, int col, double value)

{
if (row >= M->rows) return;
if (col >= M->columns) return;
M->data[row][col] = value;
}

long long TotalRAM = 0;

long long MaxRAM = 0;
int Created = 0;
long MULS = 0;
long ADDS = 0;

void RegRAM(int delta)

{
TotalRAM += delta;
if (MaxRAM < TotalRAM) MaxRAM = TotalRAM;
}

Matrix * Create(int rows,int columns)

//создать матрицу заданного размера
{
Created++;

if (rows == 0)

{
//printf("0 rows %d\n",Created);
}
if (columns == 0)
{
//printf("0 columns %d\n",Created);
}

int r,c;
Matrix * Result = malloc(sizeof(Matrix));

RegRAM(sizeof(Matrix));
if (Result == NULL)
{
printf("Out of memory\n");
getchar();
return NULL;
}

Result->rows = rows;

Result->columns = columns;
//Простая защита от выделения 0 строк
Result->realrows = Result->rows;
if (Result->rows == 0) {
Result->data = malloc(1 * sizeof(double*));
Result->realrows = 1;
} else
Result->data = malloc(Result->rows * sizeof(double*));

if (Result->data == NULL)

{
printf("Out of memory\n");
getchar();
return NULL;
}
RegRAM(Result->realrows * sizeof(double*));
for (r = 0; r < Result->realrows; r++)
{
//Простая защита от нулевых запросов выделения памяти
if (Result->columns == 0)
{
Result->data[r] = malloc(sizeof(double));
RegRAM(sizeof(double));
}
else
{
RegRAM(Result->columns*sizeof(double));
Result->data[r] = malloc(Result->columns*sizeof(double));
}
if (Result->data[r] == NULL)
{
printf("Out of memory\n");
getchar();
return NULL;
}
}

for (r = 0; r < Result->rows; r++)

for (c = 0; c < Result->columns; c++)
Result->data[r][c] = 0; //На всякий случай

return Result;

}

Matrix * Generate(int rows, int columns)

//Создать случайную матрицу
{
int r, c;
Matrix * Result = Create(rows, columns);
for (r = 0; r < rows; r++)
for (c = 0; c < columns; c++)
#ifdef DEBUG
Result->data[r][c] = rand() % 10;
#else
Result->data[r][c] = rand() % 1000 / 100.0;
#endif // DEBUG
return Result;
}

void Delete(Matrix ** M)

//Освободить память
{
if ((*M) == NULL) return;
RegRAM(-(*M)->columns*(*M)->realrows*sizeof(double) + (*M)->rows*sizeof(double*) + sizeof(Matrix));
int r;
for (r = 0; r < (*M)->realrows; r++) free((*M)->data[r]);
free(*M);
*M = NULL;
}

void Print(char * Head,Matrix * M, int width)

//Вывести на консоль или в файл
{
int r, c;
char Template[10];
char FileName[80];
FILE * F;
if (!M)
{
printf("No data %s\n",Head);
return;
}
sprintf(Template, " %d.1f", width);
Template[0] = '%';
printf("%s [%dx%d]\n", Head, M->rows, M->columns);

if (M->columns * width > 79)

{
sprintf(FileName, "%s.txt", Head);
printf("Matrix %s is too large, saved to file \"%s\"\n", Head, FileName);
F = fopen(FileName, "w");
for (r = 0; r < M->rows; r++)
{
for (c = 0; c < M->columns; c++) fprintf(F,Template, M->data[r][c]);
fprintf(F,"\n");
}
fclose(F);
return;
}

for (r = 0; r < M->rows; r++)

{
for (c = 0; c < M->columns; c++) printf(Template, M->data[r][c]);
printf("\n");
}
}

#ifdef DEBUG

#define debugPrint(a,b,c) Print(a,b,c)
#else
#define debugPrint(a,b,c)
#endif

void mulStandart(Matrix *A, Matrix *B, Matrix * Result, int size, bool Fast)

//Умножение "стандартным" способом фрагментов матриц A&B размером size
//результат помещается в result
//dra итд - смещения внутри матриц
{
//Для отчета
int r, c, k;
double S, TA, TB;
for (r = 0; r < size; r++)
for (c = 0; c < size; c++)
{
S = 0;
for (k = 0; k < size; k++)
{
if (Fast)
S+=A->data[r][k]*B->data[k][c];
else {
TA = get(A, r, k);
TB = get(B, k, c);
S += TA*TB;
}
MULS++;
ADDS++;
}
if (Fast)
Result->data[r][c] = S;
else
set(Result, r, c, S);
}
}
Matrix * MulStandart(Matrix *A, Matrix *B)
//Стандартное умножение
{
int r, c, k;
double S;
Matrix * Result;
if (!A || !B)
{
printf("No data\n");
return NULL;
}
Result = Create(A->rows, B->columns);

//Защита от некорректного размера

if (A->columns != B->rows)
{
printf("Invalid matrix sizes\n");
Delete(&Result);
return NULL;
}

for (r = 0; r < Result->rows; r++)

for (c = 0; c < Result->columns; c++)
{
S = 0;
for (k = 0; k < A->columns; k++)
{
S += A->data[r][k] * B->data[k][c];
MULS++;
ADDS++;
}
Result->data[r][c] = S;
}
return Result;
}

bool Equals(Matrix *A, Matrix *B, double epsilon)

//Проверить совпадение двух матриц
{
int r, c;
if (!A) return false;
if (!B) return false;
if (epsilon < 0) epsilon = -epsilon;
if (epsilon == 0) epsilon = 0.0000001;
if (A->rows != B->rows) return false; //Не совпадают размеры
if (A->columns != B->columns) return false; //Не совпадают размеры
for (r = 0; r < A->rows; r++)
for (c = 0; c < A->columns; c++)
if (fabs(A->data[r][c] - B->data[r][c]) > epsilon) return false; //Большое различие!
return true; //Все элементы почти одинаковы
}

//**************************************************************************//

int expand(int k)

//Найти ближайшую не меньшую k cтепень двойки
{
if (k == 0) return 0;
int result = 1;
while (result < k) result <<= 1;
return result;
}

void Copy(Matrix * From, Matrix * To, int size, int drfrom, int dcfrom, int drto, int dcto)

//Создать копию части матрицы в другой
//Размер фиксирован, матрица псевдоквадратная
{
int r, c;
double T;
for (r = 0; r < size; r++)
for (c = 0; c < size; c++)
{
T = get(From, r + drfrom, c + dcfrom);
set(To, r+drto, c+dcto, T);
}
}

void CopyTo(Matrix * From, Matrix * To, int drfrom, int dcfrom)

//Скопировать часть матрицы From в матрицу To
{
int r, c;
double T;
for (r = 0; r < To->rows; r++)
for (c = 0; c < To->columns; c++)
{
T = From->data[r+drfrom][c + dcfrom];
To->data[r][c] = T;
}
}

void CopyFrom(Matrix * From, Matrix * To, int drto, int dcto)

//Скопировать матрицу From в часть матрицы TO
{
int r, c;
double T;
for (r = 0; r < From->rows; r++)
for (c = 0; c < From->columns; c++)
{
T = From->data[r][c];
To->data[r + drto][c + dcto] = T;
}
}
void Add(Matrix * M, Matrix * To, bool Fast)
{
int r, c;
for (r = 0; r < M->rows; r++)
for (c = 0; c < M->columns; c++)
if (Fast)
To->data[r][c] += M->data[r][c];
else
set(To,r,c,get(To,r,c)+ get(M,r,c));
}

void Sub(Matrix * M, Matrix * From, bool Fast)

{
int r, c;
for (r = 0; r < M->rows; r++)
for (c = 0; c < M->columns; c++)
if (Fast)
From->data[r][c] -= M->data[r][c];
else
set(From, r, c, get(From, r, c) - get(M, r, c));
}

void mulStrassen(Matrix *A, Matrix *B, Matrix * Result, int size, bool Fast)

//Перемножить методом Штрассена матрицы с фиктивным размером size
//Функция рекурсивная
//Проведем разбиение на "подматрицы"
// A B
// A11 A12 B11 B12
// A21 A22 B21 B22
{
//При некотором size, меньшем заранее заданной величины, используется стандартный алгоритм умножения
if (size < 4) //Не столько! Гораздо больше. 4 - в демонстрационных целях
{
mulStandart(A, B, Result, size, Fast);
return;
}
//Иначе производится серия вычислений, описанных Штрассеном
//Создается (это неизбежно) несколько матриц P
//P1 = (A11+A22)*(B11+B22)
//P2 = (A21+A22) * B11
//P3 = A11 * (B12 - B22)
//P4 = A22 * (B21 - B11)
//P5 = (A11 + A12) * B22
//P6 = (A21 - A11) * (B11+B12)
//P7 = (A12 - A22) * (B21+B22)

size /= 2;

//Выделить четвертушки

//P1 = (A11+A22)*(B11+B22)

Matrix * A11 = Create(size, size); Copy(A, A11, size, 0, 0, 0, 0);
Matrix * A22 = Create(size, size); Copy(A, A22, size, size, size, 0, 0);
Matrix * T1 = Create(size, size); Copy(A11, T1, size, 0, 0, 0, 0);
Add(A22, T1, Fast);

Matrix * A12 = Create(size, size); Copy(A, A12, size, 0, size, 0, 0);

Matrix * A21 = Create(size, size); Copy(A, A21, size, size, 0, 0, 0);
Matrix * B11 = Create(size, size); Copy(B, B11, size, 0, 0, 0, 0);
Matrix * B12 = Create(size, size); Copy(B, B12, size, 0, size, 0, 0);
Matrix * B21 = Create(size, size); Copy(B, B21, size, size, 0, 0, 0);
Matrix * B22 = Create(size, size); Copy(B, B22, size, size, size, 0, 0);
Matrix * T2 = Create(size, size); Copy(B11, T2, size, 0, 0, 0, 0);
Add(B22, T2, Fast);
Matrix * P1 = Create(size, size);
mulStrassen(T1, T2, P1, size, true);
//P2 = (A21+A22) * B11
Matrix * P2 = Create(size, size);
Copy(A21, T1, size, 0, 0, 0, 0);
Add(A22, T1, Fast);
mulStrassen(T1, B11, P2, size, true);

//P3 = A11 * (B12 - B22)

Matrix * P3 = Create(size, size);
Copy(B12, T1, size, 0, 0, 0, 0);
Sub(B22, T1, Fast);
mulStrassen(A11, T1, P3, size, true);
//P4 = A22 * (B21 - B11)
Matrix * P4 = Create(size, size);
Copy(B21, T1, size, 0, 0, 0, 0);
Sub(B11, T1, Fast);
mulStrassen(A22, T1, P4, size, true);

//P5 = (A11 + A12) * B22

Matrix * P5 = Create(size, size);
Copy(A11, T1, size, 0, 0, 0, 0);
Add(A12, T1, Fast);
mulStrassen(T1, B22, P5, size, true);
//P6 = (A21 - A11) * (B11+B12)
Matrix * P6 = Create(size, size);
Copy(A21, T1, size, 0, 0, 0, 0);
Sub(A11, T1, Fast);
Copy(B11, T2, size, 0, 0, 0, 0);
Delete(&B11);
Add(B12, T2, Fast);
Delete(&A11);
Delete(&A21);
Delete(&B12);
mulStrassen(T1, T2, P6, size, true);

//P7 = (A12 - A22) * (B21+B22)

Matrix * P7 = Create(size, size);
Copy(A12, T1, size, 0, 0, 0, 0);
Sub(A22, T1, Fast);
Copy(B21, T2, size, 0, 0, 0, 0);
Add(B22, T2, Fast);
Delete(&A12);
Delete(&A22);
Delete(&B21);
Delete(&B22);
mulStrassen(T1, T2, P7, size, true);
Delete(&T2);

//Тогда результирующая матрица будет образована так:

//C11 = P1 + P4 - P5 + P7
Copy(P1,T1,size,0,0,0,0);
Add(P4, T1, Fast);
Sub(P5, T1, Fast);
Add(P7, T1, Fast);
Delete(&P7);
Copy(T1, Result, size, 0, 0, 0, 0);
Delete(&T1);

//C12 = P3 + P5

Add(P3,P5, Fast); //значение P5 больше не нужно
Copy(P5, Result, size, 0, 0, 0, size);
Delete(&P5);
//C21 = P2 + P4
Add(P2,P4, Fast); //P4 не нужен
Copy(P4, Result, size, 0, 0, size, 0);
Delete(&P4);

//C22 = P1 - P2 +P3 + P6

Sub(P2, P1, Fast);
Delete(&P2);
Add(P3, P1, Fast);
Delete(&P3);
Add(P6, P1, Fast);
Delete(&P6);
Copy(P1,Result,size,0,0,size,size);
Delete(&P1);

//Освободить память

}

Matrix * MulStrassen(Matrix *A, Matrix *B, bool Fast)

{
int ra = expand(A->rows);
int ca = expand(A->columns);
int rb = expand(B->rows);
int cb = expand(B->columns);
int size = max(max(ra, rb), max(ca, cb));
Matrix * Result = Create(A->rows, B->columns);
mulStrassen(A, B, Result, size, Fast);
return Result;
}

Matrix * MulMatrix(Matrix *A, Matrix *B)

{
//Перемножение матриц смешанным способом.
//Выделяется под матрица, годная для перемножения алгоритмом
//Штрассена, и применяется формула
//C1 C2
//C3 C4
//C1 = A1*B1 + A2*B3
//C2 = A2*B2 + A2*B1
//C3 = A3*B1 + A4*B3
//C4 = A3*B2 + A4*B4

//Определить размеры

int sizear = expand(A->rows); //Первой матрицы - строк
int sizeac = expand(A->columns);
int sizea = min(sizear, sizeac);
int sizebr = expand(B->rows);
int sizebc = expand(B->columns);
int sizeb = min(sizebr, sizebc); //Копипаст - вреден для здоровья!
int size = min(sizea, sizeb);

#ifdef DEBUG

printf("Mul.size = %d\n", size);
#endif // DEBUG

//Возможны случаи:

//1. Обе матрицы квадратные и их размеры годные для Штрассена.
//в этом случае все просто (и встречается не так редко, как может показаться)
if (sizear == A->rows && sizeac == A->columns && sizear == sizeac
&& sizebr == B->rows && sizebc == B->columns && sizebr == sizebc
&& sizear == sizebr
)
return MulStrassen(A, B, true);
//2. Матрицы маленькие. В этом случае лучше использовать стандартный способ умножения
const int Small = 2; //На самом деле - примерно 32
if (sizear <= Small || sizeac <= Small || sizebr <= Small || sizebc <= Small)
return MulStandart(A, B);
//Выбрать безопасный size
while (size>A->rows) size /= 2;
while (size>A->columns) size /= 2;
while (size>B->rows) size /= 2;
while (size>B->columns) size /= 2;
//Свести к другим методам не получилось
//Выделить матрицы
//Из первой
Matrix * A1 = Create(size, size);
CopyTo(A, A1, 0, 0);
debugPrint("A1", A1, 5);
Matrix * B1 = Create(size, size);
CopyTo(B, B1, 0, 0);
debugPrint("B1", B1, 5);
Matrix * C1 = MulMatrix(A1, B1);

Matrix * A2 = Create(size, A->columns - size);

CopyTo(A, A2, 0, size);
debugPrint("A2", A2, 5);

Matrix * B3 = Create(B->rows - size, size);

CopyTo(B, B3, size, 0);
debugPrint("B3", B3, 5);

Matrix * T1 = MulMatrix(A2, B3);

Add(T1, C1, true);
debugPrint("C1", C1, 7);
Delete(&T1);

Matrix * A3 = Create(A->rows - size, size);

CopyTo(A, A3, size, 0);
debugPrint("A3", A3, 5);

Matrix * A4 = Create(A->rows - size, A->columns - size);

CopyTo(A, A4, size, size);
debugPrint("A4", A4, 5);

Matrix * B2 = Create(size, B->columns - size);

CopyTo(B, B2, 0, size);
debugPrint("B2", B2, 5);

Matrix * B4 = Create(B->rows - size, B->columns - size);

CopyTo(B, B4, size, size);
debugPrint("B4", B4, 5);

Matrix * C2 = MulMatrix(A1, B2);

Delete(&A1);

Matrix * T2 = MulMatrix(A2, B4);

Delete(&A2);

Add(T2, C2, true);

debugPrint("C2", C2, 7);
Delete(&T2);

//Произвести вычисления

Matrix * C3 = MulMatrix(A3, B1);
Delete(&B1);
Matrix * T3 = MulMatrix(A4, B3);
Delete(&B3);

Add(T3, C3, true);

debugPrint("C3", C3, 7);
Delete(&T3);

Matrix * C4 = MulMatrix(A3, B2);

Delete(&A3);
Delete(&B2);

Matrix * T4 = MulMatrix(A4, B4);

Delete(&A4);
Delete(&B4);

Add(T4, C4, true);

debugPrint("C4", C4, 7);
Delete(&T4);

//Собираем матрицу C

Matrix * C = Create(A->rows, B->columns);
CopyFrom(C1, C, 0, 0);
debugPrint("C(C1)", C, 7);
CopyFrom(C2, C, 0,size);//Уточнить размер
debugPrint("C(C2)", C, 7);
CopyFrom(C3, C, size,0);//Уточнить размер
debugPrint("C(C3)", C, 7);
CopyFrom(C4, C, size,size);//Уточнить размер
debugPrint("C(C4)", C, 7);
Delete(&C1); Delete(&C2); Delete(&C3); Delete(&C4);
return C;
}
//**************************************************************************//

int main()

{
srand((unsigned)time(0));
int Arows=50,Acols=80,Bcols=50;
printf("A.rows A.cols B.cols= ");
scanf("%d %d %d", &Arows, &Acols, &Bcols); getchar();

if (Arows < 1) Arows = 1;

if (Acols < 1) Acols = 1;
if (Bcols < 1) Bcols = 1;
Matrix * A = Generate(Arows,Acols);
Matrix * B = Generate(Acols, Bcols); //A.cols==B.rows!
Print("A", A, 5);
Print("B", B, 5);

//Matrix * C = MulStrassen(A, B, false);

//Print("AB-Shtrassen",C,10);
//if (!Equals(C, D, 0.00001)) printf("ERROR Shtrassen\n");

clock_t T1 = clock();

Matrix * D = MulStandart(A, B);
int MULSS = MULS;
int ADDSS = ADDS;
MULS = 0;
ADDS = 0;
clock_t T2 = clock();
Matrix * Z = MulMatrix(A, B);
clock_t T3 = clock();
Print("AB-Standart",D, 10);
Print("AB-Mix", Z, 10);
printf("\n");
if (!Equals(D, Z, 0.00001)) printf("ERROR Mix\n");
else
printf("OK\n");

if (MaxRAM<1024) printf("Max RAM = %lld bytes\n", MaxRAM);

else if (MaxRAM<1024*1024) printf("Max RAM = %lld Kb\n", MaxRAM/1024);
else if (MaxRAM<1024 * 1024 * 1024) printf("Max RAM = %lld Mb\n", MaxRAM / 1024/1024);
else printf("Max RAM = %lld Gb\n", MaxRAM / 1024 / 1024 / 1024);

printf("Standart MULS %d\n",MULSS);

// printf("Standart ADDS %d\n", ADDSS);
printf("Mix MULS %d\n", MULS);
// printf("Mix ADDS %d\n", ADDS);

int dT1 = (T2 - T1) * 1000 / CLK_TCK;

int dT2 = (T3 - T2) * 1000 / CLK_TCK;

printf("Standart Time = %d ms\n", dT1);

printf("Mix Time = %d ms\n", dT2);

getchar();

}