-
Notifications
You must be signed in to change notification settings - Fork 57
/
matmul_v1.c
64 lines (61 loc) · 1.29 KB
/
matmul_v1.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
// Writen by Attractive Chaos; distributed under the MIT license
#include <stdlib.h>
#include <stdio.h>
double **mm_init(int n)
{
double **m;
int i;
m = (double**)malloc(n * sizeof(void*));
for (i = 0; i < n; ++i)
m[i] = calloc(n, sizeof(double));
return m;
}
void mm_destroy(int n, double **m)
{
int i;
for (i = 0; i < n; ++i) free(m[i]);
free(m);
}
double **mm_gen(int n)
{
double **m, tmp = 1. / n / n;
int i, j;
m = mm_init(n);
for (i = 0; i < n; ++i)
for (j = 0; j < n; ++j)
m[i][j] = tmp * (i - j) * (i + j);
return m;
}
// better cache performance by transposing the second matrix
double **mm_mul(int n, double *const *a, double *const *b)
{
int i, j, k;
double **m, **c;
m = mm_init(n); c = mm_init(n);
for (i = 0; i < n; ++i) // transpose
for (j = 0; j < n; ++j)
c[i][j] = b[j][i];
for (i = 0; i < n; ++i) {
double *p = a[i], *q = m[i];
for (j = 0; j < n; ++j) {
double t = 0.0, *r = c[j];
for (k = 0; k < n; ++k)
t += p[k] * r[k];
q[j] = t;
}
}
mm_destroy(n, c);
return m;
}
int main(int argc, char *argv[])
{
int n = 100;
double **a, **b, **m;
if (argc > 1) n = atoi(argv[1]);
n = (n/2) * 2;
a = mm_gen(n); b = mm_gen(n);
m = mm_mul(n, a, b);
fprintf(stderr, "%lf\n", m[n/2][n/2]);
mm_destroy(n, a); mm_destroy(n, b); mm_destroy(n, m);
return 0;
}