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
|
#include <stdio.h>
typedef struct {
double re, im;
} COMPLEX;
static inline void COMPLEX_zero(COMPLEX *r) {
r->re = r->im = 0.0;
}
static inline void COMPLEX_add(COMPLEX *r,
const COMPLEX *x,
const COMPLEX *y) {
double re = x->re + y->re;
double im = x->im + y->im;
r->re = re;
r->im = im;
}
static inline void COMPLEX_mul(COMPLEX *r,
const COMPLEX *x,
const COMPLEX *y) {
double re = x->re * y->re - x->im * y->im;
double im = x->im * y->re + x->re * y->im;
r->re = re;
r->im = im;
}
static inline void COMPLEX_mul_add(COMPLEX *r,
const COMPLEX *x,
const COMPLEX *y) {
double re = r->re + x->re * y->re - x->im * y->im;
double im = r->im + x->im * y->re + x->re * y->im;
r->re = re;
r->im = im;
}
void COMPLEX_mat_mul1(unsigned m, unsigned n, unsigned p,
COMPLEX * restrict c, unsigned stride_c,
const COMPLEX *a, unsigned stride_a,
const COMPLEX *b, unsigned stride_b) {
for(unsigned i=0; i<m; i++) {
for(unsigned k=0; k<p; k++) {
COMPLEX_zero(c+i*stride_c+k);
}
}
for(unsigned i=0; i<m; i++) {
for(unsigned k=0; k<p; k++) {
for(unsigned j=0; j<n; j++) {
COMPLEX_mul_add(c+i*stride_c+k, a+i*stride_a+j, b+j*stride_b+k);
}
}
}
}
int main() {
COMPLEX a = { 1, 2 };
COMPLEX b = { 7, 4 };
COMPLEX r;
COMPLEX_mul(&r, &a, &b);
printf("%g %g\n", r.re, r.im);
}
|
