Как использовать этот код C для умножения двух матриц с использованием алгоритма Штрассена?
Я искал реализацию алгоритма Штрассена в C, и я нашел этот код в конце.
Чтобы использовать multiply функция:
void multiply(int n, matrix a, matrix b, matrix c, matrix d);
который умножает две матрицы a, b и помещает результат в c (d является промежуточной матрицей). Матрицы a а также b должен иметь следующий тип:
typedef union _matrix
{
double **d;
union _matrix **p;
} *matrix;
Я выделил динамически четыре матрицы a, b, c, d (двумерные массивы двойников) и присвоили их адреса в поле _matrix.d:
#include "strassen.h"
#define SIZE 50
int main(int argc, char *argv[])
{
double ** matA, ** matB, ** matC, ** matD;
union _matrix ma, mb, mc, md;
int i = 0, j = 0, n;
matA = (double **) malloc(sizeof(double *) * SIZE);
for (i = 0; i < SIZE; i++)
matA[i] = (double *) malloc(sizeof(double) * SIZE);
// Do the same for matB, matC, matD.
ma.d = matA;
mb.d = matB;
mc.d = matC;
md.d = matD;
// Initialize matC and matD to 0.
// Read n.
// Read matA and matB.
multiply(n, &ma, &mb, &mc, &md);
return 0;
}
Этот код успешно компилируется, но завершается с n > BREAK,
strassen.c:
#include "strassen.h"
/* c = a * b */
void multiply(int n, matrix a, matrix b, matrix c, matrix d)
{
if (n <= BREAK) {
double sum, **p = a->d, **q = b->d, **r = c->d;
int i, j, k;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++) {
for (sum = 0., k = 0; k < n; k++)
sum += p[i][k] * q[k][j];
r[i][j] = sum;
}
} else {
n /= 2;
sub(n, a12, a22, d11);
add(n, b21, b22, d12);
multiply(n, d11, d12, c11, d21);
sub(n, a21, a11, d11);
add(n, b11, b12, d12);
multiply(n, d11, d12, c22, d21);
add(n, a11, a12, d11);
multiply(n, d11, b22, c12, d12);
sub(n, c11, c12, c11);
sub(n, b21, b11, d11);
multiply(n, a22, d11, c21, d12);
add(n, c21, c11, c11);
sub(n, b12, b22, d11);
multiply(n, a11, d11, d12, d21);
add(n, d12, c12, c12);
add(n, d12, c22, c22);
add(n, a21, a22, d11);
multiply(n, d11, b11, d12, d21);
add(n, d12, c21, c21);
sub(n, c22, d12, c22);
add(n, a11, a22, d11);
add(n, b11, b22, d12);
multiply(n, d11, d12, d21, d22);
add(n, d21, c11, c11);
add(n, d21, c22, c22);
}
}
/* c = a + b */
void add(int n, matrix a, matrix b, matrix c)
{
if (n <= BREAK) {
double **p = a->d, **q = b->d, **r = c->d;
int i, j;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
r[i][j] = p[i][j] + q[i][j];
} else {
n /= 2;
add(n, a11, b11, c11);
add(n, a12, b12, c12);
add(n, a21, b21, c21);
add(n, a22, b22, c22);
}
}
/* c = a - b */
void sub(int n, matrix a, matrix b, matrix c)
{
if (n <= BREAK) {
double **p = a->d, **q = b->d, **r = c->d;
int i, j;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
r[i][j] = p[i][j] - q[i][j];
} else {
n /= 2;
sub(n, a11, b11, c11);
sub(n, a12, b12, c12);
sub(n, a21, b21, c21);
sub(n, a22, b22, c22);
}
}
strassen.h:
#define BREAK 8
typedef union _matrix {
double **d;
union _matrix **p;
} *matrix;
/* Notational shorthand to access submatrices for matrices named a, b, c, d */
#define a11 a->p[0]
#define a12 a->p[1]
#define a21 a->p[2]
#define a22 a->p[3]
#define b11 b->p[0]
#define b12 b->p[1]
#define b21 b->p[2]
#define b22 b->p[3]
#define c11 c->p[0]
#define c12 c->p[1]
#define c21 c->p[2]
#define c22 c->p[3]
#define d11 d->p[0]
#define d12 d->p[1]
#define d21 d->p[2]
#define d22 d->p[3]
Мой вопрос, как использовать функцию multiply (как реализовать матрицу).
4 ответа
Как сказал Атом, нужно правильно инициализировать matrix.p для обеих матриц.
1) Прежде всего, ваш matrix это союз так p по сути становится d интерпретируется как _matrix ** что здесь не имеет смысла - вот почему он падает. Вам, вероятно, нужно сделать matrix struct вместо.
В заключение, p по определению массив подматриц, поэтому он должен быть struct _matrix * (и вам нужно malloc фактический массив при необходимости) или struct _matrix[4] (что невозможно:)).
typedef struct _matrix
{
double **d;
struct _matrix *p;
} *matrix;
2) Теперь посмотрим что p должно быть.
│
A.d -> d1 -> a[1,1] a[1,2]│a[1,3] a[1,4]
d2 -> a[2,1] a[2,2]│a[2,3] a[2,4]
─────────────────────────────
d3 -> a[3,1] a[3,2]│a[3,3] a[3,4]
d4 -> a[4,1] a[4,2]│a[4,3] a[4,4]
│
p указывает на массив matrix структуры! Особенность состоит в том, чтобы сделать dиз этих структур указывают внутрь А таким образом, что (p[k].d)[i][j] является элементом соответствующей подматрицы:
p[0].d -> p01 -> a[1,1] p[1].d -> p11 -> a[1,3]
p02 -> a[2,1] p12 -> a[2,3]
p[2].d -> p21 -> a[3,1] p[3].d -> p31 -> a[3,3]
p22 -> a[4,1] p32 -> a[4,3]
Можете ли вы теперь вывести алгоритм инициализации p для квадрата А произвольного четного размера?
И КОГДА инициализировать его в первую очередь?;)
Когда n > BREAK, алгоритм умножения матриц использует иерархическое матричное представление (поле p из union _matrix, а не поле d).
Вам необходимо настроить свой код для иерархического представления при распределении памяти и при инициализации матриц a а также b,
Ну, я не совсем уверен, что не так в коде, который вы опубликовали, но я хотел бы указать вам на эту реализацию алгоритма Штрассена, который был использован в конкурсе программирования Intel некоторое время назад: http://software.intel.com/file/21818. Может быть, вы найдете там несколько полезных советов.