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 64 65 66 67 68 69 70

#include <complex>
#include <iostream>
#include <iomanip>
#include <fftw3.h>
int xy(int x, int y, int number_of_cols) { return x + y * number_of_cols; }
void matrix_product(
fftw_complex* out,
fftw_complex const* a, int a_width, int a_height,
double const* b, int b_width, int b_height)
{
for (int x = 0; x < b_width; ++x)
for (int y = 0; y < a_height; ++y)
{
// compute inner product of xth column of b with yth row of a
fftw_complex inner_product {0.0, 0.0};
for (int i = 0; i < a_width; ++i)
{
double const p_real = a[xy(i, y, a_width)][0];
double const p_imag = a[xy(i, y, a_width)][1];
double const q = b[xy(x, i, b_width)];
inner_product[0] += p_real * q;
inner_product[1] += p_imag * q;
}
out[xy(x, y, b_width)][0] = inner_product[0];
out[xy(x, y, b_width)][1] = inner_product[1];
}
}
void print(fftw_complex* m, int width, int height)
{
std::cout << std::fixed << std::setprecision(2);
for (int y = 0; y < height; ++y)
{
for (int x = 0; x < width; ++x)
std::cout << std::setw(6) << std::complex{m[xy(x, y, width)][0], m[xy(x, y, width)][1]} << ' ';
std::cout << '\n';
}
std::cout << '\n';
}
int main()
{
fftw_complex a[] {
{1.0, 0.0}, {2.0, 0.0}, {3.0, 0.0},
{4.0, 0.0}, {5.0, 0.0}, {6.0, 0.0}
};
int const a_width = 3;
int const a_height = 2;
double b[] {
7.0, 8.0,
9.0, 10.0,
11.0, 12.0
};
int const b_width = 2;
int const b_height = 3;
int const p_width = b_width;
int const p_height = a_height;
fftw_complex p[p_width * p_height];
matrix_product(p, a, a_width, a_height, b, b_width, b_height);
print(p, p_width, p_height);
}
